- PVSM.RU - https://www.pvsm.ru -
Добрый день, уважаемые читатели.
После написания предыдущего поста [1] про анализ временных рядов на Python, я решил исправить замечания, которые были указаны в комментариях, но при их исправлении я столкнулся с рядом проблем, например при построении сезонной модели ARIMA, т.к. подобной функции а пакете statsmodels [2] я не нашел. В итоге я решил использовать для этого функции из R [3], а поиски привели меня к библиотеке rpy2 [4] которая позволяетиспользовать функции из библиотек упомянутого языка.
У многих может возникнуть вопрос «зачем это нужно?», ведь проще просто взять R и выполнить всю работу в нем. Я полность согласен с этим утверждением, но как мне кажется, если данные требуют предварительной обработки, то ее проще произвести на Python, а возможности R использовать при необходимости именно для анализа.
Кроме этого, будет показано как интегрировать результаты выдачи работы функции R в IPython Notebook. [5]
Для начала работы надо установть rpy2. Сделать это можно с помощью команды:
pip install rpy2
Надо отметить что, для работы данной библиотеки необходим установленный язык R. Скачать его можно с офф. сайта [6].
Следующиее что необходимо сделать, это добавить нужные системные переменные. Для 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))
Итак, из графика можно заметить ежегодную сезонность (52 недели) и ярко выраженный тренд. Поэтому перед построением модели нам необходимо избавиться от тренда и сезонности.
Итак для начала, прологорифмируем исходный ряд, для выравнивания значений:
w_log = log(w)
w_log.plot(figsize=(12,6))
Как видно у нас в графике присутствует сезонность и соответственно ряд не стационарен. Проверим это с помощью теста Дикки-Фулера [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))
Ряд первых разностей, по итогам теста, оказался стационарным. А график помогает нам убедиться, что тенденция отсутствует. Осталось избавиться от сезонности.
Для этого нужно взять сезонную разность, от нашего получившегося ряда. Подробнее можно прочитать тут [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).
Общий вид данной модели
В этой модели параметры обозначают следующее:
Как определять 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() [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)
В первой строке инструкция %Rpush загружает объекты для использования в R. Иструкция %R во второй строке вызывает код в формате языка R. Данная конструкция работает в IPython Notebook.
Из графиков выше можно заметить, что остатки независимы (это видно по ACF). Кроме того из графика Q-статистики можно заметить что во всех точках значение p-value больше уровня значимости, из этого можно сделать вывод, что остатки, с большой вероятностью, являются «белым шумом».
Для прогонизования нужно подгрузить библиотеку forecast [17]
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 [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/
Нажмите здесь для печати.