- PVSM.RU - https://www.pvsm.ru -

Построение модели SARIMA с помощью Python+R

Введение

Добрый день, уважаемые читатели.
После написания предыдущего поста [1] про анализ временных рядов на Python, я решил исправить замечания, которые были указаны в комментариях, но при их исправлении я столкнулся с рядом проблем, например при построении сезонной модели ARIMA, т.к. подобной функции а пакете statsmodels [2] я не нашел. В итоге я решил использовать для этого функции из R [3], а поиски привели меня к библиотеке rpy2 [4] которая позволяетиспользовать функции из библиотек упомянутого языка.
У многих может возникнуть вопрос «зачем это нужно?», ведь проще просто взять R и выполнить всю работу в нем. Я полность согласен с этим утверждением, но как мне кажется, если данные требуют предварительной обработки, то ее проще произвести на Python, а возможности R использовать при необходимости именно для анализа.
Кроме этого, будет показано как интегрировать результаты выдачи работы функции R в IPython Notebook. [5]

Установка и настройки Rpy2

Для начала работы надо установть rpy2. Сделать это можно с помощью команды:

pip install rpy2

Надо отметить что, для работы данной библиотеки необходим установленный язык R. Скачать его можно с офф. сайта [6].
Следующиее что необходимо сделать, это добавить нужные системные переменные. Для Windows добавить следующие переменные:

  • PATH — путь до R.dll и до R.exe
  • R_HOME — путь до папки в которую установлен R
  • R_USER — имя пользователя под которым загружен windows

Также производилась установка на Mac OS X, там данных манипуляций не потребовалось.

Начало работы

Итак, если вы работаете в IPython Notebook, нужно добавить инструкцию:

%load_ext rmagic

Данное расширение позволяет вызывать некторые функции R через rpy2, и выводит результат прямо в консоль IPython Notebook, что очень удобно (ниже будет показано как это сделать). Подробнее написано здесь [7].
Теперь же загрузим, нужные библиотеки:

from pandas import read_csv, DataFrame, Series
import statsmodels.api as sm
import rpy2.robjects as R
from rpy2.robjects.packages import importr
import pandas.rpy.common as com
from pandas import date_range

Теперь, как и в предыдущей статье, загрузим данные и перейдем к недельным интервалам:

dataset = read_csv('DataSets/tovar_moving.csv',';', index_col=['date'], parse_dates=['date'], dayfirst=True)
otg = dataset.qty
w = otg.resample('w', how='sum')
w.plot(figsize=(12,6))

Построение модели SARIMA с помощью Python+R
Итак, из графика можно заметить ежегодную сезонность (52 недели) и ярко выраженный тренд. Поэтому перед построением модели нам необходимо избавиться от тренда и сезонности.

Предварительный анализ данных

Итак для начала, прологорифмируем исходный ряд, для выравнивания значений:

w_log = log(w)
w_log.plot(figsize=(12,6))

Построение модели SARIMA с помощью Python+R
Как видно у нас в графике присутствует сезонность и соответственно ряд не стационарен. Проверим это с помощью теста Дикки-Фулера [8], который проверяет гипотизу о наличии единичных корней и соотвтвенно если они есть ряд будет не стационарным. Как провести данный тест с помощью библиотеки statsmodels, я показывал в прошлый раз. Сейчас я продемонстрирую как это можно сделать с помощью функции adf.test() [9] из R.
Итак, данная функция находится в R-ой библиотеке tseries [10]. Она предназначена для анализа временных рядов и устанавливается дополнительно. Загрузить нужную библиотеку можно с помощью функции importr() [11].

stats = importr('stats')
tseries = importr('tseries')

Можно заметить, что кроме tseries, мы загрузили еще и библиотеку stats. Она нам понадобитьсядля преобразования типов.
Теперь необходимо перевести данные из типа Python в тип понятный R. Сделать это можно с помощью функции convert_to_r_dataframe() [12] на вход которой подается DataFrame, а на выходе получается вектор для R.

r_df = com.convert_to_r_dataframe(DataFrame(w_log))

Итак, вектор есть следующим шагом надо перевести его в формат временного ряда. Для этого в R существует функция ts() [13], вызов ее будет выглядеть так:

y = stats.ts(r_df)

Предварительная подготовка данных закончена и мы можем вызвать нужную нам функцию:

ad = tseries.adf_test(y, alternative="stationary", k=52)

В качестве параметров ей передается временный ряд и количество лагов, для которых будет расчитываться тест. В экономических моделях принято брать данное значение равное году, а т.к. данные у нас еженедельные, а в году 52 недели, поэтому параметр имеет такое значение.
Теперь в переменной ad содержится R-объект. Его структура описана в виде списка, описание которого мне найти не удалось. Поэтому с помощью визуального анализа я написал код, который выводит результат работы функции в понятном виде:

a = ad.names[:5]
{ad.names[i]:ad[i][0] for i in xrange(len(a))}

{'alternative': 'stationary',
'method': 'Augmented Dickey-Fuller Test',
'p.value': 0.23867869477446427,
'parameter': 52.0,
'statistic': -2.8030060277420006}

Исходя из результатов теста, исходный ряд не стационарен. Т.к. гипотеза о наличии единичных корней принимается с малой вероятностью, и, соответственно, ряд не стационарен. Теперь проверим на стационарность ряд первых разностей.
Для начала получим их с помощью Python, а затем применим ADF-тест:

diff1lev = w.diff(periods=1).dropna()
print 'p.value: %f' % sm.tsa.adfuller(diff1lev, maxlag=52)[1]

p.value: 0.000000

diff1lev.plot(figsize=(12,6))

Построение модели SARIMA с помощью Python+R
Ряд первых разностей, по итогам теста, оказался стационарным. А график помогает нам убедиться, что тенденция отсутствует. Осталось избавиться от сезонности.
Для этого нужно взять сезонную разность, от нашего получившегося ряда. Подробнее можно прочитать тут [14]. Если полученный ряд будет стационармы необходимо будет вязть его превую разность и проверить ее.

diff1lev_season = diff1lev.diff(52).dropna()

Проверим ее на стационарность с помощью ADF-тест из R:

r_df = com.convert_to_r_dataframe(DataFrame(diff1lev_season))
y = stats.ts(r_df)
ad = tseries.adf_test(y, alternative="stationary", k=52)
a = ad.names[:5]
{ad.names[i]:ad[i][0] for i in xrange(len(a))}

{'alternative': 'stationary',
'method': 'Augmented Dickey-Fuller Test',
'p.value': 0.551977997289418,
'parameter': 52.0,
'statistic': -2.0581183466564776}

Итак, ряд не стационарен, возьмем его первые разности:

diff1lev_season1lev = diff1lev_season.diff().dropna()
print 'p.value: %f' % sm.tsa.adfuller(diff1lev_season1lev, maxlag=52)[1]

p.value: 0.000395

Получившийся ряд стационарен. Теперь можно перейти к построению модели

Построение модели.

Итак, предварительный анализ закончен, и мы можем перейти к построению сезонной модели ARIMA (SARIMA).
Общий вид данной модели Построение модели SARIMA с помощью Python+R
В этой модели параметры обозначают следующее:

  • Построение модели SARIMA с помощью Python+R — порядок модели Построение модели SARIMA с помощью Python+R
  • Построение модели SARIMA с помощью Python+R — порядок интегрирования исходных данных
  • Построение модели SARIMA с помощью Python+R — порядок модели Построение модели SARIMA с помощью Python+R
  • Построение модели SARIMA с помощью Python+R — порядок сезонной составляющей Построение модели SARIMA с помощью Python+R
  • Построение модели SARIMA с помощью Python+R — порядок интегрирования сезонной составляющей
  • Построение модели SARIMA с помощью Python+R — порядок сезонной составляющей Построение модели SARIMA с помощью Python+R
  • Построение модели SARIMA с помощью Python+R — размерность сезонности(месяц, квартал и т.д.)

Как определять p, d, q, я показывал в прошлый раз. Сейчас я опишу определять порядок сезонных составляющих P,D,Q.
Начнем с определения параметра D. Он определет порядок интегрированности сезонной разности, т.е. в нашем случае он равен 1. Для определения P и Q нам как и прежде надо построить коррелограммы ACF и PACF.

fig = plt.figure(figsize=(12,8))
ax1 = fig.add_subplot(211)
fig = sm.graphics.tsa.plot_acf(diff1lev_season1lev.values.squeeze(), lags=150, ax=ax1)
ax2 = fig.add_subplot(212)
fig = sm.graphics.tsa.plot_pacf(diff1lev_season1lev, lags=150, ax=ax2)

Построение модели SARIMA с помощью Python+R
Из гарфика PACF видно, что порядок AR будет p=4, а по ACF видно, что порядок MA q = 13, т.к. 13 лаг — это последний лаг отличный от 0.
Теперь перейдем к сезонным составляющим. Для их оценки надо смотреть на лаги кратные размеру сезонности, т.е., если для нашего примера, сезонность 52, то надо рассматривать лаги 52, 104, 156, ...
В нашем случае параметры P и Q будут равны 0 (это видно если посмотреть на ACF и PACF на указанных выше лагах).
В результате наших исследований мы получили модель Построение модели SARIMA с помощью Python+R
Как было указано в начале данной статьи, что найти способы построения данной модели на Python я не нашел, поэтому я принял решение воспользоваться для это функцией arima() [15] из R. В качестве параметров ей передаются порядок модели ARIMA и, при необходимости, порядок сезонной составляющей. Но перед вызовом ее необходимо подготовить некоторые данные.
Для начала переведем наш исходный набор в формат R и переведем в формат временного ряда.

r_df = com.convert_to_r_dataframe(DataFrame(w))
y = stats.ts(r_df)

Порядок модели передается в качестве вектора R, поэтому давайте создадим его:

order = R.IntVector((4,1,13))

Так же в качестве параметра сезонной составляющей передается список, который содержит ее порядок и размер периода:

season = R.ListVector({'order': R.IntVector((0,1,0)), 'period' : 52})

Теперь мы готовы к тому, чтобы построить модель:

model = stats.arima(y, order = order, seasonal=season)

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

Проверка адекватности модели.

Итак для проверки адекватности модели надо проверить соответвуют ли остатки модели «белому шуму». Проверим это проведя Q-тест Льюнга — Бокса [16] и проверим корреляцию остатков. Для этого в R существует функция tsdiag(), в качестве параметра передается модель и количество лагов для теста.
Вызвать данную функцию можно так:

%Rpush model
%R tsdiag(model, 100)

Построение модели SARIMA с помощью Python+R
В первой строке инструкция %Rpush загружает объекты для использования в R. Иструкция %R во второй строке вызывает код в формате языка R. Данная конструкция работает в IPython Notebook.
Из графиков выше можно заметить, что остатки независимы (это видно по ACF). Кроме того из графика Q-статистики можно заметить что во всех точках значение p-value больше уровня значимости, из этого можно сделать вывод, что остатки, с большой вероятностью, являются «белым шумом».

Прогнозирование

Для прогонизования нужно подгрузить библиотеку forecast [17]

forecast = importr('forecast')

Вывести результы прогнозирования можно двумя способами.

Способ 1. Это использовать возможности интеграции IPython и R. Который был показан в предыдущем разделе:

%R plot(forecast(model))

Построение модели SARIMA с помощью Python+R

Способ 2. Второй способ это сделать прогноз с помощью библиотеки forecast, а потом перевести результат в во временную серию pandas и вывести их на экран. Код который будет выполнять написан ниже:

f = forecast.forecast(model) #строим прогноз
pred = [i[1] for i in f[3].iteritems()] #парсим результат
dt = date_range(w.index[-1], periods=len(pred)+1,freq='W')[1:] #создаем индекс из дат
pr = Series(pred, index = dt)  #записываем данные в TimeSeries

Теперь выведем результат на экран

w.plot(figsize=(12,6))
pr.plot(color = 'red')

Построение модели SARIMA с помощью Python+R

Заключение

В качестве заключения хотелось бы отметить, что для более точного анализа сезонности нужно иметь данные с 7-10 сезонами. Кроме того, хотелось бы сказать большое спасибо werwooolf [18] за оказанную помощь в процессе подготовки статьи.
В данной статье я постарался показать как строить сезонную модель ARIMA, а также показал как можно использовать связку языков R и Python для анализа данных. Я думаю специалисты и по одному и подругому языку найдут как эффективно применить описанную связку.

Автор: kuznetsovin

Источник [19]


Сайт-источник PVSM.RU: https://www.pvsm.ru

Путь до страницы источника: https://www.pvsm.ru/python/53684

Ссылки в тексте:

[1] предыдущего поста: http://habrahabr.ru/post/207160/

[2] statsmodels: http://statsmodels.sourceforge.net/stable/index.html

[3] R: http://r-project.org/

[4] rpy2: http://rpy.sourceforge.net/rpy2.html

[5] IPython Notebook.: http://ipython.org/notebook.html

[6] офф. сайта: http://www.r-project.org/

[7] здесь: http://nbviewer.ipython.org/github/ipython/ipython/blob/3607712653c66d63e0d7f13f073bde8c0f209ba8/docs/examples/notebooks/rmagic_extension.ipynb

[8] теста Дикки-Фулера: http://ru.wikipedia.org/wiki/Тест_Дики_—_Фуллера

[9] adf.test(): http://www.inside-r.org/packages/cran/tseries/docs/adf.test

[10] tseries: http://cran.r-project.org/web/packages/tseries/index.html

[11] importr(): http://rpy.sourceforge.net/rpy2/doc-dev/html/robjects_rpackages.html?highlight=importr#rpy2.robjects.packages.importr

[12] convert_to_r_dataframe(): http://pandas.pydata.org/pandas-docs/stable/r_interface.html

[13] ts(): http://stat.ethz.ch/R-manual/R-patched/library/stats/html/ts.html

[14] тут: https://onlinecourses.science.psu.edu/stat510/?q=node/67

[15] arima(): http://stat.ethz.ch/R-manual/R-patched/library/stats/html/arima.html

[16] Q-тест Льюнга — Бокса: http://ru.wikipedia.org/wiki/Q-тест_Льюнга_—_Бокса

[17] forecast: http://cran.r-project.org/web/packages/forecast/index.html

[18] werwooolf: http://habrahabr.ru/users/werwooolf/

[19] Источник: http://habrahabr.ru/post/210530/