Добрый день, уважаемые читатели.

В сегодняшней статье, я попытаюсь описать процесс анализа временных рядов с помощью 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()
Otgruzka priemka
date_oper
2009-09-01 179667 276712
2009-09-02 177670 164999
2009-09-03 152112 189181
2009-09-04 142938 254581
2009-09-05 130741 192486

Итак, как можно заметить функция 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

Итак у нас теперь есть временной ряд и можно перейти к его анализу.

Анализ временного ряда

Для начала давайте посмортим график нашего ряда:

otg.plot(figsize=(12,6))
<matplotlib.axes.AxesSubplot at 0x10843ce50>

png

Из графика видно, что наш ряд имеет небольшое кол-во выбросов, которые влияют на разброс. Кроме того анализировать отгрузки за каждый день не совсем верно, т.к., например, в конце или начале недели будут дни в которые товара отгружается значительно больше, нежели в остальные. Поэтому есть смысл перейти к месячному интервалу и среднему значению отгрузок на нем, это избавит нас от выбросов и уменьшит колебания нашего ряда. В pandas для этого есть удобная функция resample(), в качестве параметров ей передается период округления и аггрегатная функция:

otg = otg.resample('W', how='mean')
otg.plot(figsize=(12,6))
<matplotlib.axes.AxesSubplot at 0x1085ddad0>

png

Как можно заметить, новый график не имеет ярких выбросов и имеет ярко выраженный тренд. Из это можно сделать вывод о том, что ряд не является стационарным[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

png

Как можно заметить из характеристик и гистограммы, ряд у нас более менее однородный и имеет относительно небольшой разброс о чем свидетельствует коэффициент вариации: $V = \frac {\sigma}{\bar{x}} $, где $\sigma$ - cреднеквадратическое отклонение, $\bar{x}$ - среднее арифметическое выборки. В нашем случае он равен:

print 'V = %f' % (itog['std']/itog['mean'])
V = 0.437022

Проведем тест Харки — Бера для определения номарльности распределения, чтобы подтвердить предположение об однородности. Для этого в существует функция 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:  {'5%': -2.8753374677799957, '1%': -3.4617274344627398, '10%': -2.5741240890815571}
есть единичные корни, ряд не стационарен

Проведенный тест подтвердил предположения о не стационарности ряда. Во многих случаях взятие разности рядов позволяет это сделать.

Если, например, первые разности ряда стационарны, то он называется интегрированным рядом первого порядка.

Итак, давайте определим порядок интегрированного ряда для нашего ряда:

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:  {'5%': -2.8755379867788462, '1%': -3.4621857592784546, '10%': -2.574231080806213}
единичных корней нет, ряд стационарен

Как видно из кода выше получившийся ряд первых разностей приблизился к стационарному. Для полной уверенности разобъем его на несколько промежутков и убедимся мат. ожидания на разных интервалах:

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:  0.693072039563

Высокое p-value дает нам возможность утверждать, что нулевая гипотеза о равенстве средних верна, что свидетельствует о стационарности ряда. Осталось убедиться в отсутствии тренда для этого построим график нашего нового ряда:

otg1diff.plot(figsize=(12,6))
<matplotlib.axes.AxesSubplot at 0x1084d64d0>

png

Тренд действительно отсутствует, таким образом ряд первых разностей является стационарным, а наш исходный ряд - интегрированным рядом первого порядка.

Построение модели временного ряда

Для моделирования будем использовать модель ARIMA, построенную для ряда первых разностей.

Итак, чтобы построить модель нам нужно знать ее порядок, состоящий из 2-х параметров:

  1. p - порядок компоненты AR
  2. d - порядок интегрированного ряда
  3. q - порядок компонетны MA

Параметр 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)

png

После изучения коррелограммы 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 отвечает за наличие константы в моделе. Выведем информацаю по получившейся модели:

print model.summary()
                             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':q_test[1], 'p-value':q_test[2]})
       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
R^2: -0.03

Среднеквадратичное отклонение[2] нашей модели:

metrics.rmse(trn,pred[1:32])
80919.057367642512

Средняя абсолютная ошибка[2] прогноза

metrics.mae(trn,pred[1:32])
63092.763277651895

Осталось нарисовать наш прогноз на графике:

otg.plot(figsize=(12,6))
pred.plot(style='r--')
<matplotlib.axes.AxesSubplot at 0x109939290>

png

Заключение

Как можно заметить из графика наша модель строит не очень хороший прогноз. Отчасти это связано с выбросами в исходных данных, которые мы не доконца убрали, а также с модулем ARIMA пакета statsmodels, т. к. он достаточно новый. Статья больше направлена на то, чтобы показать как именно можно анализировать временные ряды на python. Также хотелось бы отметить, что в рассмотреном сегодня пакет очень полно реализованы различные методы регрессионного анализ (постараюсь показать в дальнейших статьях).

В целом для небольших исследований пакет statsmodels в полне пригоден, но для серьезной научной работы все же пока сыроват и некоторые тесты и статистики в нем отсутствуют.

Ссылки

  1. И.И. Елисеева. Эконометрика
  2. Сравнение моделей временных рядов