Главная
Блог разработчиков phpBB
 
+ 17 предустановленных модов
+ SEO-оптимизация форума
+ авторизация через соц. сети
+ защита от спама

Построение модели SARIMA с поддержкой Python R

Anna | 15.06.2014 | нет комментариев

Вступление

Добродушный день, уважаемые читатели.
Позже написания предыдущего поста про обзор временных рядов на 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).
Всеобщий вид данной модели 
В этой модели параметры обозначают следующее:

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

Как определять pdq, я показывал в прошлый раз. Теперь я опишу определять порядок сезонных составляющих 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)

Порядок мод

Источник: programmingmaster.ru
Оставить комментарий
Форум phpBB, русская поддержка форума phpBB
Рейтинг@Mail.ru 2008 - 2017 © BB3x.ru - русская поддержка форума phpBB