Анализ временных рядов с помощью python
В сегодняшней статье, я попытаюсь описать процесс анализа временных рядов с помощью python и модуля statsmodels . Данный модуль предоставляет широкий набор средств и методов для проведения статистического анализа и эконометрики.
Для примера взяты реальные данные по товарообороту одного из складских комплексов Подмосковья. Приблизительный список шагов для анализа можно посмотреть здесь.
Загрузка и предварительная обработка данных
Для начала загрузим данные и посмотрим на них:
from pandas import read_csv, DataFrame import statsmodels.api as sm from statsmodels.iolib.table import SimpleTable from sklearn.metrics import r2_score import ml_metrics as metrics dataset = read_csv('tovar_moving.csv',';', index_col=['date_oper'], parse_dates=['date_oper'], dayfirst=True) dataset.head() Итак, как можно заметить функция read_csv() , в данном случем помимо указания параметров, которые задают используемые колонки и индекс, можно заметить еще 3 параметра для работы с датой. Остановимся на них поподробнее.
parse_dates задает имена столбцов, которые будут преобразованы в тип DateTime . Стоит отметить, что если в данном столбце будут пустые значения парсинг не удастся и вернется столбец типа object . Чтобы этого избежать надо добавить параметр keep_default_na =False .
Заключительный параметр dayfirst указывает функции парсинга, что первое в строке первым идет день, а не наоборот. Если не задать этот параметр, то функция может не правильно преобразовывать даты и путать месяц и день местами. Например 01.02.2013 будет преобразовано в 02-01-2013 , что будет неправильно.
Выделим в отдельную серию временной ряд со значениями отгрузкок:
otg = dataset.Otgruzka otg.head() date_oper 2009-09-01 179667 2009-09-02 177670 2009-09-03 152112 2009-09-04 142938 2009-09-05 130741 Name: Otgruzka, dtype: int64 Итак у нас теперь есть временной ряд и можно перейти к его анализу.
Анализ временного ряда
Для начала давайте посмортим график нашего ряда:
Из графика видно, что наш ряд имеет небольшое кол-во выбросов, которые влияют на разброс. Кроме того анализировать отгрузки за каждый день не совсем верно, т.к., например, в конце или начале недели будут дни в которые товара отгружается значительно больше, нежели в остальные. Поэтому есть смысл перейти к месячному интервалу и среднему значению отгрузок на нем, это избавит нас от выбросов и уменьшит колебания нашего ряда. В pandas для этого есть удобная функция resample() , в качестве параметров ей передается период округления и аггрегатная функция:
otg = otg.resample('W', how='mean') otg.plot(figsize=(12,6)) 
Как можно заметить, новый график не имеет ярких выбросов и имеет ярко выраженный тренд. Из это можно сделать вывод о том, что ряд не является стационарным [1] .
Давайте посмотрим на характеристики данного ряда:
itog = otg.describe() otg.hist() itog count 225.000000 mean 270858.285365 std 118371.082975 min 872.857143 25% 180263.428571 50% 277898.714286 75% 355587.285714 max 552485.142857 dtype: float64 
Как можно заметить из характеристик и гистограммы, ряд у нас более менее однородный и имеет относительно небольшой разброс о чем свидетельствует коэффициент вариации : $V = \frac > $, где $\sigma$ — cреднеквадратическое отклонение , $\bar$ — среднее арифметическое выборки. В нашем случае он равен:
print 'V = %f' % (itog['std']/itog['mean']) Проведем тест Харки — Бера для определения номарльности распределения, чтобы подтвердить предположение об однородности. Для этого в существует функция jarque_bera() , которая возвращает значения данной статистики:
row = [u'JB', u'p-value', u'skew', u'kurtosis'] jb_test = sm.stats.stattools.jarque_bera(otg) a = np.vstack([jb_test]) itog = SimpleTable(a, row) print itog ========================================================== JB p-value skew kurtosis ---------------------------------------------------------- 5.60453241944 0.0606724103035 0.111910758759 2.25991843713 ---------------------------------------------------------- Значение данной статистика свидетельствует о том, нулевая гипотеза о нормальности распределения отвергается с малой вероятностью ( probably > 0.05 ), и, следовательно, наш ряд имеет нормального распределения.
Функция SimpleTable() служит для оформления вывода. В нашем случае на вход ей подается массив значений (размерность не больше 2) и список с названиями столбцов или строк.
Многие методы и модели основаны на предположениях о стационарности ряда, но как было замечено ранее наш ряд таковым скорее всего не является. Поэтому для проверки проверки стационарности давайте проведем обобщенный тест Дикки-Фуллера на наличие единичных корней. Для этого в модуле statsmodels есть функция adfuller() :
test = sm.tsa.adfuller(otg) print 'adf: ', test[0] print 'p-value: ', test[1] print'Critical values: ', test[4] if test[0]> test[4]['5%']: print 'есть единичные корни, ряд не стационарен' else: print 'единичных корней нет, ряд стационарен' adf: -1.38835541357 p-value: 0.58784577297 Critical values: есть единичные корни, ряд не стационарен Проведенный тест подтвердил предположения о не стационарности ряда. Во многих случаях взятие разности рядов позволяет это сделать.
Если, например, первые разности ряда стационарны, то он называется интегрированным рядом первого порядка .
Итак, давайте определим порядок интегрированного ряда для нашего ряда:
otg1diff = otg.diff(periods=1).dropna() В коде выше функция diff() вычисляет разность исходного ряда с рядом с заданным смещением периода. Период смещения передается как параметр period . Т.к. в разности первое значение получиться неопределенным, то нам надо избавиться от него для этого и используется метод dropna() .
Проверим получившийся ряд на стационарность:
test = sm.tsa.adfuller(otg1diff) print 'adf: ', test[0] print 'p-value: ', test[1] print'Critical values: ', test[4] if test[0]> test[4]['5%']: print 'есть единичные корни, ряд не стационарен' else: print 'единичных корней нет, ряд стационарен' adf: -5.95204224907 p-value: 2.13583392404e-07 Critical values: единичных корней нет, ряд стационарен Как видно из кода выше получившийся ряд первых разностей приблизился к стационарному. Для полной уверенности разобъем его на несколько промежутков и убедимся мат. ожидания на разных интервалах:
m = otg1diff.index[len(otg1diff.index)/2+1] r1 = sm.stats.DescrStatsW(otg1diff[m:]) r2 = sm.stats.DescrStatsW(otg1diff[:m]) print 'p-value: ', sm.stats.CompareMeans(r1,r2).ttest_ind()[1] Высокое p-value дает нам возможность утверждать, что нулевая гипотеза о равенстве средних верна, что свидетельствует о стационарности ряда. Осталось убедиться в отсутствии тренда для этого построим график нашего нового ряда:
Тренд действительно отсутствует, таким образом ряд первых разностей является стационарным, а наш исходный ряд — интегрированным рядом первого порядка .
Построение модели временного ряда
Для моделирования будем использовать модель ARIMA , построенную для ряда первых разностей.
Итак, чтобы построить модель нам нужно знать ее порядок, состоящий из 2-х параметров:
Параметр d есть и он равет 1, осталось определить p и q . Для их определения нам надо изучить авторкорреляционную(ACF) и частично автокорреляционную(PACF) функции для ряда первых разностей.
ACF поможет нам определить q , т. к. по ее коррелограмме можно определить количество автокорреляционных коэффициентов сильно отличных от 0 в модели MA
PACF поможет нам определить p , т. к. по ее коррелограмме можно определить максимальный номер коэффициента сильно отличный от 0 в модели AR .
Чтобы построить соответствующие коррелограммы, в пакете statsmodels имеются следующие функции: plot_acf() и plot_pacf() . Они выводят графики ACF и PACF , у которых по оси X откладываются номера лагов, а по оси Y значения соответствующих функций. Нужно отметить, что количество лагов в функциях и определяет число значимых коэффициентов. Итак, наши функции выглядят так:
fig = plt.figure(figsize=(12,8)) ax1 = fig.add_subplot(211) fig = sm.graphics.tsa.plot_acf(otg1diff.values.squeeze(), lags=25, ax=ax1) ax2 = fig.add_subplot(212) fig = sm.graphics.tsa.plot_pacf(otg1diff, lags=25, ax=ax2) 
После изучения коррелограммы PACF можно сделать вывод, что p = 1 , т.к. на ней только 1 лаг сильно отличнен от нуля. По коррелограмме ACF можно увидеть, что q = 3 , т.к. после лага 3 значении функций резко падают.
Итак, когда известны все параметры можно построить модель, но для ее построения мы возмем не все данные, а только часть. Данные из части не попавших в модель мы оставим для проверки точности прогноза нашей модели:
src_data_model = otg[:'2013-05-26'] model = sm.tsa.ARIMA(src_data_model, order=(1,1,1), freq='W').fit(full_output=False, disp=0) Параметр trend отвечает за наличие константы в моделе. Выведем информацаю по получившейся модели:
ARIMA Model Results ============================================================================== Dep. Variable: D.y No. Observations: 194 Model: ARIMA(1, 1, 1) Log Likelihood -2326.028 Method: css-mle S.D. of innovations 38615.075 Date: Tue, 24 Dec 2013 AIC 4660.057 Time: 02:12:47 BIC 4673.128 Sample: 09-13-2009 HQIC 4665.350 - 05-26-2013 ============================================================================== coef std err z P>|z| [95.0% Conf. Int.] ------------------------------------------------------------------------------ const 1588.2266 142.728 11.128 0.000 1308.484 1867.969 ar.L1.D.y 0.6660 0.055 12.166 0.000 0.559 0.773 ma.L1.D.y -1.0000 0.014 -72.214 0.000 -1.027 -0.973 Roots ============================================================================= Real Imaginary Modulus Frequency ----------------------------------------------------------------------------- AR.1 1.5015 +0.0000j 1.5015 0.0000 MA.1 1.0000 +0.0000j 1.0000 0.0000 ----------------------------------------------------------------------------- Как видно из данной информации в нашей модели все коэффициенты значимые и можно перейти к оценке модели.
Анализ и оценка модели
Проверим остатки данной модели на соответствие “белому шуму” , а также проанализируем коррелограму остатков, так как это может нам помочь в определении важных для включения и прогнозирования элементов регрессии.
Итак первое, что мы сделаем это проведем Q-тест Льюнга — Бокса для проверки гипотезы о том, что остатки случайны, т. е. являются “белым шумом”. Данный тест проводится на остатках модели ARIMA . Таким образом, нам надо сначала получить остатки модели и построить для них ACF, а затем к получившимся коэффициентам приметить тест. С помощью statsmadels это можно сделать так:
q_test = sm.tsa.stattools.acf(model.resid, qstat=True) #свойство resid, хранит остатки модели #qstat=True, означает что применяем указынный тест к коэф-ам print DataFrame() Q-stat p-value 0 0.531426 0.466008 1 3.073217 0.215109 2 3.644229 0.302532 3 3.906326 0.418832 4 4.701433 0.453393 5 5.433745 0.489500 6 5.444254 0.605916 7 5.445309 0.709091 8 5.900762 0.749808 9 6.004928 0.814849 10 6.155966 0.862758 11 6.299958 0.900213 12 12.731542 0.468755 13 14.707894 0.398410 14 20.720607 0.145996 15 23.197433 0.108558 16 23.949801 0.120805 17 24.119236 0.151160 18 25.616184 0.141243 19 26.035165 0.164654 20 28.969880 0.114727 21 28.973660 0.145614 22 29.017716 0.179723 23 32.114006 0.124191 24 32.284805 0.149936 25 33.123395 0.158548 26 33.129059 0.192844 27 33.760488 0.208870 28 38.421053 0.113255 29 38.724226 0.132028 30 38.973426 0.153863 31 38.978172 0.184613 32 39.318954 0.207819 33 39.382472 0.241623 34 39.423763 0.278615 35 40.083689 0.293860 36 43.849515 0.203755 37 45.704476 0.182576 38 47.132911 0.174117 39 47.365305 0.197305 Значение данной статистики и p-values, свидетельствуют о том, что гипотеза о случайности остатков не отвергается, и скорее всего данный процесс предствляет “белый шум”.
Теперь давайте расчитаем коэффициент детерминации $R^2$, чтобы понять какой процент наблюдений описывает данная модель:
pred = model.predict('2013-05-26','2014-12-31', typ='levels') trn = otg['2013-05-26':] r2 = r2_score(trn, pred[1:32]) print 'R^2: %1.2f' % r2 Среднеквадратичное отклонение [2] нашей модели:
Средняя абсолютная ошибка [2] прогноза
Осталось нарисовать наш прогноз на графике:
otg.plot(figsize=(12,6)) pred.plot(style='r--') 
Заключение
Как можно заметить из графика наша модель строит не очень хороший прогноз. Отчасти это связано с выбросами в исходных данных, которые мы не доконца убрали, а также с модулем ARIMA пакета statsmodels , т. к. он достаточно новый. Статья больше направлена на то, чтобы показать как именно можно анализировать временные ряды на python. Также хотелось бы отметить, что в рассмотреном сегодня пакет очень полно реализованы различные методы регрессионного анализ (постараюсь показать в дальнейших статьях).
В целом для небольших исследований пакет statsmodels в полне пригоден, но для серьезной научной работы все же пока сыроват и некоторые тесты и статистики в нем отсутствуют.