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

Обзор временных рядов с поддержкой python

Anna | 15.06.2014 | нет комментариев
Добрый день, уважаемые читатели.
В сегодняшней статье, я попытаюсь описать процесс обзора временных рядов с поддержкой python и модуляstatsmodels. Данный модуль предоставляет широкий комплект средств и способов для проведения статистического обзора и эконометрики. Я попытаюсь показать основные этапы обзора таких рядов, в завершении мы возведем модель ARIMA.
Для примера взяты настоящие данные по товарообороту одного из складских комплексов Подмосковья.

Загрузка и заблаговременная обработка данных

Для начала загрузим данные и посмотрим на них:

from pandas import read_csv, DataFrame
import statsmodels.api as sm
from statsmodels.iolib.table import SimpleTable
from sklearn.metrics import r2_score
import ml_metrics as metrics
In [2]:
dataset = read_csv('tovar_moving.csv',';', index_col=['date_oper'], parse_dates=['date_oper'], dayfirst=True)
dataset.head()
Otgruzka priemka
date_oper
2009-09-01 179667 276712
2009-09-02 177670 164999
2009-09-03 152112 189181
2009-09-04 142938 254581
2009-09-05 130741 192486

Выходит, как дозволено подметить функция read_csv(), в данном случем помимо указания параметров, которые задают используемые колонки и индекс, дозволено подметить еще 3 параметра для работы с датой. Остановимся на них поподробнее.
parse_dates задает имена столбцов, которые будут преобразованы в тип DateTime. Стоит подметить, что если в данном столбце будут пустые значения парсинг не удастся и вернется столбец типа object. Дабы этого избежать нужно добавить параметр keep_default_na=False.
Заключительный параметр dayfirst указывает функции парсинга, что первое в строке первым идет день, а не напротив. Если не задать данный параметр, то функция может не верно преобразовывать даты и путать месяц и день местами. Скажем 01.02.2013 будет преобразовано в 02-01-2013, что будет ненормально.
Выделим в отдельную серию временной ряд со значениями отгрузкок:

otg = dataset.Otgruzka
otg.head()
date_oper
2009-09-01 179667
2009-09-02 177670
2009-09-03 152112
2009-09-04 142938
2009-09-05 130741
Name: Otgruzka, dtype: int64

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

Обзор временного ряда

Для начала давайте посмортим график нашего ряда:

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


Из графика видно, что наш ряд имеет малое кол-во выбросов, которые влияют на разброс. Помимо того исследовать отгрузки за всякий день не вовсе правильно, т.к., скажем, в конце либо начале недели будут дни в которые товара отгружается гораздо огромнее, нежели в остальные. Следственно есть толк перейти к недельному промежутку и среднему значению отгрузок на нем, это избавит нас от выбросов и уменьшит колебания нашего ряда. В pandas для этого есть комфортная функция resample(), в качестве параметров ей передается период округления и аггрегатная функция:

otg = otg.resample('W', how='mean')
otg.plot(figsize=(12,6))


Как дозволено подметить, новейший график не имеет блестящих выбросов и имеет ясно выраженный тренд. Из это дозволено сделать итог о том, что ряд не является стационарным[1].

itog = otg.describe()
otg.hist()
itog
count 225
mean 270858.285365
std 118371.082975
min 872.857143
25% 180263.428571
50% 277898.714286
75% 355587.285714
max 552485.142857
dtype: float64

Как дозволено подметить из колляций и гистограммы, ряд у нас больше менее однородный и имеет касательно маленький разброс о чем свидетельствует показатель вариации, где  — cреднеквадратическое отклонение — среднее арифметическое выборки. В нашем случае он равен:

print 'V = %f' % (itog['std']/itog['mean'])

V = 0.437022

Проведем тест Харки — Бера для определения номарльности разделения, Дабы удостоверить предположение об однородности. Для этого в существует функция jarque_bera(), которая возвращает значения данной статистики:

row =  [u'JB', u'p-value', u'skew', u'kurtosis']
jb_test = sm.stats.stattools.jarque_bera(otg)
a = np.vstack([jb_test])
itog = SimpleTable(a, row)
print itog

Значение данной статистика свидетельствует о том, нулевая догадка о нормальности разделения отвергается с малой вероятностью (probably > 0.05), и, следственно, наш ряд имеет типичного разделения.
Функция SimpleTable() служит для оформления итога. В нашем случае на вход ей подается массив значений (размерность не огромнее 2) и список с наименованиями столбцов либо строк.
Многие способы и модели основаны на предположениях о стационарн%D0%B4%D0%B5%D0%BB%D1%8C”>AR

  • d — порядок интегрированного ряда
  • q — порядок компонетны MA

Параметр d есть и он равет 1, осталось определить p и q. Для их определения нам нужно исследоватьавторкорреляционную(ACF) и Отчасти автокорреляционную(PACF) функции для ряда первых разностей.
ACF поможет нам определить q, т. к. по ее коррелограмме дозволено определить число автокорреляционных показателей крепко хороших от 0 в модели MA
PACF поможет нам определить p, т. к. по ее коррелограмме дозволено определить наивысший номер показателя крепко чудесный от 0 в модели AR.
Дабы возвести соответствующие коррелограммы, в пакете statsmodels имеются следующие функции:plot_acf() и plot_pacf(). Они выводят графики ACF и PACF, у которых по оси X откладываются номера лагов, а по оси Y значения соответствующих функций. Необходимо подметить, что число лагов в функциях и определяет число важных показателей. Выходит, наши функции выглядят так:

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


Позже постижения коррелограммы PACF дозволено сделать итог, что p = 1, т.к. на ней только 1 лаг крепко отличнен от нуля. По коррелограмме ACF дозволено увидеть, что q = 1, т.к. позже лага 1 значении функций круто падают.
Выходит, когда вестимы все параметры дозволено возвести модель, но для ее построения мы возмем не все данные, а только часть. Данные из части не попавших в модель мы оставим для проверки точности прогноза нашей модели:

src_data_model = otg[:'2013-05-26']
model = sm.tsa.ARIMA(src_data_model, order=(1,1,1), freq='W').fit(full_output=False, disp=0)

Параметр trend отвечает за присутствие константы в моделе. Выведем информацаю по получившейся модели:

print model.summary()

Как видно из данной информации в нашей модели все показатели важные и дозволено перейти к оценке модели.

Обзор и оценка модели

Проверим остатки данной модели на соответствие «белому шуму», а также проанализируем коррелограму остатков, так как это может нам подмогнуть в определении значимых для включения и прогнозирования элементов регрессии.
Выходит первое, что мы сделаем это проведем Q-тест Льюнга — Бокса для проверки догадки о том, что остатки случайны, т. е. являются «белым шумом». Данный тест проводится на остатках модели ARIMA. Таким образом, нам нужно вначале получить остатки модели и возвести для них ACF, а после этого к получившимся показателям приметить тест. С поддержкой statsmadels это дозволено сделать так:

q_test = sm.tsa.stattools.acf(model.resid, qstat=True) #_rqvmk!
 35
40.083689
0.293860
 36
43.849515
0.203755
 37
45.704476
0.182576
 38
47.132911
0.174117
 39
47.365305
0.197305


Значение данной статистики и p-values, свидетельствуют о том, что догадка о случайности остатков не отвергается, и скорее каждого данный процесс предствляет «белый шум».
Сейчас давайте расчитаем показатель детерминации , Дабы осознать какой процент слежений описывает данная модель:

pred = model.predict('2013-05-26','2014-12-31', typ='levels')
trn = otg['2013-05-26':]
r2 = r2_score(trn, pred[1:32])
print 'R^2: %1.2f' % r2


R^2: -0.03

Среднеквадратичное отклонение[2] нашей модели:

metrics.rmse(trn,pred[1:32])


80919.057367642512

Средняя безусловная оплошность[2] прогноза:

metrics.mae(trn,pred[1:32])


63092.763277651895

Осталось нарисовать наш прогноз на графике:

otg.plot(figsize=(12,6))
pred.plot(style='r--')


Завершение


Как дозволено подметить из графика наша модель строит не дюже отличный прогноз. Отчасти это связано с выбросами в начальных данных, которые мы не доконца убрали, а также с модулем ARIMA пакета statsmodels, т. к. он довольно новейший. Статья огромнее направлена на то, Дабы показать как именно дозволено исследовать временные ряды на python. Также хотелось бы подметить, что в рассмотреном сегодня пакет дюже полно реализованы разные способы регрессионного обзор (постараюсь показать в дальнейших статьях).
В целом для маленьких изысканий пакет statsmodels в полне пригоден, но для серьезной научной работы все же пока сыроват и некоторые тесты и статистики в нем отсутствуют.

Ссылки


  1. И.И. Елисеева. Эконометрика
  2. Сопоставление моделей временных рядов
Добрый день, уважаемые читатели.
В сегодняшней статье, я попытаюсь описать процесс обзора временных рядов с поддержкой python и модуляstatsmodels. Данный модуль предоставляет широкий комплект средств и способов для проведения статистического обзора и эконометрики. Я попытаюсь показать основные этапы обзора таких рядов, в завершении мы возведем модель ARIMA.
Для примера взяты настоящие данные по товарообороту одного из складских комплексов Подмосковья.

Загрузка и заблаговременная обработка данных

Для начала загрузим данные и посмотрим на них:

from pandas import read_csv, DataFrame
import statsmodels.api as sm
from statsmodels.iolib.table import SimpleTable
from sklearn.metrics import r2_score
import ml_metrics as metrics
In [2]:
dataset = read_csv('tovar_moving.csv',';', index_col=['date_oper'], parse_dates=['date_oper'], dayfirst=True)
dataset.head()
Otgruzka priemka
date_oper
2009-09-01 179667 276712
2009-09-02 177670 164999
2009-09-03 152112 189181
2009-09-04 142938 254581
2009-09-05 130741 192486

Выходит, как дозволено подметить функция read_csv(), в данном случем помимо указания параметров, которые задают используемые колонки и индекс, дозволено подметить еще 3 параметра для работы с датой. Остановимся на них поподробнее.
parse_dates задает имена столбцов, которые будут преобразованы в тип DateTime. Стоит подметить, что если в данном столбце будут пустые значения парсинг не удастся и вернется столбец типа object. Дабы этого избежать нужно добавить параметр keep_default_na=False.
Заключительный параметр dayfirst указывает функции парсинга, что первое в строке первым идет день, а не напротив. Если не задать данный параметр, то функция может не верно преобразовывать даты и путать месяц и день местами. Скажем 01.02.2013 будет преобразовано в 02-01-2013, что будет ненормально.
Выделим в отдельную серию временной ряд со значениями отгрузкок:

otg = dataset.Otgruzka
otg.head()
date_oper
2009-09-01 179667
2009-09-02 177670
2009-09-03 152112
2009-09-04 142938
2009-09-05 130741
Name: Otgruzka, dtype: int64

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

Обзор временного ряда

Для начала давайте посмортим график нашего ряда:

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


Из графика видно, что наш ряд имеет малое кол-во выбросов, которые влияют на разброс. Помимо того исследовать отгрузки за всякий день не вовсе правильно, т.к., скажем, в конце либо начале недели будут дни в которые товара отгружается гораздо огромнее, нежели в остальные. Следственно есть толк перейти к недельному промежутку и среднему значению отгрузок на нем, это избавит нас от выбросов и уменьшит колебания нашего ряда. В pandas для этого есть комфортная функция resample(), в качестве параметров ей передается период округления и аггрегатная функция:

otg = otg.resample('W', how='mean')
otg.plot(figsize=(12,6))


Как дозволено подметить, новейший график не имеет блестящих выбросов и имеет ясно выраженный тренд. Из это дозволено сделать итог о том, что ряд не является стационарным[1].

itog = otg.describe()
otg.hist()
itog
count 225
mean 270858.285365
std 118371.082975
min 872.857143
25% 180263.428571
50% 277898.714286
75% 355587.285714
max 552485.142857
dtype: float64

Как дозволено подметить из колляций и гистограммы, ряд у нас больше менее однородный и имеет касательно маленький разброс о чем свидетельствует показатель вариации, где  — cреднеквадратическое отклонение — среднее арифметическое выборки. В нашем случае он равен:

print 'V = %f' % (itog['std']/itog['mean'])

V = 0.437022

Проведем тест Харки — Бера для определения номарльности разделения, Дабы удостоверить предположение об однородности. Для этого в существует функция jarque_bera(), которая возвращает значения данной статистики:

row =  [u'JB', u'p-value', u'skew', u'kurtosis']
jb_test = sm.stats.stattools.jarque_bera(otg)
a = np.vstack([jb_test])
itog = SimpleTable(a, row)
print itog

Значение данной статистика свидетельствует о том, нулевая догадка о нормальности разделения отвергается с малой вероятностью (probably > 0.05), и, следственно, наш ряд имеет типичного разделения.
Функция SimpleTable() служит для оформления итога. В нашем случае на вход ей подается массив значений (размерность не огромнее 2) и список с наименованиями столбцов либо строк.
Многие способы и модели основаны на предположениях о стационарн%D0%B4%D0%B5%D0%BB%D1%8C”>AR

  • d — порядок интегрированного ряда
  • q — порядок компонетны MA

Параметр d есть и он равет 1, осталось определить p и q. Для их определения нам нужно исследоватьавторкорреляционную(ACF) и Отчасти автокорреляционную(PACF) функции для ряда первых разностей.
ACF поможет нам определить q, т. к. по ее коррелограмме дозволено определить число автокорреляционных показателей крепко хороших от 0 в модели MA
PACF поможет нам определить p, т. к. по ее коррелограмме дозволено определить наивысший номер показателя крепко чудесный от 0 в модели AR.
Дабы возвести соответствующие коррелограммы, в пакете statsmodels имеются следующие функции:plot_acf() и plot_pacf(). Они выводят графики ACF и PACF, у которых по оси X откладываются номера лагов, а по оси Y значения соответствующих функций. Необходимо подметить, что число лагов в функциях и определяет число важных показателей. Выходит, наши функции выглядят так:

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


Позже постижения коррелограммы PACF дозволено сделать итог, что p = 1, т.к. на ней только 1 лаг крепко отличнен от нуля. По коррелограмме ACF дозволено увидеть, что q = 1, т.к. позже лага 1 значении функций круто падают.
Выходит, когда вестимы все параметры дозволено возвести модель, но для ее построения мы возмем не все данные, а только часть. Данные из части не попавших в модель мы оставим для проверки точности прогноза нашей модели:

src_data_model = otg[:'2013-05-26']
model = sm.tsa.ARIMA(src_data_model, order=(1,1,1), freq='W').fit(full_output=False, disp=0)

Параметр trend отвечает за присутствие константы в моделе. Выведем информацаю по получившейся модели:

print model.summary()

Как видно из данной информации в нашей модели все показатели важные и дозволено перейти к оценке модели.

Обзор и оценка модели

Проверим остатки данной модели на соответствие «белому шуму», а также проанализируем коррелограму остатков, так как это может нам подмогнуть в определении значимых для включения и прогнозирования элементов регрессии.
Выходит первое, что мы сделаем это проведем Q-тест Льюнга — Бокса для проверки догадки о том, что остатки случайны, т. е. являются «белым шумом». Данный тест проводится на остатках модели ARIMA. Таким образом, нам нужно вначале получить остатки модели и возвести для них ACF, а после этого к получившимся показателям приметить тест. С поддержкой statsmadels это дозволено сделать так:

q_test = sm.tsa.stattools.acf(model.resid, qstat=True) #_rqvmk!
 35
40.083689
0.293860
 36
43.849515
0.203755
 37
45.704476
0.182576
 38
47.132911
0.174117
 39
47.365305
0.197305


Значение данной статистики и p-values, свидетельствуют о том, что догадка о случайности остатков не отвергается, и скорее каждого данный процесс предствляет «белый шум».
Сейчас давайте расчитаем показатель детерминации , Дабы осознать какой процент слежений описывает данная модель:

pred = model.predict('2013-05-26','2014-12-31', typ='levels')
trn = otg['2013-05-26':]
r2 = r2_score(trn, pred[1:32])
print 'R^2: %1.2f' % r2


R^2: -0.03

Среднеквадратичное отклонение[2] нашей модели:

metrics.rmse(trn,pred[1:32])


80919.057367642512

Средняя безусловная оплошность[2] прогноза:

metrics.mae(trn,pred[1:32])


63092.763277651895

Осталось нарисовать наш прогноз на графике:

otg.plot(figsize=(12,6))
pred.plot(style='r--')


Завершение


Как дозволено подметить из графика наша модель строит не дюже отличный прогноз. Отчасти это связано с выбросами в начальных данных, которые мы не доконца убрали, а также с модулем ARIMA пакета statsmodels, т. к. он довольно новейший. Статья огромнее направлена на то, Дабы показать как именно дозволено исследовать временные ряды на python. Также хотелось бы подметить, что в рассмотреном сегодня пакет дюже полно реализованы разные способы регрессионного обзор (постараюсь показать в дальнейших статьях).
В целом для маленьких изысканий пакет statsmodels в полне пригоден, но для серьезной научной работы все же пока сыроват и некоторые тесты и статистики в нем отсутствуют.

Ссылки


  1. И.И. Елисеева. Эконометрика
  2. Сопоставление моделей временных рядов

 

 

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