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

Введение

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

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

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

pip install rpy2

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

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

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

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

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

%load_ext rmagic

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

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))

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

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

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

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

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

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

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

r_df = com.convert_to_r_dataframe(DataFrame(w_log))

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

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))

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

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).
Общий вид данной модели
В этой модели параметры обозначают следующее:

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

Как определять 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)

Из гарфика PACF видно, что порядок AR будет p=4, а по ACF видно, что порядок MA q = 13, т.к. 13 лаг — это последний лаг отличный от 0.
Теперь перейдем к сезонным составляющим. Для их оценки надо смотреть на лаги кратные размеру сезонности, т.е., если для нашего примера, сезонность 52, то надо рассматривать лаги 52, 104, 156, ...
В нашем случае параметры P и Q будут равны 0 (это видно если посмотреть на ACF и PACF на указанных выше лагах).
В результате наших исследований мы получили модель
Как было указано в начале данной статьи, что найти способы построения данной модели на Python я не нашел, поэтому я принял решение воспользоваться для это функцией arima() из 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-тест Льюнга — Бокса и проверим корреляцию остатков. Для этого в R существует функция tsdiag(), в качестве параметра передается модель и количество лагов для теста.
Вызвать данную функцию можно так:

%Rpush model
%R tsdiag(model, 100)

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

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

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

forecast = importr('forecast')

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

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

%R plot(forecast(model))

Способ 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')

Заключение

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

Автор: kuznetsovin

Источник