Sarima как подобрать параметры

How can I select the best SARIMA model

The aim of this note is to show, using a real data, how to select the best a SARIMA model for a given time series.

I won’t remind here the theory behind SARIMA models. You can either come to my class to learn more about time series �� or read the following books

Robert H. Shumway and David S. Stoffer, Time Series Analysis and Its Applications With R Examples, Springer, 2016
Avishek Pal and PKS Prakash, Practical Time Series Analysis, Birmingham — Mumbai, 2017.

I will use in this tutorial

the data departures available in fpp2 package and I will the series permanent . It’s about Overseas departures from Australia and I will be using the monthly number permanent departures from January 1976 — November 2016.
R packages such as fpp2 , forecast , urca and ggplot2

I start then by importing the data from fpp2 package and representing it in a graph

We can quickly notice that there’s a trend-cycle behavior and some kind of seasonality. I will firstly use auto.arima , from forecast package, command to check quickly if this kind of data can be fitted using SARIMA models. I the estimated AIC is non-negative and if the log likelihood is negative we can conclude that SARIMA models can be good for this data. I will also use auto.sarima command to get an idea about the values of $d$ , $D$ and $T$ .

I can first notice that I can continue using SARIMA models that fits with this data. I won’t use the previous estimated model by auto.arima because I’m not sure that I’m selecting the best model using the auto.arima command. I will differentiate the series until I obtain stationary. Once this step is done, I will run a selection procedure that we provide the best model adjusting my data that I can use it to make predictions.

I will first represent the process $(1-B)(1-B^12)X_t$ obtained by differentiating my original time series $X_t$ .

I can notice from the previous graph that $(1-B)(1-B^12)X_t$ has a stationary behavior. I will confirm this using ADF test and KPSS test too. Both tests can be performed using respectively the commands ur.df and ur.kpss .

Since we notice in the previous graph that the time series $(1-B)(1-B^12)X_t$ has no trend and no drift. I will then choose the type="none" for the ADF test and the type=mu for the KPSS test.

First, I notice that the regression in the ADF test is significant ; a small p-value, an $R^2$ close to 1 and the coefficients are significant. The value of test-statistic is -11.4782. It’s very smaller that the critical value corresponding to the p-value 1%. Then the p-value of the hypothesis test: \[H_0\,:\, \phi=0\;\mbox< vs >\;$H_1\,:\, \phi<0\] Hence $H_0$ is rejected and then I can conclude that $(1-B)(1-B^12)X_t$ can be considered as a stationary process.

Let’s now perform a KPSS test

The value of the test-statistic is .0221 smaller than the critical value corresponding to the p-value 10%. We can also conclude that $(1-B)(1-B^12)X_t$ can be considered as a stationary process.

I can then conclude that $d=D=1$ and $T=12$ . Let’s notice that $T=12$ is expected since we’re dealing with monthly data.

In the next step I will determine the possible values of $p$ , $P$ , $q$ and $Q$ . I will ACF to determine the couples $(q,Q)$ and I will use PACF to determine the couples $(p,P)$ .

I can deduce from the previous graph that the ACF are equal to zero from the order 14 and these ACF are decreasing to zero. Then $q+12Q\leq 13$ . Let’s recall that $q+12Q$ is the degree of the MA operator.

Let’s now draw the PACF

I will then have the same conclusion as in the ACF and I will choose the couples $p+12P$ such that $p+12P\leq 13$ .

I will now start the selection procedure. I will first construct a matrix with the possibles values of $p$ , $d$ , $q$ , $P$ , $D$ , $Q$ and $T$ . Each row corresponds to each case. In the total we have to estimate 256. It’s a lot but it just takes about 10 to 15 minutes to get all results. I will then eliminate the models with negative AIC and then I will eliminate the models that generate not white noise residuals.

Construction of the matrix of the parameters $p$ , $d$ , $q$ , $P$ , $D$ , $Q$ and $T$ :

Estimating the models now. You can notice that I have added the command try in the loop. Indeed I didn’t want that the loops stops when the estimation algorithm doesn’t converge. It will just provide a void object.

I will now check the residuals, I had constructed a test that I use to eliminate the models that not-white-noise residual. Indeed, I perform a Portemanteau test with lags from 1 to 10. The model is rejected if at least one p-value from the later test is smaller than the significance level 5%. I will be using the command Box.test.2 from caschrono package.

Hence the object aa will tell us if the white-noise test was successful or not.

I will now create a vector aic that contains the AIC of the estimated models

Русские Блоги

AR(p) Модель авторегрессии, то есть использовать себя, чтобы вернуться к себе. Основное предположение состоит в том, что текущее значение последовательности зависит от исторической ценности последовательности. p указывает, сколько исторических значений используется для возврата прогнозируемого значения.

Чтобы определить начальный p, нужно посмотреть на PACF Изобразите и найдите самое большое значимое запаздывание по времени, а остальные запаздывания не значимы после p.

MA(q) Модель скользящего среднего предназначена для моделирования ошибки временного ряда и предположения, что текущая ошибка зависит от ошибки с запаздыванием. допустимый ACF Найдите начальное значение на карте.

I(d) Указывает, что порядок интегрирования равен D. Интегрирование происходит потому, что мы сначала дифференцируем временной ряд d раз, чтобы сделать ряд стабильным. Например, парабола дифференцируется дважды, чтобы получить устойчивое (стационарное) ускорение, правильное ускорение оценивается, а интегрирование выполняется дважды, чтобы восстановить исходный временной ряд.

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

S(s) Используется для моделирования сезонности последовательности, s представляет продолжительность сезона.

С учетом сезонности требуются три дополнительных параметра (P, D, Q).

P представляет собой порядок сезонной авторегрессии, выведенный из PACF. В отличие от малого p, необходимо учитывать временной лаг, кратный длине сезона. Например, если продолжительность сезона 24, то на диаграмме pacf следует проверить интенсивность лагов 24,48,72. Если pacf последовательности с лагом 48 является значимым, то P равно 2.

Q и P берутся аналогично, но выбираются через диаграмму ACF.

D представляет порядок сезонной разницы, которая обычно составляет 0 или 1. Сезонная разница составляет 1.

Хорошо, после долгого разговора вы не поняли, тогда мы используем SARIMA для моделирования ��

Импортировать пакет

Функция автокорреляции и функция частичной автокорреляции

Прочитать данные

Временной ряд кликов по объявлению

Всего 216 = 24 × 9 216 = 24 \times 9 2 1 6 = 2 4 × 9 Данные, 9 дней, 24 часа в сутки

Нарисуйте функции автокорреляции и частичной автокорреляции последовательности

Видно, что данные по количеству кликов по рекламе демонстрируют сильную сезонность, с продолжительностью сезона 1 день и 24 часа, поэтому сначала мы делаем разницу и удаляем сезонность

После удаления сезонности вы чувствуете себя намного лучше? ��

Но на графике автокорреляции и графике частичной автокорреляции все еще слишком много значительных временных лагов.

Давайте сделаем еще одно отличие

Выбор параметра SARIMA

p равно 4, потому что 4-ступенчатая задержка PACF выдающаяся, а следующие не значимы.
d равно 1, потому что первое различие сделано
q Согласно ACF, это должно быть около 4
s 24, без сомнения
P может быть 2, потому что гистерезис 24 (1 с) и 48 (2 с) несколько важен для PACF.
D равно 1, потому что мы сделали сезонную разницу
Q может быть 1, 24-е (1 с) отставание ACF очевидно, а 48-е отставание — нет.

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

Но достоверно одно: d = 1, D = 1, s = 24

Все комбинации параметров-кандидатов перечислены ниже, имеется 36 групп параметров-кандидатов.

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

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

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

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

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

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

Читать:

Mozilla background task что это

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

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

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

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

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

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

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

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

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

<'alternative': 'stationary',
‘method’: ‘Augmented Dickey-Fuller Test’,
‘p.value’: 0.23867869477446427,
‘parameter’: 52.0,
‘statistic’: -2.8030060277420006>

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

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

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

<'alternative': 'stationary',
‘method’: ‘Augmented Dickey-Fuller Test’,
‘p.value’: 0.551977997289418,
‘parameter’: 52.0,
‘statistic’: -2.0581183466564776>

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

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

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

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

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

Из гарфика 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, поэтому давайте создадим его:

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

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

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

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

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

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

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

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

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

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

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