Расчет абсолютной ошибки временного ряда осуществляется для оценки полученной модели - Исправление ошибок и поиск оптимальных решений проблем

Временные ряды в эконометрических исследованиях

Большинство
эконометрических моделей строится как
динамические эконометрические модели.
Это означает, что моделирование
причинно-следственных связей между
переменными осуществляется во времени,
а исходные данные представлены в форме
временных рядов.

Временной ряд х_t
(t=1;n)
– ряд значений какого-либо показателя
за несколько последовательных промежутков
времени.

Каждый временной
ряд х_t
складывается
из следующих основных составляющих
(компонентов):

1) Тенденции,
характеризующей общее направление
динамики изучаемого явления. Аналитически
тенденция выражается некоторой функцией
времени, называемой трендом (Т).

2) Циклической или
периодической составляющей, характеризующей
циклические или периодические колебания
изучаемого явления. Колебания представляют
собой отклонения фактических уровней
ряда от тренда. Объем продаж некоторых
товаров подвержен сезонным колебаниям.
Сезонные
колебания
(S)
–
периодические колебания, которые имеют
определенный и постоянный период равный
годовому промежутку. Конъюнктурные
колебания (К)
связаны с большими экономическими
циклами, период таких колебаний –
несколько лет.

3) Случайной
составляющей, которая является результатом
воздействия множества случайных факторов
(Е).

Тогда уровень ряда
можно представить как функцию от этих
составляющих (компонентов):
=f(T,
K, S, E).

В зависимости от
взаимосвязи между составляющими может
быть построена либо аддитивная модель:
=T+K+S+E,
либо мультипликативная модель:=T·K·S·E
ряда динамики.

Для определения
состава компонентов (структуры временного
ряда) в
модели временного ряда строят
автокорреляционную функцию.

Автокорреляция
– корреляционная
связь между последовательными уровнями
одного и того же ряда динамики (сдвинутыми
на определенный промежуток времени L
— лаг). То есть, автокорреляция — это связь
между рядом: x₁,
x₂,
… x_n-_l
и рядом x_1+l,
x_2+l,
…,x_n,
где L—
положительное целое число. Автокорреляция
может быть измерена коэффициентом
автокорреляции:

где
,

–средний уровень
ряда (x₁₊_L,
x₂₊_L,…,
x_n
),

средний уровень
ряда (x₁,
x₂,…,
x_n-L
),

_t
, _t-L
– средние квадратические отклонения,
для рядов (x_1+L,
x_2+L,…,
x_n
) и (x₁,
x₂,…,
x_n-L
) соответственно.

Лаг (сдвиг во
времени) определяет порядок коэффициента
автокорреляции. Если L=1,
то имеем коэффициент автокорреляции
1-ого порядка r_t,t-1,
если L=2,
то коэффициент автокорреляции 2-ого
порядка r_t,t-₂и т.д.
Следует учитывать, что с увеличением
лага на единицу, число пар значений, по
которым рассчитывается коэффициент
автокорреляции уменьшается на 1. Поэтому
обычно рекомендуют максимальный порядок
коэффициента автокорреляции равный
n/4.

Рассчитав несколько
коэффициентов автокорреляции, можно
определить лаг (L),
при котором автокорреляция (r_t,t-L)
наиболее высокая, выявив тем самым
структуру
временного ряда.
Если наиболее высоким оказывается
значение r_t,t-₁,
то исследуемый ряд содержит только
тенденцию. Если наиболее высоким оказался
r_t,t-L,
то ряд содержит колебания периодом L.
Если ни один из r_t,t-L
не является
значимым, можно сделать одно из двух
предположений:

— либо ряд не
содержит тенденции и циклических
колебаний, а его уровень определяется
только случайной компонентой;

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

Последовательность
коэффициентов автокорреляции 1, 2 и т.д.
порядков называют автокорреляционной
функцией временного ряда. График
зависимости значений коэффициентов
автокорреляции от величины лага (порядка
коэффициента автокорреляции) называют
коррелограммой.

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

Рассмотрим на
примере как построить коррелограмму,
чтобы определяется структуру временного
ряда.

Пусть нам даны
поквартальные данные
об объеме выпуска некоторого товара
некоторой фирмой –х
(усл.ед.) за 3 года:

1993	1994	1995
1	2	3	4	1	2	3	4	1	2	3	4
410	560	715	500	520	740	975	670	705	950	1200	900

Чтобы построить
коррелогорамму для нашего примера,
исходный ряд динамики дополним рядами
из уровней этого ряда, сдвинутыми во
времени (таблица 6).

Таблица 6

t	1	2	3	4	5	6	7	8	9	10	11	12
х_t	—	560	715	500	520	740	975	670	705	950	1200	900	r_t,t-1=0,537
x_t-1	—	410	560	715	500	520	740	975	670	705	950	1200
х_t	—	—	715	500	520	740	975	670	705	950	1200	900	r_t,t-2=0,085
х_t-2	—	—	410	560	715	500	520	740	975	670	705	950
х_t	—	—	—	500	520	740	975	670	705	950	1200	900	r_t,t-3=0,445
х_t-3	—	—	—	410	560	715	500	520	740	975	670	705
х_t	—	—	—	—	520	740	975	670	705	950	1200	900	r_t,t-4=0,990
х_t-4	—	—	—	—	410	560	715	500	520	740	975	670
х_t	—	—	—	—	—	740	975	670	705	950	1200	900	r_t,t-5=0,294
х_t-5	—	—	—	—	—	410	560	715	500	520	740	975

Рассчитаем
коэффициенты корреляции:

1-ого порядка для
рядов х_t
и х_t_-1,

2-ого порядка для
рядов х_t
и х_t_-2,

3-его порядка для
рядов х_t
и х_t_-3,

4-ого порядка для
рядов х_t
и х_t_-4,

5-ого порядка для
рядов х_t
и х_t_-5

Результаты
расчетов представлены в таблице 7.

Таблица 7

Лаг (порядок) – L	r_t,t-L	Коррелограмма
1	0,537	****
2	0,085	*
3	0,445	***
4	0,990	*****
5	0,294	**

Вывод: в данном
ряду динамики имеется тенденция (т.к.
r_t,t-1=0,537
→1) и периодические колебания с периодом
(L)
равным 4, т.е. имеют место сезонные
колебания (т.к. r_t,t-4=0,99
→1).

Построение
модели временного ряда с сезонными
колебаниями (аддитивная модель).

Процесс построения
модели временного ряда (х),
содержащего n
уровней
некоторого показателя за Z
лет, с L
сезонными колебаниями включает следующие
шаги:

1) Выравнивание
исходного ряда методом скользящей
средней
(х^c).
Произведем выравнивание исходного ряда
взятого из примера, рассмотренного
выше, методом скользящей средней с
периодом усреднения равным 3. Результаты
представлены в таблице 9 (столбец 4).

2) Расчет
значений сезонной составляющей
S_i,
i=1;L,
где L–
число сезонов в году. Для нашего примера
L=4
(сезоны — кварталы).

Расчет значений
сезонных составляющих осуществляется
после устранения тенденции из исходных
уровней ряда: x—
x^c
(столбец 5,
таблица 9). Для дальнейшего расчета S_i
построим отдельную таблицу. Строки
данной таблицы соответствуют сезонам,
столбцы — годам. В теле таблицы находятся
значения: x—
x^c.
По этим данным рассчитываются средние
оценки сезонных составляющих каждой
строке (S^c_i).
Если сумма всех средних оценок равна
нулю (),
то данные средние и будут окончательными
значениями сезонных составляющих
(S_i=S^c_i).
Если их сумма не равна нулю, то
рассчитываются скорректированные
значения сезонных составляющих вычитанием
из средней оценки величины равной
отношению суммы средних оценок к их
общему числу ().
Для нашего примера расчет значений S_i
представлен в таблице 8.

Таблица 8

Номер сезона	Год 1	Год 2	Год 3	Средняя оценка сезонной составляющей	Скорректированная оценка сезонной составляющей S_i
1	—	-66,67	-70,00	-68,33	-67,15
2	-1,67	-5,00	-1,67	-2,78	-1,60
3	123,33	180,00	183,33	162,22	163,40
4	-78,33	-113,33	—	-95,83	-94,66
Итого				-4,72	0

3) Устранение
влияния сезонной составляющей из
исходного ряда динамики :
x^S= x—S_i.
Результаты расчета x^S
для нашего примера представлены в
столбце 6 таблицы 9.

4) Аналитическое
выравнивание уровней
x^S
(построение тренда):
.

Расчет параметров
при аналитическом выравнивании чаще
всего производится с помощью метода
наименьших квадратов (МНК). При этом
поиск параметров для линейного уравнения
тренда можно упростить, если отсчет
времени производить так, чтобы сумма
показателей времени изучаемого ряда
динамики была равна нулю. Для этого
вводится новая условная переменная
времени t^y,
такая, что t^y
=0. Уравнение тренда при этом будет
следующим:
.

При нечетном числе
уровней ряда динамики для получения 
t^y=0
уровень, находящийся в середине ряда,
принимается за условное начало отсчета
времени (периоду или моменту времени,
соответствующему данному уровню
присваивается нулевое значение). Даты
времени, расположенные левее этого
уровня, обозначаются натуральными
числами со знаком минус (-1 –2 –3 …), а
даты времени, расположенные правее
этого уровня – натуральными числами
со знаком плюс (1 2 3 …).

Если число уровней
ряда четное, периоды времени левой
половины ряда (до середины) нумеруются
–1, -3, -5 и т.д. А периоды правой половины
— +1, +3, +5 и.т.д. При этом t^y
будет равна 0.

Система нормальных
уравнений (соответствующих МНК)
преобразуется к виду:

Отсюда параметры
уравнения рассчитываются по формулам:

Интерпретация
параметров линейного уравнения тренда
:

—
уровень ряда за период времени t^у=0;

—
средний абсолютный прирост уровня ряда
за единичный промежуток времени.

В нашем примере
четное число уровней ряда: n=12. Следовательно,
условная переменная времени для 6-ого
элемента ряда будет равна –1, а для 7-ого
+1. Значения переменной _i^y
содержатся во 2-ом столбце таблицы 9.

Параметры линейного
тренда будут:
=14257,5/572=24,93;=8845/12=737,08.
Это значит, что с каждым кварталом объем
выпуска товара в среднем увеличивается
на 2∙28,7 усл.ед. А средний за период с
1993 по 1995гг объем выпуска составил 738,75
усл.ед.

Рассчитаем значения
трендовой компоненты по формуле
(столбец 7
таблицы 9).

5) Учет сезонной
составляющей в выровненных уровнях
ряда (=T+S).
Результаты расчета для нашего примера
представлены в столбце 8 таблицы 9.

6) Расчет абсолютной
ошибки
временного ряда (Е=x—)
осуществляется для оценки качества
полученной модели. Результаты расчета
для нашего примера представлены в
столбце 9 таблицы 9.

Таблица 9

T	t^у	x	x^c	x- x^c	x^s	T		E
1	2	3	4	5	6	7	8	9
1	-11	410	—	—	477,15	462,90	395,75	14,25
2	-9	560	561,67	-1,67	561,60	512,75	511,15	48,85
3	-7	715	591,67	123,33	551,60	562,60	726,00	-11,01
4	-5	500	578,33	-78,33	594,65	612,45	517,80	-17,80
5	-3	520	586,67	-66,67	587,15	662,31	595,15	-75,15
6	-1	740	745,00	-5,00	741,60	712,16	710,56	29,44
7	1	975	795,00	180,00	811,60	762,00	925,41	49,59
8	3	670	783,33	-113,33	764,65	811,86	717,21	-47,21
9	5	705	775,00	-70,00	772,15	861,71	794,56	-89,56
10	7	950	951,67	-1,67	951,60	911,56	909,97	40,03
11	9	1200	1016,67	183,33	1036,60	961,41	1124,82	75,18
12	11	900	—	—	994,65	1011,27	916,61	-16,61
Итого	8845			8845,00	8845,00	8845,00	16,61

Значимость
параметров линейного уравнения тренда
(Т)
определяется на основе t-критерия
Стьюдента также как и в линейном парном
регрессионном анализе.

Прогнозирование
по аддитивной модели.

Пусть требуется
дать прогноз уровня временного ряда на
период (n+1).
Точечный прогноз значения уровня
временного ряда х_n₊₁
в аддитивной модели есть сумма трендовой
компоненты и сезонной компоненты
(соответствующей i–ому
сезону прогноза):
=T_n₊₁+S_i.

Для построения
доверительного интервала прогноза
нужно рассчитать среднюю ошибку прогноза:

^р
=,

где h—
число параметров в уравнении тренда;

t^yp
– значение условной переменной времени
для периода прогнозирования.

Затем рассчитаем
предельную ошибку прогноза: ^р
=t_a·_р,

где t_a—
коэффициент доверия, определяемый по
таблицам Стьюдента по уровню значимости
α и числу степеней свободы равным (n—h).

Окончательно
получим: (-^р;

+^р).

Задание №
4

На
основе данных, приведенных в таблице
10 и соответствующих Вашему варианту
(таблица 11), постройте модель временного
ряда. Для этого требуется:

Построить
коррелограмму и определить имеет ли
ряд тенденцию и сезонные колебания.
Провести сглаживание
ряда скользящей средней и рассчитать
значения сезонной составляющей.
Построить уравнения
тренда и сделать выводы.
На основе полученной
модели сделать прогноз на следующие
два квартала с учетом выявленной
сезонности.

Таблица 10

Основные показатели
развития производственной фирмы

за период с 2002 по
2007 гг. (по сопоставимой оценке)

N наблюдения	Год	Квартал	Объем производства продукции, млн.руб.	Среднегодовая стоимость основных производственных фондов, млн.руб.	Среднегодовая численность ППП, чел.,	Дебиторская задолженность, млн.руб.	Среднегодовая стоимость оборотных средств, млн.руб.	Балансовая прибыль, млн. руб.	Чистая прибыль, млн. руб.
А	Б	В	1	2	3	4	5	6	7
1	2002	1	1065	1062	713	25	837	94	36
2	2	851	682	507	27	685	78	27
3	3	531	726	361	34	837	87	22
4		4	922	1153	557	44	1161	75	29
5	2003	1	1095	1213	607	42	1151	84	34
6	2	986	898	598	39	822	63	28
7	3	822	794	368	48	1383	86	30
8		4	1137	1441	646	60	884	82	35
9	2004	1	1301	1600	693	63	1309	78	40
10	2	1038	967	718	40	1028	72	33
11	3	780	1246	363	48	1771	84	33
12	4	1435	1458	639	71	1310	102	40
13	2005	1	1593	1412	708	87	1372	112	36
14	2	1658	891	614	65	1272	92	27
15	3	1363	1061	348	67	1821	99	30
16	4	1737	1287	636	76	1571	113	36
17	2006	1	1719	1635	825	101	1758	95	36
18		2	1521	1166	622	84	1505	79	28
19		3	1049	1230	514	73	2109	112	28
20	4	1790	1514	703	93	1787	116	28
21	2007	1	2016	1642	797	96	2197	90	39

Таблица 11

Номера наблюдений
и показатель,

соответствующие
варианту контрольной работы

Номер варианта	Номер начального наблюдения	Номер конечного наблюдения	Номер показателя из табл.4	Номер варианта	Номер начального наблюдения	Номер конечного наблюдения	Номер показателя из табл.4
1	2	3	4	5	6	7	8
1	1	12	1	51	1	12	6
2	2	13	1	52	2	13	6
3	3	14	1	53	3	14	6
4	4	15	1	54	4	15	6
5	5	16	1	55	5	16	6
6	6	17	1	56	6	17	6
7	7	18	1	57	7	18	6
8	8	19	1	58	8	19	6
9	9	20	1	59	9	20	6
10	10	21	1	60	10	21	6
11	1	12	2	61	1	12	7
12	2	13	2	62	2	13	7
13	3	14	2	63	3	14	7
14	4	15	2	64	4	15	7
15	5	16	2	65	5	16	7
16	6	17	2	66	6	17	7
17	7	18	2	67	7	18	7
18	8	19	2	68	8	19	7
19	9	20	2	69	9	20	7
20	10	21	2	70	10	21	7
21	1	12	3	71	1	12	1
22	2	13	3	72	2	13	2
23	3	14	3	73	3	14	3
24	4	15	3	74	4	15	4
25	5	16	3	75	5	16	5
26	6	17	3	76	6	17	6
27	7	18	3	77	7	18	7
28	8	19	3	78	8	19	1
29	9	20	3	79	9	20	2
30	10	21	3	80	10	21	3
31	1	12	4	81	1	12	4
32	2	13	4	82	2	13	5
33	3	14	4	83	3	14	6
34	4	15	4	84	4	15	7
35	5	16	4	85	5	16	1
36	6	17	4	86	6	17	2

Окончание табл.
11

1	2	3	4	5	6	7	8
37	7	18	4	87	7	18	3
38	8	19	4	88	8	19	4
39	9	20	4	89	9	20	5
40	10	21	4	90	10	21	6
41	1	12	5	91	1	12	7
42	2	13	5	92	2	13	1
43	3	14	5	93	3	14	2
44	4	15	5	94	4	15	3
45	5	16	5	95	5	16	4
46	6	17	5	96	6	17	5
47	7	18	5	97	7	18	6
48	8	19	5	98	8	19	7
49	9	20	5	99	9	20	1
50	10	21	5	0	10	21	2

Соседние файлы в папке Эконометрика (Задание)

#
13.02.2015139.26 Кб67_318_эконометрика.xls
#
#
#
13.02.201512.54 Mб12Практикум по эконометрике — Елисеева И.И..pdf

Источник

Анализ временных рядов

Общее введение
Две основные цели
Идентификация модели
временных рядов
- Систематическая
  составляющая и случайный шум
- Два общих типа компонент
  временных рядов
- Анализ тренда
- Анализ сезонности
АРПСС (Бокс и Дженкинс) и
автокорреляции
- Общее введение
- Два основных процесса
- Модель АРПСС
- Идентификация
- Оценивание параметров
- Оценивание модели
Прерванные временные ряды
Экспоненциальное
сглаживание
- Общее введение
- Простое экспоненциальное
  сглаживание
- Выбор лучшего значения
  параметра a (альфа)
- Индексы качества подгонки
- Сезонная и несезонная модели
  с трендом или без тренда
Сезонная декомпозиция (метод
Census I)
- Общее введение
- Вычисления
Сезонная корректировка X-11
(метод Census II)
- Сезонная корректировка:
  основные идеи и термины
- Метод Census II
- Таблицы результатов
  корректировки X-11
- Подробное описание всех
  таблиц результатов, вычисляемых в методе X-11
Анализ распределенных лагов
- Общая цель
- Общая модель
- Распределенный лаг Алмона
Одномерный анализ Фурье
Кросс-спектральный анализ
- Общее введение
- Основные понятия и принципы
- Результаты для каждой
  переменной
- Кросс-периодограмма,
  кросс-плотность, квадратурная плотность и
  кросс-амплитуда
- Квадрат когерентности,
  усиление и фазовый сдвиг
- Как создавались данные для
  примера
Спектральный анализ —
Основные понятия и принципы
- Частота и период
- Общая структура модели
- Простой пример
- Периодограмма
- Проблема рассеяния
- Добавление констант во
  временной ряд (пэддинг)
- Косинус-сглаживание
- Окна данных и оценки
  спектральной плотности
- Подготовка данных к анализу
- Результаты для случая, когда в
  ряде отсутствует периодичность
Быстрое преобразование Фурье
- Общее введение
- Вычисление БПФ во временных
  рядах

В следующих разделах мы вначале представим
обзор методов, используемых для идентификации
моделей временных рядов (таких как сглаживание,
подгонка и автокорреляции). Затем опишем общий
класс моделей, которые могут быть использованы
для описания рядов и построения прогнозов
(модели авторегрессии и скользящего среднего).
Наконец, расскажем о некоторых простых, но часто
используемых методах, основанных на линейной
регрессии. За дальнейшей информацией обратитесь
к соответствующим разделам.

Общее введение

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

Подробное обсуждение этих методов можно найти
в следующих работах: Anderson (1976), Бокс и Дженкинс
(1976), Kendall (1984), Kendall and Ord (1990), Montgomery, Johnson, and Gardiner (1990),
Pankratz (1983), Shumway (1988), Vandaele (1983), Walker (1991), Wei (1989).

Две основные цели

Существуют две основные цели анализа временных
рядов: (1) определение природы ряда и (2)
прогнозирование (предсказание будущих значений
временного ряда по настоящим и прошлым
значениям). Обе эти цели требуют, чтобы модель
ряда была идентифицирована и, более или менее,
формально описана. Как только модель определена,
вы можете с ее помощью интерпретировать
рассматриваемые данные (например, использовать в
вашей теории для понимания сезонного изменения
цен на товары, если занимаетесь экономикой). Не
обращая внимания на глубину понимания и
справедливость теории, вы можете
экстраполировать затем ряд на основе найденной
модели, т.е. предсказать его будущие значения.

Идентификация модели временных
рядов

Систематическая
составляющая и случайный шум
Два общих типа компонент
временных рядов
Анализ тренда
Анализ сезонности

За более полной информацией о простых
автокорреляциях (обсуждаемых в этом разделе) и
других автокорреляциях, см. Anderson (1976), Box and Jenkins
(1976), Kendall (1984), Pankratz (1983), and Vandaele (1983). См. также:

АРПСС (Бокс и Дженкинс) и
автокорреляции
Прерванные временные ряды
Экспоненциальное
сглаживание
Сезонная декомпозиция (метод
Census I)
Сезонная корректировка X-11
(метод Census II)
Таблицы результатов
корректировки X-11
Анализ распределенных лагов
Одномерный анализ Фурье
Кросс-спектральный анализ
Основные понятия и принципы
Быстрое преобразование Фурье

Систематическая составляющая и
случайный шум

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

Два общих типа компонент
временных рядов

Большинство регулярных составляющих временных
рядов принадлежит к двум классам: они являются
либо трендом, либо сезонной составляющей. Тренд
представляет собой общую систематическую
линейную или нелинейную компоненту, которая
может изменяться во времени. Сезонная
составляющая — это периодически повторяющаяся
компонента. Оба эти вида регулярных компонент
часто присутствуют в ряде одновременно.
Например, продажи компании могут возрастать из
года в год, но они также содержат сезонную
составляющую (как правило, 25% годовых продаж
приходится на декабрь и только 4% на август).

Эту общую модель можно понять на
«классическом» ряде — Ряд G (Бокс и
Дженкинс, 1976, стр. 531), представляющем месячные
международные авиаперевозки (в тысячах) в
течение 12 лет с 1949 по 1960 (см. файл Series_g.sta).
График месячных перевозок ясно показывает почти
линейный тренд, т.е. имеется устойчивый рост
перевозок из года в год (примерно в 4 раза больше
пассажиров перевезено в 1960 году, чем в 1949). В то же
время характер месячных перевозок повторяется,
они имеют почти один и тот же характер в каждом
годовом периоде (например, перевозок больше в
отпускные периоды, чем в другие месяцы). Этот
пример показывает довольно определенный тип
модели временного ряда, в которой амплитуда
сезонных изменений увеличивается вместе с
трендом. Такого рода модели называются моделями
с мультипликативной сезонностью.

Анализ тренда

Не существует «автоматического» способа
обнаружения тренда в временном ряде. Однако если
тренд является монотонным (устойчиво возрастает
или устойчиво убывает), то анализировать такой
ряд обычно нетрудно. Если временные ряды
содержат значительную ошибку, то первым шагом
выделения тренда является сглаживание.

Сглаживание. Сглаживание всегда включает
некоторый способ локального усреднения данных,
при котором несистематические компоненты
взаимно погашают друг друга. Самый общий метод
сглаживания — скользящее среднее, в котором
каждый член ряда заменяется простым или
взвешенным средним n соседних членов, где n
— ширина «окна» (см. Бокс и Дженкинс, 1976; Velleman
and Hoaglin, 1981). Вместо среднего можно использовать
медиану значений, попавших в окно. Основное
преимущество медианного сглаживания, в
сравнении со сглаживанием скользящим средним,
состоит в том, что результаты становятся более
устойчивыми к выбросам (имеющимся внутри окна).
Таким образом, если в данных имеются выбросы
(связанные, например, с ошибками измерений), то
сглаживание медианой обычно приводит к более
гладким или, по крайней мере, более
«надежным» кривым, по сравнению со
скользящим средним с тем же самым окном. Основной
недостаток медианного сглаживания в том, что при
отсутствии явных выбросов, он приводит к более
«зубчатым» кривым (чем сглаживание
скользящим средним) и не позволяет использовать
веса.

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

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

Анализ сезонности

Периодическая и сезонная зависимость
(сезонность) представляет собой другой общий тип
компонент временного ряда. Это понятие было
проиллюстрировано ранее на примере
авиаперевозок пассажиров. Можно легко видеть,
что каждое наблюдение очень похоже на соседнее;
дополнительно, имеется повторяющаяся сезонная
составляющая, это означает, что каждое
наблюдение также похоже на наблюдение, имевшееся
в том же самом месяце год назад. В общем,
периодическая зависимость может быть формально
определена как корреляционная зависимость
порядка k между каждым i-м элементом ряда и
(i-k)-м элементом (Kendall, 1976). Ее можно измерить с
помощью автокорреляции (т.е. корреляции между
самими членами ряда); k обычно называют лагом
(иногда используют эквивалентные термины:
сдвиг, запаздывание). Если ошибка измерения не
слишком большая, то сезонность можно определить
визуально, рассматривая поведение членов ряда
через каждые k временных единиц.

Автокорреляционная коррелограмма. Сезонные
составляющие временного ряда могут быть найдены
с помощью коррелограммы. Коррелограмма
(автокоррелограмма) показывает численно и
графически автокорреляционную функцию (AКФ),
иными словами коэффициенты автокорреляции (и их
стандартные ошибки) для последовательности
лагов из определенного диапазона (например, от 1
до 30). На коррелограмме обычно отмечается
диапазон в размере двух стандартных ошибок на
каждом лаге, однако обычно величина
автокорреляции более интересна, чем ее
надежность, потому что интерес в основном
представляют очень сильные (а, следовательно,
высоко значимые) автокорреляции (см. Элементарные
понятия статистики).

Исследование коррелограмм. При изучении
коррелограмм следует помнить, что
автокорреляции последовательных лагов
формально зависимы между собой. Рассмотрим
следующий пример. Если первый член ряда тесно
связан со вторым, а второй с третьим, то первый
элемент должен также каким-то образом зависеть
от третьего и т.д. Это приводит к тому, что
периодическая зависимость может существенно
измениться после удаления автокорреляций
первого порядка, т.е. после взятия разности с
лагом 1).

Частные автокорреляции. Другой полезный
метод исследования периодичности состоит в
исследовании частной автокорреляционной
функции (ЧАКФ), представляющей собой
углубление понятия обычной автокорреляционной
функции. В ЧАКФ устраняется зависимость между
промежуточными наблюдениями (наблюдениями внутри
лага). Другими словами, частная автокорреляция на
данном лаге аналогична обычной автокорреляции,
за исключением того, что при вычислении из нее
удаляется влияние автокорреляций с меньшими
лагами (см. Бокс и Дженкинс, 1976; см. также McDowall,
McCleary, Meidinger, and Hay, 1980). На лаге 1 (когда нет
промежуточных элементов внутри лага), частная
автокорреляция равна, очевидно, обычной
автокорреляции. На самом деле, частная
автокорреляция дает более «чистую» картину
периодических зависимостей.

Удаление периодической зависимости. Как
отмечалось выше, периодическая составляющая для
данного лага k может быть удалена взятием
разности соответствующего порядка. Это означает,
что из каждого i-го элемента ряда вычитается (i-k)-й
элемент. Имеются два довода в пользу таких
преобразований.

Во-первых, таким образом можно определить
скрытые периодические составляющие ряда.
Напомним, что автокорреляции на
последовательных лагах зависимы. Поэтому
удаление некоторых автокорреляций изменит
другие автокорреляции, которые, возможно,
подавляли их, и сделает некоторые другие
сезонные составляющие более заметными.

Во-вторых, удаление сезонных составляющих
делает ряд стационарным,
что необходимо для применения АРПСС
и других методов, например, спектрального
анализа.

АРПСС

Общее введение
Два основных процесса
Модель АРПСС
Идентификация
Оценивание параметров
Оценивание модели

Дополнительная информация о методах Анализа
временных рядов дана также в следующих
разделах:

Идентификация модели
временных рядов
Прерванные временные ряды
Экспоненциальное
сглаживание
Сезонная декомпозиция (метод
Census I)
Сезонная корректировка X-11
(метод Census II)
Таблицы результатов
корректировки X-11
Анализ распределенных лагов
Одномерный анализ Фурье
Кросс-спектральный анализ
Основные понятия и принципы
Быстрое преобразование Фурье

Общее введение

Процедуры оценки параметров и прогнозирования,
описанные в разделе Идентификация
модели временных рядов, предполагают, что
математическая модель процесса известна. В
реальных данных часто нет отчетливо выраженных
регулярных составляющих. Отдельные наблюдения
содержат значительную ошибку, тогда как вы
хотите не только выделить регулярные компоненты,
но также построить прогноз. Методология АРПСС,
разработанная Боксом и Дженкинсом (1976), позволяет
это сделать. Данный метод чрезвычайно популярен
во многих приложениях, и практика подтвердила
его мощность и гибкость (Hoff, 1983; Pankratz, 1983; Vandaele, 1983).
Однако из-за мощности и гибкости, АРПСС — сложный
метод. Его не так просто использовать, и
требуется большая практика, чтобы овладеть им.
Хотя часто он дает удовлетворительные
результаты, они зависят от квалификации
пользователя (Bails and Peppers, 1982). Следующие разделы
познакомят вас с его основными идеями. Для
интересующихся кратким, рассчитанным на
применение, (нематематическим) введением в АРПСС,
рекомендуем книгу McCleary, Meidinger, and Hay (1980).

Два основных процесса

Процесс авторегрессии. Большинство
временных рядов содержат элементы, которые
последовательно зависят друг от друга. Такую
зависимость можно выразить следующим
уравнением:

x_t =
+ ₁*x_(t-1) + ₂*x_(t-2) + ₃*x_(t-3) + … +

Здесь:
—
константа (свободный член),
₁,
₂,
₃
— параметры авторегрессии.

Вы видите, что каждое наблюдение есть сумма
случайной компоненты (случайное воздействие, ) и линейной
комбинации предыдущих наблюдений.

Требование стационарности. Заметим, что
процесс авторегрессии будет стационарным
только, если его параметры лежат в определенном
диапазоне. Например, если имеется только один
параметр, то он должен находиться в интервале -1<<+1. В противном случае,
предыдущие значения будут накапливаться и
значения последующих x_t могут быть
неограниченными, следовательно, ряд не будет стационарным.
Если имеется несколько параметров
авторегрессии, то можно определить аналогичные
условия, обеспечивающие стационарность (см.
например, Бокс и Дженкинс, 1976; Montgomery, 1990).

Процесс скользящего среднего. В отличие от
процесса авторегрессии, в процессе скользящего
среднего каждый элемент ряда подвержен
суммарному воздействию предыдущих ошибок. В
общем виде это можно записать следующим образом:

x_t = µ + _t — ₁*_(t-1) — ₂*_(t-2) — ₃*_(t-3) — …

Здесь:
µ —
константа,
₁,
₂,
₃ —
параметры скользящего среднего.

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

Обратимость. Не вдаваясь в детали, отметим,
что существует «двойственность» между
процессами скользящего среднего и авторегрессии
(см. например, Бокс и Дженкинс, 1976; Montgomery, Johnson, and
Gardiner, 1990). Это означает, что приведенное выше
уравнение скользящего среднего можно переписать
(обратить) в виде уравнения авторегрессии
(неограниченного порядка), и наоборот. Это так
называемое свойство обратимости. Имеются
условия, аналогичные приведенным выше условиям стационарности,
обеспечивающие обратимость модели.

Модель АРПСС

Модель авторегрессии и скользящего среднего. Общая
модель, предложенная Боксом и Дженкинсом (1976)
включает как параметры авторегрессии, так и
параметры скользящего среднего. Именно, имеется
три типа параметров модели: параметры
авторегрессии (p), порядок разности (d), параметры
скользящего среднего (q). В обозначениях Бокса
и Дженкинса модель записывается как АРПСС (p, d, q).
Например, модель (0, 1, 2) содержит 0
(нуль) параметров авторегрессии (p) и 2
параметра скользящего среднего (q), которые
вычисляются для ряда после взятия разности с
лагом 1.

Идентификация. Как отмечено ранее, для
модели АРПСС необходимо, чтобы ряд был стационарным,
это означает, что его среднее постоянно, а
выборочные дисперсия и автокорреляция не
меняются во времени. Поэтому обычно необходимо
брать разности ряда до тех пор, пока он не станет
стационарным
(часто также применяют логарифмическое
преобразование для стабилизации дисперсии).
Число разностей, которые были взяты, чтобы
достичь стационарности, определяются параметром
d (см. предыдущий раздел). Для того чтобы
определить необходимый порядок разности, нужно
исследовать график ряда и автокоррелограмму.
Сильные изменения уровня (сильные скачки вверх
или вниз) обычно требуют взятия несезонной
разности первого порядка (лаг=1). Сильные
изменения наклона требуют взятия разности
второго порядка. Сезонная составляющая требует
взятия соответствующей сезонной разности (см.
ниже). Если имеется медленное убывание
выборочных коэффициентов автокорреляции в
зависимости от лага, обычно берут разность
первого порядка. Однако следует помнить, что для
некоторых временных рядов нужно брать разности
небольшого порядка или вовсе не брать их.
Заметим, что чрезмерное количество взятых
разностей приводит к менее стабильным оценкам
коэффициентов.

На этом этапе (который обычно называют идентификацией
порядка модели, см. ниже) вы также должны
решить, как много параметров авторегрессии (p)
и скользящего среднего (q) должно
присутствовать в эффективной и экономной модели
процесса. (Экономность модели означает, что в
ней имеется наименьшее число параметров и
наибольшее число степеней свободы среди всех
моделей, которые подгоняются к данным). На
практике очень редко бывает, что число
параметров p или q больше 2 (см. ниже более
полное обсуждение).

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

Константа в моделях АРПСС. Дополнительно
модели АРПСС могут содержать константу,
интерпретация которой зависит от подгоняемой
модели. Именно, если (1) в модели нет параметров
авторегрессии, то константа есть среднее значение ряда, если (2)
параметры авторегрессии имеются, то константа
представляет собой свободный член. Если бралась
разность ряда, то константа представляет собой
среднее или свободный член преобразованного
ряда. Например, если бралась первая разность
(разность первого порядка), а параметров
авторегрессии в модели нет, то константа
представляет собой среднее значение
преобразованного ряда и, следовательно, коэффициент
наклона линейного тренда исходного.

Идентификация

Число оцениваемых параметров. Конечно, до
того, как начать оценивание, вам необходимо
решить, какой тип модели будет подбираться к
данным, и какое количество параметров
присутствует в модели, иными словами, нужно
идентифицировать модель АРПСС. Основными
инструментами идентификации порядка модели
являются графики, автокорреляционная функция
(АКФ), частная автокорреляционная функция (ЧАКФ).
Это решение не является простым и требуется
основательно поэкспериментировать с
альтернативными моделями. Тем не менее,
большинство встречающихся на практике временных
рядов можно с достаточной степенью точности
аппроксимировать одной из 5 основных моделей (см.
ниже), которые можно идентифицировать по виду
автокорреляционной (АКФ) и частной
автокорреляционной функции (ЧАКФ). Ниже дается
список этих моделей, основанный на рекомендациях
Pankratz (1983); дополнительные практические советы
даны в Hoff (1983), McCleary and Hay (1980), McDowall, McCleary, Meidinger, and Hay
(1980), and Vandaele (1983). Отметим, что число параметров
каждого вида невелико (меньше 2), поэтому нетрудно
проверить альтернативные модели.

Один параметр (p): АКФ — экспоненциально
убывает; ЧАКФ — имеет резко выделяющееся значение
для лага 1, нет корреляций на других лагах.
Два параметра авторегрессии (p): АКФ имеет
форму синусоиды или экспоненциально убывает;
ЧАКФ имеет резко выделяющиеся значения на лагах 1,
2, нет корреляций на других лагах.
Один параметр скользящего среднего (q): АКФ
имеет резко выделяющееся значение на лаге 1,
нет корреляций на других лагах. ЧАКФ
экспоненциально убывает.
Два параметра скользящего среднего (q): АКФ
имеет резко выделяющиеся значения на лагах 1, 2,
нет корреляций на других лагах. ЧАКФ имеет форму
синусоиды или экспоненциально убывает.
Один параметр авторегрессии (p) и один параметр
скользящего среднего (q): АКФ экспоненциально
убывает с лага 1; ЧАКФ — экспоненциально
убывает с лага 1.

Сезонные модели. Мультипликативная сезонная
АРПСС представляет естественное развитие и
обобщение обычной модели АРПСС на ряды, в которых
имеется периодическая сезонная компонента. В
дополнении к несезонным параметрам, в модель
вводятся сезонные параметры для определенного
лага (устанавливаемого на этапе идентификации
порядка модели). Аналогично параметрам простой
модели АРПСС, эти параметры называются: сезонная
авторегрессия (ps), сезонная разность (ds) и
сезонное скользящее среднее (qs). Таким
образом, полная сезонная АРПСС может быть
записана как АРПСС (p,d,q)(ps,ds,qs).
Например, модель (0,1,2)(0,1,1) включает 0
регулярных параметров авторегрессии, 2
регулярных параметра скользящего среднего и 1
параметр сезонного скользящего среднего. Эти
параметры вычисляются для рядов, получаемых
после взятия одной разности с лагом 1 и далее
сезонной разности. Сезонный лаг, используемый
для сезонных параметров, определяется на этапе
идентификации порядка модели.

Общие рекомендации относительно выбора
обычных параметров (с помощью АКФ и ЧАКФ)
полностью применимы к сезонным моделям. Основное
отличие состоит в том, что в сезонных рядах АКФ и
ЧАКФ имеют существенные значения на лагах,
кратных сезонному лагу (в дополнении к
характерному поведению этих функций,
описывающих регулярную (несезонную) компоненту
АРПСС).

Оценивание параметров

Существуют различные методы оценивания
параметров, которые дают очень похожие оценки, но
для данной модели одни оценки могут быть более
эффективны, а другие менее эффективны. В общем, во
время оценивания порядка модели используется
так называемый квазиньютоновский алгоритм
максимизации правдоподобия (вероятности)
наблюдения значений ряда по значениям
параметров (см. Нелинейное
оценивание). Практически это требует
вычисления (условных) сумм квадратов (SS)
остатков модели. Имеются различные способы
вычисления суммы квадратов остатков SS; вы
можете выбрать: (1) приближенный метод
максимального правдоподобия МакЛеода и Сейлза
(1983), (2) приближенный метод максимального
правдоподобия с итерациями назад, (3)точный метод
максимального правдоподобия по Meларду (1984).

Сравнение методов. В общем, все методы дают
очень похожие результаты. Также все методы
показали примерно одинаковую эффективность на
реальных данных. Однако метод 1 (см. выше) —
самый быстрый, и им можно пользоваться для
исследования очень длинных рядов (например,
содержащих более 30,000 наблюдений). Метод Меларда
(номер 3) может оказаться неэффективным, если
оцениваются параметры сезонной модели с большим
сезонным лагом (например, 365 дней). С другой
стороны, вы можете использовать вначале
приближенный метод максимального правдоподобия
(для того, чтобы найти прикидочные оценки
параметров), а затем точный метод; обычно
требуется только несколько итераций точного
метода (номер 3, выше), чтобы получить
окончательные оценки.

Стандартные ошибки оценок. Для всех оценок
параметров вычисляются так называемые асимптотические
стандартные ошибки, для вычисления которых
используется матрица частных производных
второго порядка, аппроксимируемая конечными
разностями (см. также раздел Нелинейное
оценивание).

Штраф. Процедура оценивания минимизирует
(условную) сумму квадратов остатков модели. Если
модель не является адекватной, может случиться
так, что оценки параметров на каком-то шаге
станут неприемлемыми — очень большими (например,
не удовлетворяют условию стационарности). В
таком случае, SS будет приписано очень большое
значение (штрафное значение). Обычно это
«заставляет» итерационный процесс удалить
параметры из недопустимой области. Однако в
некоторых случаях и эта стратегия может
оказаться неудачной, и вы все равно увидите на
экране (во время процедуры оценивания) очень
большие значения SS на серии итераций. В таких
случаях следует с осторожностью оценивать
пригодность модели. Если модель содержит много
параметров и, возможно, имеется интервенция (см.
ниже), то следует несколько раз испытать процесс
оценивания с различными начальными. Если модель
содержит много параметров и, возможно,
интервенцию (см. ниже), вам следует повторить
процедуру с различными начальными значениями
параметров.

Оценивание модели

Оценки параметров. Если значения
вычисляемой t статистики не значимы,
соответствующие параметры в большинстве случаев
удаляются из модели без ущерба подгонки.

Другой критерий качества. Другой обычной
мерой надежности модели является сравнение
прогноза, построенного по урезанному ряду с
«известными (исходными) данными».

Однако качественная модель должна не только
давать достаточно точный прогноз, но быть
экономной и иметь независимые остатки,
содержащие только шум без систематических
компонент (в частности, АКФ остатков не должна
иметь какой-либо периодичности). Поэтому
необходим всесторонний анализ остатков. Хорошей
проверкой модели являются: (a) график остатков и
изучение их трендов, (b) проверка АКФ остатков (на
графике АКФ обычно отчетливо видна
периодичность).

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

Ограничения. Следует напомнить, что модель
АРПСС является подходящей только для рядов,
которые являются стационарными
(среднее, дисперсия и автокорреляция примерно
постоянны во времени); для нестационарных рядов
следует брать разности. Рекомендуется иметь, как
минимум, 50 наблюдений в файле исходных данных.
Также предполагается, что параметры модели
постоянны, т.е. не меняются во времени.

Прерванные временные ряды

Обычный вопрос, возникающий при анализе
временных рядов, состоит в следующем,
воздействует или нет внешнее событие на
последовательность наблюдений. Например,
привела ли новая экономическая политика к росту
экономики, как обещалось; изменил ли новый закон
интенсивность преступлений и т.д. В общем, нужно
оценивать воздействия одного или нескольких
дискретных событий на значения ряда. Этот вид
анализа прерванных временных рядов подробно
описан в книге McDowall, McCleary, Meidinger, and Hay (1980).
Различают следующие три типа воздействий: (1)
устойчивое скачкообразное, (2) устойчивое
постепенное, (3) скачкообразное временное. См.
также следующие разделы:

Идентификация модели
временных рядов
АРПСС
Экспоненциальное
сглаживание
Сезонная декомпозиция (метод
Census I)
Сезонная корректировка X-11
(метод Census II)
Таблицы результатов
корректировки X-11
Анализ распределенных лагов
Одномерный анализ Фурье
Кросс-спектральный анализ
Основные понятия и принципы
Быстрое преобразование Фурье

Экспоненциальное сглаживание

Общее введение
Простое экспоненциальное
сглаживание
Выбор лучшего значения
параметра a (альфа)
Индексы качества подгонки
Сезонная и несезонная модели
с трендом или без тренда

См. также:

Идентификация модели
временных рядов
АРПСС (Бокс и Дженкинс) и
автокорреляции
Прерванные временные ряды
Сезонная декомпозиция (метод
Census I)
Сезонная корректировка X-11
(метод Census II)
Таблицы результатов
корректировки X-11
Анализ распределенных лагов
Одномерный анализ Фурье
Кросс-спектральный анализ
Основные понятия и принципы
Быстрое преобразование Фурье

Общее введение

Экспоненциальное сглаживание — это очень
популярный метод прогнозирования многих
временных рядов. Исторически метод был
независимо открыт Броуном и Холтом. Броун служил
на флоте США во время второй мировой войны, где
занимался обнаружением подводных лодок и
системами наведения. Позже он применил открытый
им метод для прогнозирования спроса на запасные
части. Свои идеи он описал в книге, вышедшей в
свет в 1959 году. Исследования Холта были
поддержаны Департаментом военно-морского флота
США. Независимо друг от друга, Броун и Холт
открыли экспоненциальное сглаживание для
процессов с постоянным трендом, с линейным
трендом и для рядов с сезонной составляющей.

Gardner (1985), предложил «единую» классификацию
методов экспоненциального сглаживания.
Превосходное введение в эти методы можно найти в
книгах Makridakis, Wheelwright, and McGee (1983), Makridakis and Wheelwright (1989),
Montgomery, Johnson, and Gardiner (1990).

Простое экспоненциальное
сглаживание

Простая и прагматически ясная модель
временного ряда имеет следующий вид: X_t = b + _t, где b —
константа и
(эпсилон) — случайная ошибка. Константа b относительно
стабильна на каждом временном интервале, но
может также медленно изменяться со временем.
Один из интуитивно ясных способов выделения b
состоит в том, чтобы использовать сглаживание
скользящим средним, в котором последним
наблюдениям приписываются большие веса, чем
предпоследним, предпоследним большие веса, чем
пред-предпоследним и т.д. Простое
экспоненциальное именно так и устроено. Здесь
более старым наблюдениям приписываются
экспоненциально убывающие веса, при этом, в
отличие от скользящего среднего, учитываются все
предшествующие наблюдения ряда, а не те, что
попали в определенное окно. Точная формула
простого экспоненциального сглаживания имеет
следующий вид:

S_t = *X_t + (1-)*S_t-1

Когда эта формула применяется рекурсивно, то
каждое новое сглаженное значение (которое
является также прогнозом) вычисляется как
взвешенное среднее текущего наблюдения и
сглаженного ряда. Очевидно, результат
сглаживания зависит от параметра (альфа). Если равно 1, то
предыдущие наблюдения полностью игнорируются.
Если равно 0, то
игнорируются текущие наблюдения. Значения между 0, 1 дают
промежуточные результаты.

Эмпирические исследования Makridakis и др. (1982;
Makridakis, 1983) показали, что весьма часто простое
экспоненциальное сглаживание дает достаточно
точный прогноз.

Выбор лучшего значения
параметра (альфа)

Gardner (1985) обсуждает различные теоретические и
эмпирические аргументы в пользу выбора
определенного параметра сглаживания. Очевидно,
из формулы, приведенной выше, следует, что должно попадать в интервал между 0
(нулем) и 1 (хотя Brenner et al., 1968, для дальнейшего
применения анализа АРПСС считают, что 0<<2). Gardner (1985)
сообщает, что на практике обычно рекомендуется
брать меньше .30.
Однако в исследовании Makridakis et al., (1982), большее .30, часто
дает лучший прогноз. После обзора литературы,
Gardner (1985) приходит к выводу, что лучше оценивать
оптимально по данным
(см. ниже), чем просто «гадать» или
использовать искусственные рекомендации.

Оценивание лучшего значения с помощью данных. На
практике параметр сглаживания часто ищется с поиском
на сетке. Возможные значения параметра
разбиваются сеткой с определенным шагом.
Например, рассматривается сетка значений от = 0.1 до = 0.9, с шагом 0.1.
Затем выбирается ,
для которого сумма квадратов (или средних
квадратов) остатков (наблюдаемые значения минус
прогнозы на шаг вперед) является минимальной.

Индексы качества подгонки

Самый прямой способ оценки прогноза,
полученного на основе определенного значения — построить график
наблюдаемых значений и прогнозов на один шаг
вперед. Этот график включает в себя также остатки
(отложенные на правой оси Y). Из графика ясно
видно, на каких участках прогноз лучше или хуже.

Такая визуальная проверка точности прогноза
часто дает наилучшие результаты. Имеются также
другие меры ошибки, которые можно использовать
для определения оптимального параметра (см. Makridakis, Wheelwright, and McGee,
1983):

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

Средняя абсолютная ошибка. Средняя
абсолютная ошибка (САО) вычисляется как среднее абсолютных
ошибок. Если она равна 0 (нулю), то имеем
совершенную подгонку (прогноз). В сравнении со
средней квадратической ошибкой, эта мера
«не придает слишком большого значения»
выбросам.

Сумма квадратов ошибок (SSE),
среднеквадратическая ошибка. Эти величины
вычисляются как сумма (или среднее) квадратов
ошибок. Это наиболее часто используемые индексы
качества подгонки.

Относительная ошибка (ОО). Во всех
предыдущих мерах использовались действительные
значения ошибок. Представляется естественным
выразить индексы качества подгонки в терминах относительных
ошибок. Например, при прогнозе месячных продаж,
которые могут сильно флуктуировать (например, по
сезонам) из месяца в месяц, вы можете быть вполне
удовлетворены прогнозом, если он имеет точность
?10%. Иными словами, при прогнозировании
абсолютная ошибка может быть не так интересна
как относительная. Чтобы учесть относительную
ошибку, было предложено несколько различных
индексов (см. Makridakis, Wheelwright, and McGee, 1983). В первом
относительная ошибка вычисляется как:

ОО_t = 100*(X_t — F_t )/X_t

где X_t — наблюдаемое
значение в момент времени t, и F_t — прогноз (сглаженное
значение).

Средняя относительная ошибка (СОО).
Это значение вычисляется как среднее
относительных ошибок.

Средняя абсолютная относительная ошибка
(САОО). Как и в случае с обычной средней
ошибкой отрицательные и положительные
относительные ошибки будут подавлять друг друга.
Поэтому для оценки качества подгонки в целом (для
всего ряда) лучше использовать среднюю абсолютную
относительную ошибку. Часто эта мера более
выразительная, чем среднеквадратическая ошибка.
Например, знание того, что точность прогноза ±5%,
полезно само по себе, в то время как значение 30.8
для средней квадратической ошибки не может быть
так просто проинтерпретировано.

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

Первое сглаженное значение S₀.
Если вы взгляните снова на формулу простого
экспоненциального сглаживания, то увидите, что
следует иметь значение S₀ для
вычисления первого сглаженного значения
(прогноза). В зависимости от выбора параметра (в частности, если близко к 0), начальное
значение сглаженного процесса может оказать
существенное воздействие на прогноз для многих
последующих наблюдений. Как и в других
рекомендациях по применению экспоненциального
сглаживания, рекомендуется брать начальное
значение, дающее наилучший прогноз. С другой
стороны, влияние выбора уменьшается с длиной
ряда и становится некритичным при большом числе
наблюдений.

Сезонная и несезонная модели с
трендом или без тренда

В дополнение к простому экспоненциальному
сглаживанию, были предложены более сложные
модели, включающие сезонную компоненту и
трендом. Общая идея таких моделей состоит в том,
что прогнозы вычисляются не только по предыдущим
наблюдениям (как в простом экспоненциальном
сглаживании), но и с некоторыми задержками, что
позволяет независимо оценить тренд и сезонную
составляющую. Gardner (1985) обсудил различные модели в
терминах сезонности (отсутствует, аддитивная
сезонность, мультипликативная) и тренда
(отсутствует, линейный тренд, экспоненциальный,
демпфированный).

Аддитивная и мультипликативная
сезонность. Многие временные ряды имеют
сезонные компоненты. Например, продажи игрушек
имеют пики в ноябре, декабре и, возможно, летом,
когда дети находятся на отдыхе. Эта
периодичность имеет место каждый год. Однако
относительный размер продаж может слегка
изменяться из года в год. Таким образом, имеет
смысл независимо экспоненциально сгладить
сезонную компоненту с дополнительным
параметром, обычно обозначаемым как (дельта). Сезонные
компоненты, по природе своей, могут быть
аддитивными или мультипликативными. Например, в
течение декабря продажи определенного вида
игрушек увеличиваются на 1 миллион долларов
каждый год. Для того чтобы учесть сезонное
колебание, вы можете добавить в прогноз на каждый
декабрь 1 миллион долларов (сверх
соответствующего годового среднего). В этом
случае сезонность — аддитивная. Альтернативно,
пусть в декабре продажи увеличились на 40%, т.е. в 1.4
раза. Тогда, если общие продажи малы, то
абсолютное (в долларах) увеличение продаж в
декабре тоже относительно мало (процент роста
константа). Если в целом продажи большие, то
абсолютное (в долларах) увеличение продаж будет
пропорционально больше. Снова, в этом случае
продажи увеличатся в определенное число раз, и
сезонность будет мультипликативной (в данном
случае мультипликативная сезонная составляющая
была бы равна 1.4). На графике различие между двумя
видами сезонности состоит в том, что в аддитивной
модели сезонные флуктуации не зависят от
значений ряда, тогда как в мультипликативной
модели величина сезонных флуктуаций зависит от
значений временного ряда.

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

Аддитивная модель:

Прогноз_t = S_t + I_t-p

Мультипликативная модель:

Прогноз_t = S_t*I_t-p

В этой формуле S_t
обозначает (простое) экспоненциально сглаженное
значение ряда в момент t, и I_t-p обозначает сглаженный
сезонный фактор в момент t
минус p (p —
длина сезона). Таким образом, в сравнении с
простым экспоненциальным сглаживанием, прогноз
«улучшается» добавлением или умножением
сезонной компоненты. Эта компонента оценивается
независимо с помощью простого экспоненциального
сглаживания следующим образом:

Аддитивная модель:

I_t = I_t-p + *(1-)*e_t

Мультипликативная модель:

I_t = I_t-p + *(1-)*e_t/S_t

Обратите внимание, что предсказанная сезонная
компонента в момент t
вычисляется, как соответствующая компонента на
последнем сезонном цикле плюс ошибка (e_t,
наблюдаемое минус прогнозируемое значение в
момент t). Ясно, что параметр принимает значения
между 0 и 1. Если он равен нулю, то сезонная
составляющая на следующем цикле та же, что и на
предыдущем. Если
равен 1, то сезонная составляющая
«максимально» меняется на каждом шаге из-за
соответствующей ошибки (множитель (1-) не
рассматривается из-за краткости введения). В
большинстве случаев, когда сезонность
присутствует, оптимальное значение лежит между 0 и 1.

Линейный, экспоненциальный,
демпфированный тренд. Возвращаясь к примеру
с игрушками, мы можем увидеть наличие линейного
тренда (например, каждый год продажи
увеличивались на 1 миллион), экспоненциального
(например, каждый год продажи возрастают в 1.3
раза) или демпфированного тренда (в первом году
продажи возросли на 1 миллион долларов; во втором
увеличение составило только 80% по сравнению с
предыдущим, т.е. на $800,000; в следующем году вновь
увеличение было только на 80%, т.е. на $800,000 * .8 = $640,000
и т.д.). Каждый тип тренда по-своему проявляется в
данных. В целом изменение тренда — медленное в
течение времени, и опять (как и сезонную
компоненту) имеет смысл экспоненциально
сгладить его с отдельным параметром
[обозначаемым (гамма)
— для линейного и экспоненциального тренда, (фи) — для
демпфированного тренда].

Параметры сглаживания (линейный и экспоненциальный тренд) и (демпфированный тренд). Аналогично
сезонной компоненте компонента тренда
включается в процесс экспоненциального
сглаживания. Сглаживание ее производится в
каждый момент времени независимо от других
компонент с соответствующими параметрами. Если равно 0, то тренд
постоянен для всех значений временного ряда (и
для всех прогнозов). Если равно 1, то тренд «максимально»
определяется ошибками наблюдений. Параметр учитывает, как сильно
изменяется тренд, т.е. как быстро он
«демпфируется» или, наоборот, возрастает.

Сезонная декомпозиция (метод Census I)

Общее введение
Вычисления

См. также:

Идентификация модели
временных рядов
АРПСС (Бокс и Дженкинс) и
автокорреляции
Прерванные временные ряды
Экспоненциальное
сглаживание
Сезонная корректировка X-11
(метод Census II)
Таблицы результатов
корректировки X-11
Анализ распределенных лагов
Одномерный анализ Фурье
Кросс-спектральный анализ
Основные понятия и принципы
Быстрое преобразование Фурье

Общее введение

Предположим, что у вас имеются ежемесячные
данные о пассажиропотоке на международных
авиалиниях за 12 лет (см. Бокс и Дженкинс, 1976). Если
изобразить эти данные на графике, то будет хорошо
видно, что (1) объем пассажиропотока имеет во
времени возрастающий линейный тренд, и (2) в ряде
имеется ежегодно повторяющаяся закономерность — сезонность
(большинство перевозок приходится на летние
месяцы, кроме того, имеется пик меньшей высоты в
районе декабрьских каникул). Цель сезонной
декомпозиции и корректировки как раз и состоит в
том, чтобы отделить эти компоненты, то есть
разложить ряд на составляющую тренда, сезонную
компоненту и оставшуюся нерегулярную
составляющую. «Классический» прием,
позволяющий выполнить такую декомпозицию,
известен как метод Census I. Этот метод
описывается и обсуждается в работах Makridakis,
Wheelwright, and McGee (1983) и Makridakis and Wheelwright (1989).

Общая модель. Основная идея сезонной
декомпозиции проста. В общем случае временной
ряд типа того, который описан выше, можно
представить себе состоящим из четырех различных
компонент: (1) сезонной компоненты (обозначается S_t,
где t обозначает момент времени), (2) тренда (T_t),
(3) циклической компоненты (C_t) и (4)
случайной, нерегулярной компоненты или
флуктуации (I_t). Разница между
циклической и сезонной компонентой состоит в
том, что последняя имеет регулярную (сезонную)
периодичность, тогда как циклические факторы
обычно имеют более длительный эффект, который к
тому же меняется от цикла к циклу. В методе Census I
тренд и циклическую компоненту обычно
объединяют в одну тренд-циклическую компоненту
(TC_t). Конкретные функциональные
взаимосвязи между этими компонентами могут
иметь самый разный вид. Однако, можно выделить
два основных способа, с помощью которых они могут
взаимодействовать: аддитивно и мультипликативно:

Аддитивная модель:

X_t = TC_t + S_t + I_t

Мультипликативная модель:

X_t = T_t*C_t*S_t*I_t

Здесь X_t обозначает
значение временного ряда в момент времени t. Если имеются какие-то априорные сведения
о циклических факторах, влияющих на ряд
(например, циклы деловой конъюнктуры), то можно
использовать оценки для различных компонент для
составления прогноза будущих значений ряда.
(Однако для прогнозирования предпочтительнее экспоненциальное
сглаживание, позволяющее учитывать сезонную
составляющую и тренд.)

Аддитивная и мультипликативная
сезонность. Рассмотрим на примере различие
между аддитивной и мультипликативной сезонными
компонентами. График объема продаж детских
игрушек, вероятно, будет иметь ежегодный пик в
ноябре-декабре, и другой — существенно меньший по
высоте — в летние месяцы, приходящийся на
каникулы. Такая сезонная закономерность будет
повторяться каждый год. По своей природе
сезонная компонента может быть аддитивной или
мультипликативной. Так, например, каждый год
объем продаж некоторой конкретной игрушки может
увеличиваться в декабре на 3 миллиона долларов.
Поэтому вы можете учесть эти сезонные изменения, прибавляя
к своему прогнозу на декабрь 3 миллиона. Здесь
мы имеем аддитивную сезонность. Может
получиться иначе. В декабре объем продаж
некоторой игрушки может увеличиваться на 40%, то
есть умножаться на множитель 1.4. Это значит,
например, что если средний объем продаж этой
игрушки невелик, то абсолютное (в денежном
выражении) увеличение этого объема в декабре
также будет относительно небольшим (но в
процентном исчислении оно будет постоянным);
если же игрушка продается хорошо, то и абсолютный
(в долларах) рост объема продаж будет
значительным. Здесь опять, объем продаж
возрастает в число раз, равное определенному множителю,
а сезонная компонента, по своей природе,
мультипликативная компонента (в данном случае
равная 1.4). Если перейти к графикам временных
рядов, то различие между этими двумя видами
сезонности будет проявляться так: в аддитивном
случае ряд будет иметь постоянные сезонные
колебания, величина которых не зависит от общего
уровня значений ряда; в мультипликативном случае
величина сезонных колебаний будет меняться в
зависимости от общего уровня значений ряда.

Аддитивный и мультипликативный тренд-цикл. Рассмотренный
пример можно расширить, чтобы проиллюстрировать
понятия аддитивной и мультипликативной
тренд-циклических компонент. В случае с
игрушками, тренд «моды» может привести к
устойчивому росту продаж (например, это может
быть общий тренд в сторону игрушек
образовательной направленности). Как и сезонная
компонента, этот тренд может быть по своей
природе аддитивным (продажи ежегодно
увеличиваются на 3 миллиона долларов) или
мультипликативным (продажи ежегодно
увеличиваются на 30%, или возрастают в 1.3 раза).
Кроме того, объем продаж может содержать
циклические компоненты. Повторим еще раз, что
циклическая компонента отличается от сезонной
тем, что она обычно имеет большую временную
протяженность и проявляется через неравные
промежутки времени. Так, например, некоторая
игрушка может быть особенно «горячей» в
течение летнего сезона (например, кукла,
изображающая персонаж популярного мультфильма,
которая к тому же агрессивно рекламируется). Как
и в предыдущих случаях, такая циклическая
компонента может изменять объем продаж
аддитивно, либо мультипликативно.

Вычисления

В вычислительном отношении процедура метода Сезонной
декомпозиции (Census I) следует стандартным
формулам, см. Makridakis, Wheelwright, and McGee (1983) или Makridakis and
Wheelwright (1989).

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

Отношения или разности. После взятия
скользящих средних вся сезонная (т.е. внутри
сезона) изменчивость будет исключена, и поэтому
разность (в случае аддитивной модели) или
отношение (для мультипликативной модели) между
наблюдаемым и сглаженным рядом будет выделять
сезонную составляющую (плюс нерегулярную
компоненту). Более точно, ряд скользящих средних
вычитается из наблюдаемого ряда (в аддитивной
модели) или же значения наблюдаемого ряда
делятся на значения скользящих средних (в
мультипликативной модели).

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

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

Получающийся в результате ряд называется
сезонной корректировкой ряда (из ряда убрана
сезонная составляющая)..

Тренд-циклическая компонента.
Напомним, что циклическая компонента отличается
от сезонной компоненты тем, что
продолжительность цикла, как правило, больше, чем
один сезонный период, и разные циклы могут иметь
разную продолжительность. Приближение для
объединенной тренд-циклической компоненты можно
получить, применяя к ряду с сезонной поправкой
процедуру 5-точечного (центрированного)
взвешенного скользящего среднего с весами 1, 2, 3, 2,
1.

Случайная или нерегулярная компонента.
На последнем шаге выделяется случайная или
нерегулярная компонента (погрешность) путем
вычитания из ряда с сезонной поправкой
(аддитивная модель) или делением этого ряда
(мультипликативная модель) на тренд-циклическую
компоненту.

Сезонная корректировка X-11 (метод Census II)

Общие идеи, лежащие в основе сезонной
декомпозиции и корректировки, изложены в
разделе, посвященном методу сезонной
корректировки Census I (см. Сезонная
декомпозиция (метод Census I)). Метод Census II (2)
является развитием и уточнением обычного метода
корректировки. На протяжении многих лет
различные варианты метода Census II развивались в
Бюро Переписи США (US Census Bureau); один из вариантов
этого метода, получивший широкую известность и
наиболее часто применяемый в государственных
органах и сфере бизнеса, называется «вариант X-11
метода Census II» (см. Shiskin, Young, and Musgrave, 1967).
Впоследствии этот усовершенствованный вариант
метода Census II стал называться просто X-11.
Помимо документации, которую можно получить из
Census Bureau, подробное описание метода дано в работах
Makridakis, Wheelwright and McGee (1983), Makridakis and Wheelwright (1989).

За дополнительной информацией обратитесь к
следующим разделам:

Сезонная корректировка:
основные идеи и термины
Метод Census II
Таблицы результатов
корректировки X-11
Подробное описание всех
таблиц результатов, вычисляемых в методе X-11

За дальнейшей информацией обратитесь к Анализу временных рядов и
следующим разделам:

Идентификация модели
временных рядов
АРПСС (Бокс и Дженкинс) и
автокорреляции
Прерванные временные ряды
Экспоненциальное
сглаживание
Сезонная декомпозиция (метод
Census I)
Анализ распределенных лагов
Одномерный анализ Фурье
Кросс-спектральный анализ
Основные понятия и принципы
Быстрое преобразование Фурье

Сезонная корректировка: основные идеи и
термины

Предположим, что у вас имеются ежемесячные
данные о пассажиропотоке на международных
авиалиниях за 12 лет (см. Бокс и Дженкинс, 1976). Если
изобразить эти данные на графике, то будет хорошо
видно, что (1) объем пассажиропотока имеет во
времени возрастающий линейный тренд, и что (2) в
ряде имеется ежегодно повторяющаяся
закономерность — сезонность (большинство
перевозок приходится на летние месяцы, кроме
того, имеется пик меньшей высоты в районе
декабрьских каникул). Цель сезонной декомпозиции
и корректировки как раз и состоит в том, чтобы
отделить эти компоненты, то есть разложить ряд на
составляющую тренда, сезонную компоненту и
оставшуюся нерегулярную составляющую.
«Классический» прием, позволяющий выполнить
такую декомпозицию, известен как метод Census I
(см. раздел Census I). Этот метод
описывается и обсуждается в работах Makridakis,
Wheelwright, and McGee (1983) и Makridakis and Wheelwright (1989).

Аддитивная модель:

X_t = TC_t + S_t + I_t

Мультипликативная модель:

X_t = T_t*C_t*S_t*I_t

Здесь X_t обозначает
значение временного ряда в момент времени t.

Если имеются какие-то априорные сведения о
циклических факторах, влияющих на ряд (например,
циклы деловой конъюнктуры), то можно
использовать оценки для различных компонент для
составления прогноза будущих значений ряда.
(Однако для прогнозирования предпочтительнее экспоненциальное
сглаживание, позволяющее учитывать сезонную
составляющую и тренд.)

Аддитивный и мультипликативный тренд-цикл.
Рассмотренный пример можно расширить, чтобы
проиллюстрировать понятия аддитивной и
мультипликативной тренд-циклических компонент.
В случае с игрушками, тренд «моды» может
привести к устойчивому росту продаж (например,
это может быть общий тренд в сторону игрушек
образовательной направленности). Как и сезонная
компонента, этот тренд может быть по своей
природе аддитивным (продажи ежегодно
увеличиваются на 3 миллиона долларов) или
мультипликативным (продажи ежегодно
увеличиваются на 30%, или возрастают в 1.3 раза).
Кроме того, объем продаж может содержать
циклические компоненты. Повторим еще раз, что
циклическая компонента отличается от сезонной
тем, что она обычно имеет большую временную
протяженность и проявляется через неравные
промежутки времени. Так, например, некоторая
игрушка может быть особенно «горячей» в
течение летнего сезона (например, кукла,
изображающая персонаж популярного мультфильма,
которая к тому же агрессивно рекламируется). Как
и в предыдущих случаях, такая циклическая
компонента может изменять объем продаж
аддитивно, либо мультипликативно.

Метод Census II

Основной метод сезонной декомпозиции и
корректировки, рассмотренный в разделе Сезонная корректировка:
основные идеи и термины, может быть
усовершенствован различными способами. На самом
деле, в отличие от многих методов моделирования
временных рядов (в частности, АРПСС),
которые основаны на определенной теоретической
модели, вариант X-11 метода Census II представляет
собой просто результат многочисленных специально
разработанных приемов и усовершенствований,
которые доказали свою работоспособность в
многолетней практике решения реальных задач (см.
Burman, 1979, Kendall and Ord, 1990, Makridakis and Wheelwright, 1989; Wallis, 1974).
Некоторые из наиболее важных усовершенствований
перечислены ниже.

Поправка на число рабочих дней. В
месяцах разное число дней и разное число рабочих
дней. Если мы анализируем, например, цифры
ежемесячной выручки парка аттракционов, то
разница в числе суббот и воскресений (пиковые
дни) в разных месяцах существенным образом
скажется на различиях в ежемесячных показателях
дохода. Вариант X-11 метода Census II дает
пользователю возможность проверить,
присутствует ли во временном ряду этот эффект
числа рабочих дней, и если да, то внести
соответствующие поправки.

Выбросы. Большинство реальных
временных рядов содержит выбросы, то есть резко
выделяющиеся наблюдения, вызванные какими-то
исключительными событиями. Например, забастовка
персонала может сильно повлиять на месячные или
годовые показатели выпуска продукции фирмы.
Такие выбросы могут исказить оценки сезонной
компоненты и тренда. В процедуре X-11
предусмотрены корректировки на случай появления
выбросов, основанные на использовании
«принципов статистического контроля»:
значения, выходящие за определенный диапазон
(который определяется в терминах, кратных сигма,
т.е. стандартных отклонений), могут быть
преобразованы или вовсе пропущены, и только
после этого будут вычисляться окончательные
оценки параметров сезонности.

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

Критерии и итоговые статистики. Помимо
оценки основных компонент ряда, можно вычислить
различные сводные статистики. Например, можно
сформировать таблицы дисперсионного анализа для
проверки значимости фактора сезонной
изменчивости и ряда и фактора рабочих дней (см.
выше), процедура метода X-11 вычисляет также
ежемесячные относительные изменения в случайной
и тренд-циклической компонентах. С увеличением
продолжительности временного промежутка,
измеряемого в месяцах или, в случае квартального
варианта метода X-11 — в кварталах года,
изменения в тренд-циклической компоненте, вообще
говоря, будут нарастать, в то время как изменения
случайной составляющей должны оставаться
примерно на одном уровне. Средняя длина
временного интервала, на котором изменения
тренд-циклической компоненты становятся
примерно равными изменениям случайной
компоненты, называется месяцем (кварталом)
циклического доминирования, или сокращенно
МЦД (соответственно КЦД). Например, если МЦД равно
двум, то на сроках более двух месяцев
тренд-циклическая компонента станет
доминировать над флуктуациями нерегулярной
(случайной) компоненты. Эти и другие результаты
более подробно будут обсуждаться далее.

Таблицы результатов
корректировки X-11

Вычисления, которые производятся в процедуре
X-11, лучше всего обсуждать в контексте
таблиц результатов, которые при этом выдаются.
Процедура корректировки разбивается на семь
этапов, которые обычно обозначаются буквами A
— G.

Априорная корректировка (помесячная сезонная
корректировка). Перед тем, как к временному
ряду, содержащему ежемесячные значения, будет
применяться какая-либо сезонная корректировка,
могут быть произведены различные корректировки,
заданные пользователем. Можно ввести еще один
временной ряд, содержащий априорные
корректирующие факторы; значения этого ряда
будут вычитаться из исходного ряда (аддитивная
модель), или же значения исходного ряда будут
поделены на значения корректирующего ряда
(мультипликативная модель). В случае
мультипликативной модели пользователь может
также определить свои собственные поправочные
коэффициенты (веса) на число рабочих дней. Эти
веса будут использоваться для корректировки
ежемесячных наблюдений, так чтобы учитывалось
число рабочих дней в этом месяце.
Предварительное оценивание вариации числа
рабочих дней (месячный вариант X-11) и весов. На
следующем шаге вычисляются предварительные
поправочные коэффициенты на число рабочих дней
(только в месячном варианте X-11) и веса,
позволяющие уменьшить эффект выбросов.
Окончательное оценивание вариации числа
рабочих дней и нерегулярных весов (месячный
вариант X-11). Поправки и веса, вычисленные в
пункте B, используются для построения
улучшенных оценок тренд-циклической и сезонной
компонент. Эти улучшенные оценки используются
для окончательного вычисления факторов числа
рабочих дней (в месячном варианте X-11) и весов.
Окончательное оценивание сезонных факторов,
тренд-циклической, нерегулярной и сезонно
скорректированной компонент ряда. Окончательные
значения факторов рабочих дней и весов,
вычисленные в пункте C, используются для
вычисления окончательных оценок для компонент
ряда.
Модифицированные ряды: исходный, сезонно
скорректированный и нерегулярный. Исходный и
окончательный сезонно скорректированный ряды, а
также нерегулярная компонента модифицируются
путем сглаживания выбросов. Полученные в
результате этого, модифицированные ряды
позволяют пользователю проверить устойчивость
сезонной корректировки.
Месяц (квартал) циклического доминирования
(МЦД, КЦД), скользящее среднее и сводные
показатели. IНа этом этапе вычислений
рассчитываются различные сводные
характеристики (см. далее), позволяющие
пользователю исследовать относительную
важность разных компонент, среднюю флуктуацию от
месяца к месяцу (от квартала к кварталу), среднее
число идущих подряд изменений в одну сторону и
др.
Графики. Наконец, вы можете построить
различные графики итоговых результатов.
Например, можно построить окончательно
скорректированный ряд в хронологическом порядке
или по месяцам (см. ниже).

Подробное описание всех таблиц
результатов, вычисляемых в методе X-11

На каждом из этапов A — G (см. раздел Таблицы результатов
корректировки X-11) вычислялись различные
таблицы результатов. Обычно все они нумеруются, а
также им приписывается буква, соответствующая
этапу анализа. Например, таблица B 11 содержит
предварительно сезонно скорректированный ряд; C
11 — это более точно сезонно скорректированный
ряд, а D 11 — окончательный сезонно
скорректированный ряд. Далее приводится
перечень всех таблиц. Таблицы, помеченные
звездочкой (*), недоступны (или неприменимы) при
анализе квартальных показателей. Кроме того, в
случае квартальной корректировки некоторые из
описанных ниже вычислений несколько
видоизменяются. Так, например, для вычисления
сезонных факторов вместо 12-периодного (т.е.
12-месячного) скользящего среднего используется
4-периодное (4-квартальное) скользящее среднее;
предварительная тренд-циклическая компонента
вычисляется по центрированному 4-периодному
скользящему среднему, а окончательная оценка
тренд-циклической компоненты вычисляется по
5-точечному среднему Хендерсона.

В соответствии со стандартом метода X-11,
принятым Бюро переписи США, предусмотрены три
степени подробности вывода: Стандартный (17 — 27
таблиц), Длинный (27 — 39 таблиц) и Полный (44 —
59 таблиц). Имеется также возможность выводить
только таблицы результатов, выбранные
пользователем. В следующих далее описаниях
таблиц, буквы С, Д и П рядом с
названием таблицы указывают, какие таблицы
выводятся и/или распечатываются в
соответствующем варианте вывода. (Для графиков
предусмотрены два уровня подробности вывода: Стандартный
и Все.)

Щелкните на имени таблицы для
получения информации о ней.

* A 1.
Исходный ряд (С)

*
A 2. Априорные месячные поправки (С)

*
A 3. Исходный ряд, скорректированный с помощью
априорных месячных поправок (С)

* A
4. Априорные поправки на рабочие дни (С)

B
1. Ряд после априорной корректировки либо
исходный ряд (С)

B 2.
Тренд-цикл (Д)

B
3. Немодифицированные S-I разности или отношения (П)

B
4. Значения для замены выбросов S-I разностей
(отношений) (П)

B 5.
Сезонная составляющая (П)

B 6.
Сезонная корректировка ряда (П)

B 7.
Тренд-цикл (Д)

B 8.
Немодифицированные S-I разности (отношения) (П)

B
9. Значения для замены выбросов S-I разностей
(отношений) (П)

B 10.
Сезонная составляющая (Д)

B 11.
Сезонная корректировка ряда (П)

B 12. (не используется)

B 13.
Нерегулярная составляющая ряда (Д)

Таблицы B 14 — B 16, B 18 и B 19: Поправка на число
рабочих дней. Эти таблицы доступны только при
анализе ежемесячных данных. Число разных дней
недели (понедельников, вторников и т.д.)
колеблется от месяца к месяцу. Бывают ряды, в
которых различия в числе рабочих дней в месяце
могут давать заметный разброс ежемесячных
показателей (например, месячный доход парка
аттракционов сильно зависит от того, сколько в
этом месяце было выходных дней). Пользователь
имеет возможность определить начальные веса для
каждого дня недели (см. A
4), и/или эти веса могут быть оценены по данным
(пользователь также может сделать использование
этих весов условным, т.е. только в тех случаях,
когда они объясняют значительную часть
дисперсии).

*
B 14. Выбросы нерегулярной составляющей,
исключенные из регрессии рабочих дней (Д)

* B
15. Предварительная регрессия рабочих дней (Д)

*
B 16. Поправки на число рабочих дней, полученные из
коэффициентов регрессии (П)

B
17. Предварительные веса нерегулярной компоненты (Д)

*
B 18. Поправки на число рабочих дней, полученные из
комбинированных весов дней недели (П)

*
B 19. Исходный ряд с поправками на рабочие дни и
априорную вариацию (П)

C
1. Исходный ряд, модифицированный с помощью
предварительных весов, с поправкой на рабочие
дни и априорную вариацию (Д)

C 2.
Тренд-цикл (П)

C 3. (не используется)

C
4. Модифицированные S-I разности (отношения) (П)

C 5.
Сезонная составляющая (П)

C 6.
Сезонная корректировка ряда (П)

C 7.
Тренд-цикл (Д)

C 8. (не используется)

C
9. Модифицированные S-I разности (отношения) (П)

C 10.
Сезонная составляющая (Д)

C 11.
Сезонная корректировка ряда (П)

C 12. (не используется)

C 13.
Нерегулярная составляющая (С)

Таблицы C 14 — C 16, C 18 и C 19: Поправка на число
рабочих дней. Эти таблицы доступны только при
анализе ежемесячных данных и если при этом
требуется поправка на различное число рабочих
дней. В этом случае поправки на число рабочих
дней вычисляются по уточненным значениям
сезонно скорректированных рядов аналогично
тому, как это делалось в пункте B (B 14 — B 16,
B 18, B 19).

*
C 14. Выбросы нерегулярной составляющей,
исключенные из регрессии рабочих дней (С)

* C
15. Регрессия рабочих дней — окончательный вариант
(С)

*
C 16. Поправки на число рабочих дней, полученные из
коэффициентов регрессии, — окончательный вариант
(С)

C
17. Окончательные веса нерегулярной компоненты (С)

*
C 18. Поправки на число рабочих дней, полученные из
комбинированных весов дней недели —
окончательный вариант (С)

*
C 19. Исходный ряд с поправками на рабочие дни и
априорную вариацию (С)

D
1. Исходный ряд, модифицированный с помощью
окончательных весов, с поправкой на рабочие дни и
априорную вариацию (Д)

D 2.
Тренд-цикл (П)

D 3. (не используется)

D
4. Модифицированные S-I разности (отношения) (П)

D 5.
Сезонная составляющая (П)

D 6.
Сезонная корректировка ряда (П)

D 7.
Тренд-цикл (Д)

D
8. Немодифицированные S-I разности (отношения) —
окончательный вариант (С)

D
9. Окончательные значения для замены выбросов S-I
разностей (отношений) (С)

D 10.
Сезонная составляющая — окончательный вариант (С)

D
11. Сезонная корректировка ряда — окончательный
вариант (С)

D 12.
Тренд-циклическая компонента — окончательный
вариант (С)

D 13.
Нерегулярная составляющая — окончательный
вариант (С)

E 1.
Модифицированный исходный ряд (С)

E 2.
Модифицированный ряд с сезонной поправкой (С)

E 3.
Модифицированная нерегулярная составляющая (С)

E
4. Разности (отношения) годовых сумм (С)

E
5. Разности (относительные изменения) исходного
ряда (С)

E
6. Разности (относительные изменения)
окончательного варианта ряда с сезонной
поправкой (С)

F 1.
МЦД (КЦД) скользящее среднее (С)

F 2.
Сводные показатели (С)

G 1. График (С)

G 2. График (С)

G 3. График (В)

G 4. График (В)

Анализ распределенных лагов

Общая цель
Общая модель
Распределенный лаг Алмона

За дальнейшей информацией обратитесь к Анализу временных рядов и
следующим разделам:

Идентификация модели
временных рядов
АРПСС (Бокс и Дженкинс) и
автокорреляции Вводный обзор АРПСС
Прерванные временные ряды
Экспоненциальное
сглаживание
Сезонная декомпозиция (метод
Census I)
Сезонная корректировка X-11
(метод Census II)
Таблицы результатов
корректировки X-11
Одномерный анализ Фурье
Кросс-спектральный анализ
Основные понятия и принципы
Быстрое преобразование Фурье

Общая цель

Анализ распределенных лагов — это специальный
метод оценки запаздывающей зависимости между
рядами. Например, предположим, вы производите
компьютерные программы и хотите установить
зависимость между числом запросов, поступивших
от покупателей, и числом реальных заказов. Вы
могли бы записывать эти данные ежемесячно в
течение года и затем рассмотреть зависимость
между двумя переменными: число запросов и число
заказов зависит от запросов, но зависит с
запаздыванием. Однако очевидно, что запросы
предшествуют заказам, поэтому можно ожидать, что
число заказов. Иными словами, в зависимости между
числом запросов и числом продаж имеется

временной сдвиг (лаг) (см. также автокорреляции и
кросскорреляции).

Такого рода зависимости с запаздыванием
особенно часто возникают в эконометрике.
Например, доход от инвестиций в новое
оборудование отчетливо проявится не сразу, а
только через определенное время. Более высокий
доход изменяет выбор жилья людьми; однако эта
зависимость, очевидно, тоже проявляется с
запаздыванием. [Подобные задачи возникают в
страховании, где временной ряд клиентов и ряд
денежных поступлений сдвинуты друг относительно
друга].

Во всех этих случаях, имеется независимая или объясняющая
переменная, которая воздействует на зависимые
переменные с некоторым запаздыванием (лагом).
Метод распределенных лагов позволяет
исследовать такого рода зависимость.

Подробные обсуждения зависимостей с
распределенными лагами имеются в
эконометрических учебниках, например, в Judge, Griffith,
Hill, Luetkepohl, and Lee (1985), Maddala (1977), and Fomby, Hill, and Johnson (1984).
Ниже дается краткое описание этих методов.
Предполагается, что вы знакомы с понятием
корреляции (см. Основные
статистики и таблицы), кросскорреляции и
основными идеями множественной регрессии (см. Множественная регрессия).

Общая модель

Пусть y — зависимая переменная, a независимая
или объясняющая x. Эти переменные измеряются
несколько раз в течение определенного отрезка
времени. В некоторых учебниках по эконометрике
зависимая переменная называется также эндогенной
переменной, a зависимая или объясняемая
переменная экзогенной переменной.
Простейший способ описать зависимость между
этими двумя переменными дает следующее линейное
уравнение:

Y_t = _i*x_t-i

В этом уравнении значение зависимой переменной
в момент времени t является
линейной функцией переменной x,
измеренной в моменты t, t-1,
t-2 и т.д. Таким образом, зависимая
переменная представляет собой линейные функции x и x, сдвинутых на
1, 2, и т.д. временные периоды. Бета коэффициенты
(i) могут
рассматриваться как параметры наклона в этом
уравнении. Будем рассматривать это уравнение как
специальный случай уравнения линейной регрессии
(см. раздел Множественная
регрессия). Если коэффициент переменной с
определенным запаздыванием (лагом) значим, то
можно заключить, что переменная y
предсказывается (или объясняется) с
запаздыванием.

Распределенный лаг Алмона

Обычная проблема, возникающая в множественной
регрессии, состоит в том, что соседние значения x
сильно коррелируют. В самом крайнем случае, это
приводит к тому, что корреляционная матрица не
будет обратимой и коэффициенты бета не могут
быть вычислены. В менее экстремальных ситуациях
вычисления этих коэффициентов и их стандартные
ошибки становятся ненадежными из-за
вычислительных ошибок (ошибок округления). В
контексте множественной регрессии эта проблема
хорошо известна как проблема мультиколлинеарности
(см. раздел Множественная
регрессия).

Алмон (1965) предложил специальную процедуру,
которая в данном случае уменьшает
мультиколлинеарность. Именно, пусть каждый
неизвестный коэффициент записан в виде:

_i =
₀ + ₁*i + … + _q*i^q

Алмон показал, что во многих случаях (в
частности, чтобы избежать мультиколлинеарности)
легче оценить коэффициенты альфа, чем
непосредственно коэффициенты бета. Такой
метод оценивания коэффициентов бета
называется полиномиальной аппроксимацией.

Неправильная спецификация. Общая
проблема полиномиальной аппроксимации, состоит
в том, что длина лага и степень полинома
неизвестны заранее. Последствия
неправильного определения (спецификации) этих
параметров потенциально серьезны (в силу
смещения, возникающего в оценках при
неправильном задании параметров). Этот вопрос
подробно обсуждается в книгах Frost (1975), Schmidt and Waud
(1973), Schmidt and Sickles (1975) и Trivedi and Pagan (1979).

Одномерный анализ Фурье

В спектральном анализе исследуются
периодические модели данных. Цель анализа —
разложить комплексные временные ряды с
циклическими компонентами на несколько основных
синусоидальных функций с определенной длиной
волн. Термин «спектральный» — своеобразная
метафора для описания природы этого анализа.
Предположим, вы изучаете луч белого солнечного
света, который, на первый взгляд, кажется
хаотически составленным из света с различными
длинами волн. Однако, пропуская его через призму,
вы можете отделить волны разной длины или
периодов, которые составляют белый свет.
Фактически, применяя этот метод, вы можете теперь
распознавать и различать разные источники света.
Таким образом, распознавая существенные
основные периодические компоненты, вы узнали
что-то об интересующем вас явлении. В сущности,
применение спектрального анализа к временным
рядам подобно пропусканию света через призму. В
результате успешного анализа можно обнаружить
всего несколько повторяющихся циклов различной
длины в интересующих вас временных рядах,
которые, на первый взгляд, выглядят как случайный
шум.

Наиболее известный пример применения
спектрального анализа — циклическая природа
солнечных пятен (например, см. Блумфилд, 1976 или
Шамвэй, 1988). Оказывается, что активность
солнечных пятен имеет 11-ти летний цикл. Другие
примеры небесных явлений, изменения погоды,
колебания в товарных ценах, экономическая
активность и т.д. также часто используются в
литературе для демонстрации этого метода. В
отличие от АРПСС или метода экспоненциального
сглаживания (см. разделы АРПСС
и Экспоненциальное
сглаживание), цель спектрального анализа —
распознать сезонные колебания различной длины, в
то время как в предшествующих типах анализа,
длина сезонных компонент обычно известна (или
предполагается) заранее и затем включается в
некоторые теоретические модели скользящего
среднего или автокорреляции.

Классический текст по спектральному анализу —
Bloomfield (1976); однако другие подробные обсуждения
могут быть найдены в Jenkins and Watts (1968), Brillinger (1975), Brigham
(1974), Elliott and Rao (1982), Priestley (1981), Shumway (1988) или Wei (1989).

За дальнейшей информацией обратитесь к Анализу временных рядов и
следующим разделам:

Основные понятия и принципы
Быстрое преобразование Фурье
Идентификация модели
временных рядов
АРПСС (Бокс и Дженкинс) и
автокорреляции Вводный обзор АРПСС
Прерванные временные ряды
Анализ распределенных лагов
Сезонная декомпозиция (метод
Census I)
Экспоненциальное
сглаживание
Кросс-спектральный анализ

Кросс-спектральный анализ

Общее введение
Основные понятия и принципы
Результаты для каждой
переменной
Кросс-периодограмма,
кросс-плотность, квадратурная плотность и
кросс-амплитуда
Квадрат когерентности,
усиление и фазовый сдвиг
Как создавались данные для
примера

За дальнейшей информацией обратитесь к Анализу временных рядов и
следующим разделам:

Идентификация модели
временных рядов
АРПСС (Бокс и Дженкинс) и
автокорреляции Вводный обзор АРПСС
Прерванные временные ряды
Экспоненциальное
сглаживание
Сезонная декомпозиция (метод
Census I)
Сезонная корректировка X-11
(метод Census II)
Таблицы результатов
корректировки X-11
Анализ распределенных лагов
Одномерный анализ Фурье
Основные понятия и принципы
Быстрое преобразование Фурье

Общее введение

Кросс-спектральный анализ развивает Одномерный анализ Фурье и
позволяет анализировать одновременно два ряда.
Мы предполагаем, что вы уже прочитали введение к
разделу одномерного спектрального анализа.
Подробное обсуждение кросс-спектрального
анализа можно найти в книгах Bloomfield (1976), Jenkins and Watts
(1968), Brillinger (1975), Brigham (1974), Elliott and Rao (1982), Priestley (1981),
Shumway (1988), or Wei (1989).

Периодичность ряда на определенных
частотах. Наиболее известный пример
применения спектрального анализа — циклическая
природа солнечных пятен (например, см. Блумфилд,
1976 или Шамвэй, 1988). Оказывается, что активность
солнечных пятен имеет 11-ти летний цикл. Другие
примеры небесных явлений, изменения погоды,
колебания в товарных ценах, экономическая
активность и т.д. также часто используются в
литературе для демонстрации этого метода.

Основные понятия и принципы

Простой пример. Рассмотрим следующие два
ряда с 16 наблюдениями:

	ПЕРЕМ1	ПЕРЕМ2
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16	1.000 1.637 1.148 -.058 -.713 -.383 .006 -.483 -1.441 -1.637 -.707 .331 .441 -.058 -.006 .924	-.058 -.713 -.383 .006 -.483 -1.441 -1.637 -.707 .331 .441 -.058 -.006 .924 1.713 1.365 .266

С первого взгляда нелегко рассмотреть
взаимосвязь между двумя рядами. Тем не менее, как
показано ниже, ряды создавались так, что содержат
две сильно коррелируемые периодичности. Далее
показаны части таблицы результатов из
кросс-спектрального анализа (спектральные
оценки были сглажены окном Парзена ширины 3).

Незавмсимая (X): ПЕРЕМ1 Зависимая (Y): ПЕРЕМ2
Частота	Период	X плотность	Y плотность	Кросс плотность	Кросс квадр.	Кросс амплит.
0.000000 .062500 .125000 .187500 .250000 .312500 .375000 .437500 .500000	16.00000 8.00000 5.33333 4.00000 3.20000 2.66667 2.28571 2.00000	.000000 8.094709 .058771 3.617294 .333005 .091897 .052575 .040248 .037115	.024292 7.798284 .100936 3.845154 .278685 .067630 .036056 .026633 0.000000	-.00000 2.35583 -.04755 -2.92645 -.26941 -.07435 -.04253 -.03256 0.00000	0.00000 -7.58781 .06059 2.31191 .14221 .02622 .00930 .00342 0.00000	.000000 7.945114 .077020 3.729484 .304637 .078835 .043539 .032740 0.000000

Результаты для каждой
переменной

Полная таблица результатов содержит все
спектральные статистики, вычисленные для
каждого ряда, как описано в разделе Одномерный анализ Фурье.
Взглянув на приведенные выше результаты,
очевидно, что оба ряда имеют основные
периодичности на частотах .0625 и .1875.

Кросс-периодограмма,
кросс-плотность, квадратурная плотность и
кросс-амплитуда

Аналогично результатам для одной переменной,
полная итоговая таблица результатов также
покажет значения периодограммы для
кросс-периодограммы. Однако кросс-спектр состоит
из комплексных чисел,
которые могут быть разделены на действительную и
мнимую части. Они могут быть сглажены для
вычисления оценок кросс-плотности и
квадратурной плотности (квадр-плотность для
краткости), соответственно. (Причины сглаживания
и различные функции весов для сглаживания
обсуждаются в разделе Одномерный
анализ Фурье.) Квадратный корень из суммы
квадратов значений кросс-плотности и
квадр-плотности называется кросс-амплитудой.
Кросс-амплитуда может интерпретироваться как
мера ковариации между соответствующими
частотными компонентами двух рядов. Таким
образом из результатов, показанных в таблице
результатов выше, можно заключить, что частотные
компоненты .0625 и .1875 двух рядов взаимосвязаны.

Квадрат когерентности, усиление
и фазовый сдвиг

Существуют дополнительные статистики, которые
будут показаны в полной итоговой таблице
результатов.

Квадрат когерентности. Можно
нормировать значения кросс-амплитуды, возведя их
в квадрат и разделив на произведение оценок
спектральной плотности каждого ряда. Результат
называется квадратом когерентности, который
может быть проинтерпретирован как квадрат
коэффициента корреляции (см. раздел Корреляции); т.е.
значение когерентности — это квадрат корреляции
между циклическими компонентами двух рядов
соответствующей частоты. Однако значения
когерентности не следует объяснять таким
образом; например, когда оценки спектральной
плотности обоих рядов очень малы, могут
получиться большие значения когерентности
(делитель в выражении когерентности может быть
очень маленьким), даже если нет существенных
циклических компонент в каждом ряду
соответствующей частоты.

Усиление. Значение усиления в анализе
вычисляется делением значения кросс-амплитуды
на оценки спектральной плотности одного или двух
рядов. Следовательно, может быть вычислено два
значения усиления, которые могут
интерпретироваться как стандартные
коэффициенты регрессии, соответствующей
частоты, полученные методом наименьших
квадратов.

Фазовый сдвиг. В заключение, оценки
фазового сдвига вычисляются как арктангенс (tan**-1)
коэффициента пропорциональности оценки
квадр-плотности и оценки кросс-плотности. Оценки
фазового сдвига (обычно обозначаемые греческой
буквой y) измеряют, насколько каждая частотная
компонента одного ряда опережает частотные
компоненты другого.

Как создавались данные для
примера

Теперь вернемся к примеру данных, приведенному
выше. Большие оценки спектральной плотности для
обоих рядов и значения кросс-амплитуды для
частот = 0.0625 и = .1875 предполагают две
существенных синхронных периодичности с этими
частотами в обоих рядах. Фактически, два ряда
создавались как:

v1 = cos(2**.0625*(v0-1))
+ .75*sin(2**.2*(v0-1))

v2 = cos(2**.0625*(v0+2)) +
.75*sin(2**.2*(v0+2))

(где v0 — номер наблюдения).
Действительно, анализ, представленный в этом
обзоре, очень хорошо воспроизводит
периодичность, заложенную в данные.

Спектральный анализ — Основные понятия и
принципы

Частота и период
Общая структура модели
Простой пример
Периодограмма
Проблема рассеяния
Добавление констант во
временной ряд (пэддинг)
Косинус-сглаживание
Окна данных и оценки
спектральной плотности
Подготовка данных к анализу
Результаты для случая, когда в
ряде отсутствует периодичность

За дальнейшей информацией обратитесь к Анализу временных рядов и
следующим разделам:

Идентификация модели
временных рядов
АРПСС (Бокс и Дженкинс) и
автокорреляции Вводный обзор АРПСС
Прерванные временные ряды
Экспоненциальное
сглаживание
Сезонная декомпозиция (метод
Census I)
Сезонная корректировка X-11
(метод Census II)
Таблицы результатов
корректировки X-11
Анализ распределенных лагов
Одномерный анализ Фурье
Кросс-спектральный анализ
Быстрое преобразование Фурье

Частота и период

Длина волны функций синуса или косинуса, как
правило, выражается числом циклов (периодов) в
единицу времени (Частота), часто обозначается
греческой буквой ню (; в некоторых учебниках также
используют f). Например, временной ряд,
состоящий из количества писем, обрабатываемых
почтой, может иметь 12 циклов в году. Первого числа
каждого месяца отправляется большое количество
корреспонденции (много счетов приходит именно
первого числа каждого месяца); затем, к середине
месяца, количество корреспонденции уменьшается;
и затем вновь возрастает к концу месяца. Поэтому
каждый месяц колебания в количестве
корреспонденции, обрабатываемой почтовым
отделением, будут проходить полный цикл. Таким
образом, если единица анализа — один год, то будет равно 12 (поскольку
имеется 12 циклов в году). Конечно, могут быть и
другие циклы с различными частотами. Например,
годичные циклы (=1)
и, возможно, недельные циклы (=52 недели в год).

Период Т функций синуса или косинуса
определяется как продолжительность по времени
полного цикла. Таким образом, это обратная
величина к частоте: T = 1/. Возвратимся к примеру с почтой из
предыдущего абзаца, здесь месячный цикл будет
равен 1/12 = 0.0833 года. Другими словами, это период
составляет 0.0833 года.

Общая структура модели

Как было отмечено ранее, цель спектрального
анализа — разложить ряд на функции синусов и
косинусов различных частот, для определения тех,
появление которых особенно существенно и
значимо. Один из возможных способов сделать это —
решить задачу линейной множественной
регрессии (см. раздел Множественная
регрессия), где зависимая переменная
-наблюдаемый временной ряд, а независимые
переменные или регрессоры: функции синусов всех
возможных (дискретных) частот. Такая модель
линейной множественной регрессии может быть
записана как:

x_t = a₀ + [a_k*cos(_k*t) + b_k*sin(_k*t)] (для k = 1 до q)

Следующее общее понятие классического
гармонического анализа в этом уравнении — (лямбда) -это
круговая частота, выраженная в радианах в
единицу времени, т.е. = 2**_k, где — константа пи =
3.1416 и _k = k/q. Здесь
важно осознать, что вычислительная задача
подгонки функций синусов и косинусов разных длин
к данным может быть решена с помощью
множественной линейной регрессии. Заметим, что
коэффициенты a_k при
косинусах и коэффициенты b_k
при синусах — это коэффициенты регрессии,
показывающие степень, с которой соответствующие
функции коррелируют с данными [заметим, что сами
синусы и косинусы на различных частотах не
коррелированы или, другим языком, ортогональны.
Таким образом, мы имеем дело с частным случаем
разложения по ортогональным полиномам.] Всего
существует q различных синусов
и косинусов (см. также Множественная регрессия);
интуитивно ясно, что число функций синусов и
косинусов не может быть больше числа данных в
ряде. Не вдаваясь в подробности, отметим, если n —
количество данных, то будет n/2+1 функций
косинусов и n/2-1 функций синусов. Другими
словами, различных синусоидальных волн будет
столько же, сколько данных, и вы сможете
полностью воспроизвести ряд по основным
функциям. (Заметим, если количество данных в ряде
нечетно, то последнее наблюдение обычно
опускается. Для определения синусоидальной
функции нужно иметь, по крайней мере, две точки:
высокого и низкого пика.)

В итоге, спектральный анализ определяет
корреляцию функций синусов и косинусов
различной частоты с наблюдаемыми данными. Если
найденная корреляция (коэффициент при
определенном синусе или косинусе) велика, то
можно заключить, что существует строгая
периодичность на соответствующей частоте в
данных.

Комплексные числа (действительные и мнимые
числа). Во многих учебниках по спектральному
анализу структурная модель, показанная выше,
представлена в комплексных числах; т.е. параметры
оцениваемого процесса описаны с помощью
действительной и мнимой части преобразования
Фурье. Комплексное число состоит из
действительного и мнимого числа. Мнимые числа, по
определению, — это числа, умноженные на константу i,
где i определяется как квадратный корень из -1.
Очевидно, корень квадратный из -1 не существует в
обычном сознании (отсюда термин мнимое число);
однако арифметические операции над мнимыми
числами могут производиться естественным
образом [например, (i*2)**2= -4]. Полезно представление
действительных и мнимых чисел, образующих
двумерную координатную плоскость, где
горизонтальная или X-ось представляет все
действительные числа, а вертикальная или Y-ось
представляет все мнимые числа. Комплексные числа
могут быть представлены точками на двумерной
плоскости. Например, комплексное число 3+i*2 может
быть представлено точкой с координатами {3,2} на
этой плоскости. Можно также представить
комплексные числа как углы; например, можно
соединить точку, соответствующую комплексному
числу на плоскости с началом координат
(комплексное число 0+i*0), и измерить угол наклона
этого вектора к горизонтальной оси. Таким
образом интуитивно ясно, каким образом формула
спектрального разложения, показанная выше, может
быть переписана в комплексной области. В таком
виде математические вычисления часто более
изящны и проще в выполнении, поэтому многие
учебники предпочитают представление
спектрального анализа в комплексных числах.

Простой пример

Шамвэй (1988) предлагает следующий простой пример
для объяснения спектрального анализа. Создадим
ряд из 16 наблюдений, полученных из уравнения,
показанного ниже, а затем посмотрим, каким
образом можно извлечь из него информацию.
Сначала создадим переменную и определим ее как:

x = 1*cos(2**.0625*(v0-1))
+ .75*sin(2**.2*(v0-1))

Эта переменная состоит из двух основных
периодичностей — первая с частотой =.0625 (или периодом 1/=16; одно наблюдение
составляет 1/16-ю длины полного цикла, или весь
цикл содержит каждые 16 наблюдений) и вторая с
частотой =.2 (или
периодом 5). Коэффициент при косинусе (1.0) больше
чем коэффициент при синусе (.75). Итоговая таблица
результатов спектрального анализа показана
ниже.

	Спектральный анализ: ПЕРЕМ1 (shumex.sta) Число наблюдений: 16
t	Час- тота	Период	Косинус корэфф.	Синус корэфф.	Периодо- грамма
0 1 2 3 4 5 6 7 8	.0000 .0625 .1250 .1875 .2500 .3125 .3750 .4375 .5000	16.00 8.00 5.33 4.00 3.20 2.67 2.29 2.00	.000 1.006 .033 .374 -.144 -.089 -.075 -.070 -.068	0.000 .028 .079 .559 -.144 -.060 -.031 -.014 0.000	.000 8.095 .059 3.617 .333 .092 .053 .040 .037

Теперь рассмотрим столбцы таблицы результатов.
Ясно, что наибольший коэффициент при косинусах
расположен напротив частоты .0625. Наибольший
коэффициент при синусах соответствует частоте
.1875. Таким образом, эти две частоты, которые были
«внесены» в данные, отчетливо проявились.

Периодограмма

Функции синусов и косинусов независимы (или
ортогональны); поэтому можно просуммировать
квадраты коэффициентов для каждой частоты, чтобы
вычислить периодограмму. Более часто,
значения периодограммы вычисляются как:

P_k = синус-коэффициент_k²
+ косинус-коэффициент_k² * N/2

где P_k — значения
периодограммы на частоте _k , и N — общая
длина ряда. Значения периодограммы можно
интерпретировать как дисперсию (вариацию) данных
на соответствующей частоте. Обычно значения
периодограммы изображаются в зависимости от
частот или периодов.

Проблема рассеяния

В примере, приведенном выше, функция синуса с
частотой 0.2 была «вставлена» в ряд. Однако
из-за того, что длина ряда равна 16, ни одна из
частот, полученных в таблице результатов, не
совпадает в точности с этой частотой. На практике
в этих случаях часто оказывается, что
соответствующая частота «рассеивается» на
близкие частоты. Например, могут быть найдены
большие значения периодограммы для двух близких
частот, когда в действительности существует
только одна основная функция синуса или косинуса
с частотой, которая попадает на одну из этих
частот или лежит между найденными частотами.
Существует три подхода к решению проблемы
рассеяния:

При помощи добавление констант во временной
ряда ряда можно увеличить частоты,
Применяя сглаживание ряда перед анализом,
можно уменьшить рассеяние или
Применяя сглаживание периодограммы, можно
идентифицировать основные частотные области или
(спектральные плотности), которые
существенно влияют на циклическое поведение
ряда.

Ниже смотрите описание каждого из этих
подходов.

Добавление констант во
временной ряд (пэддинг)

Так как частотные величины вычисляются как N/t,
можно просто добавить в ряд константы (например,
нули), и таким образом получить увеличение
частот. Фактически, если вы добавите в файл
данных, описанный в примере выше, десять нулей,
результаты не изменятся; т.е. наибольшие пики
периодограммы будут находиться по-прежнему на
частотах близких к .0625 и .2. (Добавление констант
во временной ряд также часто желательно для
увеличения вычислительной эффективности; см.
ниже.)

Косинус-сглаживание

Так называемый процесс косинус-сглаживания —
рекомендуемое преобразование ряда,
предшествующее спектральному анализу. Оно
обычно приводит к уменьшению рассеяния в
периодограмме. Логическое обоснование этого
преобразования подробно объясняется в книге
Bloomfield (1976, стр. 80-94). По существу, количественное
отношение (p) данных в начале и в конце ряда
преобразуется при помощи умножения на веса:

w_t = 0.5*{1-cos[*(t — 0.5)/m]} (для t=0 до m-1)
w_t = 0.5*{1-cos[*(N — t +
0.5)/m]} (для t=N-m до N-1)

где m выбирается так, чтобы 2*m/N было равно коэффициенту
пропорциональности сглаживаемых данных (p).

Окна данных и оценки
спектральной плотности

На практике, при анализе данных обычно не очень
важно точно определить частоты основных функций
синусов или косинусов. Скорее, т.к. значения
периодограммы — объект существенного случайного
колебания, можно столкнуться с проблемой многих
хаотических пиков периодограммы. В этом случае
хотелось бы найти частоты с большими спектральными
плотностями, т.е. частотные области, состоящие
из многих близких частот, которые вносят
наибольший вклад в периодическое поведение
всего ряда. Это может быть достигнуто путем
сглаживания значений периодограммы с помощью
преобразования взвешенного скользящего
среднего. Предположим, ширина окна скользящего
среднего равна m (должно быть нечетным
числом); следующие наиболее часто используемые
преобразования (заметим: p = (m-1)/2).

Окно Даниэля (равные веса). Окно
Даниэля (Daniell, 1946) означает простое (с равными
весами) сглаживание скользящим средним значений
периодограммы; т.е. каждая оценка спектральной
плотности вычисляется как среднее m/2
предыдущих и последующих значений
периодограммы.

Окно Тьюки. В окне Тьюки (Blackman and Tukey, 1958)
или Тьюки-Ханна (Hanning) (названное в честь Julius Von Hann),
для каждой частоты веса для взвешенного
скользящего среднего значений периодограммы
вычисляются как:

w_j = 0.5 + 0.5*cos(*j/p) (для j=0 до p)
w_-j = w_j (для j 0)

Окно Хемминга. В окне Хемминга
(названного в честь R. W. Hamming) или Тьюки-Хемминга
(Blackman and Tukey, 1958), для каждой частоты, веса для
взвешенного скользящего среднего значений
периодограммы вычисляются как:

w_j = 0.54 + 0.46*cos(*j/p) (для j=0 до p)
w_-j = w_j (для j 0)

Окно Парзена. В окне Парзена (Parzen, 1961),
для каждой частоты, веса для взвешенного
скользящего среднего значений периодограммы
вычисляются как:

w_j = 1-6*(j/p)² + 6*(j/p)³    (для
j = 0 до p/2)
w_j = 2*(1-j/p)³    (для j = p/2 + 1 до p)
w_-j = w_j    (для j 0)

Окно Бартлетта. В окне Бартлетта (Bartlett,
1950) веса вычисляются как:

w_j = 1-(j/p) (для j = 0 до p)
w_-j = w_j (для j 0)

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

Подготовка данных к анализу

Теперь рассмотрим несколько других
практических моментов спектрального анализа.
Обычно, полезно вычесть среднее из значений ряда
и удалить тренд (чтобы добиться стационарности)
перед анализом. Иначе периодограмма и
спектральная плотность «забьются» очень
большим значением первого коэффициента при
косинусе (с частотой 0.0). По существу, среднее — это
цикл частоты 0 (нуль) в единицу времени; т.е.
константа. Аналогично, тренд также не
представляет интереса, когда нужно выделить
периодичность в ряде. Фактически оба этих
эффекта могут заслонить более интересные
периодичности в данных, поэтому и среднее, и
(линейный) тренд следует удалить из ряда перед
анализом. Иногда также полезно сгладить данные
перед анализом, чтобы убрать случайный шум,
который может засорять существенные
периодические циклы в периодограмме.

Результаты для случая, когда в
ряде отсутствует периодичность

В заключение, зададим вопрос: что, если
повторяющихся циклов в данных нет, т.е. если
каждое наблюдение совершенно независимо от всех
других наблюдений? Если распределение
наблюдений соответствует нормальному, такой
временной ряд может быть белым шумом (подобный
белый шум можно услышать, настраивая радио). Если
исходный ряд — белый шум, то значения
периодограммы будут иметь экспоненциальное
распределение. Таким образом, проверкой на
экспоненциальность значений периодограммы
можно узнать, отличается ли исходный ряд от
белого шума. Пользователь может также построить
одновыборочную статистику d статистику
Колмогорова-Смирнова (cм. также раздел Непараметрическая статистика и
распределения).

Проверка, что шум — белый в ограниченной
полосе частот. Заметим, что также можно
получить значения периодограммы для
ограниченной частотной области. Снова, если
введенный ряд — белый шум с соответствующими
частотами (т.е. если нет существенных
периодических циклов этих частот), то
распределение значений периодограммы должно
быть снова экспоненциальным.

Быстрое преобразование Фурье (БПФ)

Общее введение
Вычисление БПФ во временных
рядах

За дальнейшей информацией обратитесь к Анализу временных рядов и
следующим разделам:

Идентификация модели
временных рядов
АРПСС (Бокс и Дженкинс) и
автокорреляции Вводный обзор АРПСС
Прерванные временные ряды
Экспоненциальное
сглаживание
Сезонная декомпозиция (метод
Census I)
Сезонная корректировка X-11
(метод Census II)
Таблицы результатов
корректировки X-11
Анализ распределенных лагов
Одномерный анализ Фурье
Кросс-спектральный анализ
Основные понятия и принципы

Общее введение

Интерпретация результатов спектрального
анализа обсуждается в разделе Основные
понятия и принципы, однако там мы не
обсуждали вычислительные проблемы, которые в
действительности очень важны. До середины 1960-х
для представления спектрального разложения
использовались точные формулы для нахождения
параметров синусов и косинусов. Соответствующие
вычисления требовали как минимум N**2 (комплексных)
умножений. Таким образом, даже сегодня
высокоскоростному компьютеру потребовалось бы
очень много времени для анализа даже небольшого
временного ряда (для 8,000 наблюдений
потребовалось бы по меньшей мере 64 миллиона
умножений).

Ситуация кардинально изменилась с открытием
так называемого алгоритма
быстрого преобразования Фурье, или БПФ для
краткости. Достаточно сказать, что при
применении алгоритма БПФ время выполнения
спектрального анализа ряда длины N стало
пропорционально N*log2(N) что конечно
является огромным прогрессом.

Однако недостаток стандартного алгоритма БПФ
состоит в том, что число данных ряда должно быть
равным степени 2 (т.е. 16, 64, 128, 256, …). Обычно это
приводит к необходимости добавлять нули во
временной ряд, который, как описано выше, в
большинстве случаев не меняет характерные пики
периодограммы или оценки спектральной
плотности. Тем не менее, в некоторых случаях,
когда единица времени значительна, добавление
констант во временной ряд может сделать
результаты более громоздкими.

Вычисление БПФ во временных
рядах

Выполнение быстрого преобразования Фурье
чрезвычайно эффективно. На большинстве
стандартных компьютеров, ряд с более чем 100,000
наблюдений легко анализируется. Однако
существует несколько моментов, которые надо
помнить при анализе рядов большого размера.

Как упоминалось ранее, для применения
стандартного (и наиболее эффективного) алгоритма
БПФ требуется, чтобы длина исходного ряда была
равна степени 2. Если это не так, должны быть
проведены дополнительные вычисления. Будут
использоваться простые точные вычислительные
формулы, пока исходный ряд относительно мал, и
вычисления можно выполнить за относительно
короткое время. Для длинных временных рядов,
чтобы применить алгоритм БПФ, используется
основной подход, описанный Monro и Branch (1976). Этот
метод требует значительно больше памяти; однако
ряд рассматриваемой длины может анализироваться
все еще очень быстро, даже если число наблюдений
не является степенью 2.

Для временных рядов, длина которых не равна
степени 2, мы можем дать следующие рекомендации:
если размер исходного ряда не превосходит
средний размер (т.е. имеется только несколько
тысяч наблюдений), не стоит беспокоиться. Анализ
займет несколько секунд. Для анализа средних и
больших рядов (например, содержащих свыше 100,000
наблюдений), добавьте в ряд константы (например
нули) до тех пор, пока длина ряда не станет
степенью 2 и затем примените косинус-сглаживание
ряда в разведочной части анализа ваших данных.

Дополнительная информация по методам анализа данных, добычи данных,
визуализации и прогнозированию содержится на
Портале StatSoft (http://www.statsoft.ru/home/portal/default.asp)
и в Углубленном Учебнике StatSoft (Учебник с формулами).

Все права на материалы электронного учебника принадлежат компании StatSoft

Источник

Рост точности прогноза — это точка роста оборотных средств при неизменном объеме. Чем меньше ошибка прогноза, тем меньше денег необходимо на обслуживание модели прогноза.

Дополнительные оборотные средства за счет повышения точности прогноза мы получим, если будем использовать модель прогнозирования, которая дает наименьшую среднюю абсолютную ошибку прогноза.

В данной статье мы рассмотрим:

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

По ходу статьи мы разберем

Что такое ошибка прогноза;
Как рассчитывается среднее абсолютное отклонение;

и рассчитаем:

Прогноз с помощью модели «Скользящей средней к 4-м месяцам с аддитивной сезонностью»;
Прогноз с помощью модели «Логарифмический тренд с сезонностью»;
Ошибку прогноза для каждой модели;
Среднее абсолютное отклонение и для каждой модели.

А также сравним модели, опираясь на среднее абсолютное отклонение и оценим экономию оборотных средств за счет использования более точной модели прогнозирования.

Скачайте файл с примером

Что такое ошибка прогноза?

Ошибкой прогноза продаж является разность между фактическими продажами и прогнозом продаж.

Чем меньше ошибка прогноза, тем более точные решения мы приминаем в закупках, производстве, планировании … а следовательно более эффективно распределяем оборотные средства и повышаем оборачиваемость товаров.

Существует несколько методов оценки ошибок. Большинство этих методов состоит в усреднении некоторых функций от разностей между действительными значениями и их прогнозами.

Ошибку прогноза (et) для каждого момента времени во временном ряду мы можем вычистить по формуле:

et = Yt — Y^t ,

где

Yt — действительное значение временного ряда в момент t — в наших примерах объем продаж,
Y^t — прогноз значения Yt — в наших примерах прогноз объема продаж.

Ошибка MAD — среднее абсолютное отклонение

Среднее абсолютное отклонение (MAD) измеряет точность прогноза, усредняя величины ошибок прогноза (абсолютные значения каждой ошибки). Чаще всего MAD используют, когда ошибку прогноза необходимо измерить в тех же единицах, что и исходные значения временного ряда.

Формула вычисления ошибки:

Yt — действительное значение временного ряда в момент t,
Y^t — прогноз значения Yt,
n — номера периодов

Среднее абсолютное отклонение — средняя ошибка (разность между фактом продаж и прогнозом продаж) по модулю.

Рассчитаем прогноз и оценим следующие модели:

Скользящей средней к 4-м месяцам с аддитивной сезонностью;
Логарифмический тренд с сезонностью.

Скачайте файл с примером

Для оценки ошибки модели «Скользящей средней к 4-м месяцам с аддитивной сезонностью» рассчитаем:

Скользящую среднюю к 4-м месяцам;
Разность между значениями ряда и средними значениями к 4-м месяцам (пункт 1);
Усредним разность ряда и средней для каждого месяца получим сезонность в абсолютной величание — аддитивную сезонность;
Продлим значения ряда с помощью скользящей средней к 4-м месяцам и скорректируем её аддитивной сезонностью;
Модель прогноза для каждого момента времени t;
Ошибку прогноза;
Среднее абсолютное отклонение.

1. Скользящую среднюю к 4-м месяцам для каждого момента времени во временном ряду начиная с 5-го периода:

2. Разность между значениями ряда и средними значениями к 4-м месяцам для каждого момента времени t (пункт 1):

3. Усредним разность ряда и средней для каждого месяца получим сезонность в абсолютной величание — аддитивную сезонность.

Для этого вначале выделим номера месяцев с помощью функции Excel =месяц(дата). Для этого проверяем являются ли наши даты «январь 2010 г.», датой, если нет, то переводим в дату и используя функцию Excel =месяц(дата), получаем номера месяцев:

Получаем ряд с пронумерованными месяцами:

Далее усредняем отклонения ряда от средней для каждого месяца, получаем 12 значений аддитивной сезонности.

Для этого используем формулы Excel:

СУММЕСЛИ($D$7:$AY$7 (диапазон с номерами месяцев);D7 (номер месяца, для которого мы рассчитываем сезонность);$D$6:$AY$6 (разность между рядом и средней))

СЧЁТЕСЛИ($D$7:$AY$7(диапазон с номерами месяцев);D7(номер месяца, для которого мы рассчитываем сезонность))

Обязательно фиксируем ссылки на диапазоны с «Номерами месяцев» и «разность между рядом и средней». Подробнее об этом в статье » Как зафиксировать ссылку в Excel»

Подробнее о формулах Excel СУММЕСЛИ и СЧЁТЕСЛИ читайте с статье «Формулы Excel «СУММЕСЛИ» и «СЧЕТЕСЛИ» при расчете сезонности»

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

mad 6

4. Продлим значения ряда с помощью скользящей средней к 4-м месяцам и скорректируем её аддитивной сезонностью.

Скорректируем скользящую рассчитанной аддитивной сезонностью.

Для этого к прогнозному среднему прибавим аддитивную сезонность. Сезонность для каждого месяца подтянем с помощью функции Excel ГПР.

Подробнее об этом читайте с статье «ГПР в Excel на примере скользящей средней».

Прогноз = средние продажи за последние 4 месяца + сезонность:

=СРЗНАЧ(AV4:AY4(средние продажи за 4 последних месяца))+ГПР(AZ3 (искомый номер месяца);$D$8:$O$9 (зафиксированная ссылка на таблицу с сезонностью);2 (номер строки);0)

5. Рассчитаем модель прогноза для каждого момента времени t.

К скользящей средней прибавим аддитивную сезонность начиная с 5 периода:

6. Рассчитаем значение ошибки для каждого месяца.

Для этого из объема продаж вычтем значение прогнозной модели:

7. Определим среднее абсолютное отклонение.

Для каждого момента времени t рассчитаем ошибку по модулю с помощью формулы Excel =ABS(H11 (ссылка на ошибку)):

Среднее абсолютное отклонение равно средней ошибке по модулю:

Среднее абсолютное отклонение для модели скользящей средней к 4-м месяцам с аддитивной сезонностью у нас равно 55 475

Теперь рассчитаем прогноз с помощью «Логарифмического тренда с сезонностью».

Выделим логарифмический тренд;
Рассчитаем сезонность;
Рассчитаем значение модели;
Рассчитаем ошибку прогноза и Среднее абсолютное отклонение.

Скачайте файл с примером

1. Выделим логарифмический тренд.

О всех возможных способах выделения логарифмического тренда в Excel вы можете узнать в нашей статье «5 способов расчета логарифмического тренда в Excel. + О логарифмическом тренде и его применении».

Рассчитаем значения тренда с помощью функции =ПРЕДСКАЗ(LN(D2(номер периода));$D$4:$AY$4 (зафиксированная ссылка на диапазон с объемами продаж);LN($D$2:$AY$2 (зафиксированная ссылка на диапазон с номерами периодов)))

2. Рассчитываем отклонения объемов продаж от тренда (объем продаж делим на значения тренда):

3. Определяем сезонность с помощью формул Excel =СУММЕСЛИ() и =СЧЁТЕСЛИ()

Подробнее о формулах Excel СУММЕСЛИ и СЧЁТЕСЛИ читайте с статье «Формулы Excel «СУММЕСЛИ» и «СЧЕТЕСЛИ» при расчете сезонности»:

Т.к. полученная сезонность в среднем равна 1, то нормирующий коэффициент вводить не нужно, и среднее отклонение у нас будет равно сезонности по месяцам.

4. Определим значения модели прогноза для каждого момента времени t, для этого значения тренда умножим на сезонность. Сезонность подтянем с помощью функции ГПР (см. статью «Функция ГПР в Excel»):

5. Определим ошибку прогноза для каждого момента времени t. Для этого из объема продаж вычтем значение модели прогноза для каждого момента времени t:

6. Рассчитаем ошибку по модулю с помощью функции =ABS(D21″ссылка на ошибку»):

7. Получим среднее абсолютное отклонение по модулю — среднее значение ошибки по модулю:

Среднее абсолютное отклонение для модели «Логарифмического тренда с сезонностью» у нас равно 70 412

Оценим эффективность использования в рамках года одной модели относительно другой.

Среднее абсолютное отклонение для модели

«Логарифмического тренда с сезонностью» = 70 412 руб.
«Скользящей средней к 4-м месяцам с аддитивной сезонностью» = 55 475 руб.

Итак модель скользящей средней делает более точный прогноз по сравнение с логарифмическим трендом для этого ряда в месяц на 14 937 руб. = 70 412 руб. — 55 475 руб.

В результате для нас это означает экономию оборотных средств на 14 937 руб. в месяц на обслуживание модели и 179 242 руб. в год, т.е. 14 937 руб. в месяц = 14 937 руб. * 12 месяцев =179 242 руб.

Т.е. в год мы получаем дополнительные оборотные средства в размере 179 242 руб.

Вот так вот за счет оценки точности прогноза и использования модели, которая дает меньшую ошибку прогноза, вы получаете дополнительные оборотные средства — 179 242 руб. в год.

Скачайте файл с примером

Коллеги, эти 2 модели я выбрал наугад, давайте теперь оценим модель, которую автоматически подберет Forecast4AC PRO. И оценим, какой эффект в год нам это даст.

В настройках программы во вкладке «Доп. возможности» ставим галочку «MAD — Среднее абсолютное отклонение» и отключаем модели экспоненциального сглаживания (т.к. они дают ошибку для данного ряда больше чем скользящая средняя и трендовые модели):

Сохраняем и рассчитываем прогноз с помощью Forecast4AC PRO с автоматическим выбором модели.

Автоматически программа выбрала модель Средняя за 2 предыдущих периода + Сезонность относительно средней в абсолютной величине (т.е. аддитивная сезонность). Среднее абсолютное отклонение для этой модели у нас получилось равным 39 882 руб.

Экономия оборотных средств модели Forecast4AC PRO относительно модели скользящей к 4-м месяцам в год руб.:	187 115 руб.
Экономия оборотных средств модели Forecast4AC PRO относительно модели Логарифмический тренд с сезонностью в год руб.:	366 357 руб.

Оцените точность моделей прогнозирования, которые вы используете сейчас, рассчитайте проноз с помощью Forecast4AC PRO и оцените сумму оборотных средств, которую вы можете сэкономить за счет использования нашей программы.

Точных вам прогнозов!

Присоединяйтесь к нам!

Скачивайте бесплатные приложения для прогнозирования и бизнес-анализа:

Novo Forecast Lite — автоматический расчет прогноза в Excel.
4analytics — ABC-XYZ-анализ и анализ выбросов в Excel.
Qlik Sense Desktop и QlikView Personal Edition — BI-системы для анализа и визуализации данных.

Тестируйте возможности платных решений:

Novo Forecast PRO — прогнозирование в Excel для больших массивов данных.

Получите 10 рекомендаций по повышению точности прогнозов до 90% и выше.

Зарегистрируйтесь и скачайте решения

Статья полезная? Поделитесь с друзьями

Источник

Оценка ошибки прогнозирования временного ряда

Работая с научными публикациями, сталкиваюсь с различными показателями ошибок прогнозирования временных рядов. Среди всех встречающихся оценок ошибки прогнозирования стоит отметить две, которые в настоящее время, являются самыми популярными: MAE и MAPE.
Пусть ошибка есть разность:
,
где Z(t) – фактическое значение временного ряда, а – прогнозное.
Тогда формулы для оценок ошибки прогнозирования временных рядов для N отчетов можно записать в следующем виде.

MAPE – средняя абсолютная ошибка в процентах

Данная оценка применяется для временных рядов, фактические значения которых значительно больше 1. Например, оценки ошибки прогнозирования энергопотребления почти во всех статьях приводятся как значения MAPE.

Если же фактические значения временного ряда близки к 0, то в знаменателе окажется очень маленькое число, что сделает значение MAPE близким к бесконечности – это не совсем корректно. Например, фактическая цена РСВ = 0.01 руб/МВт.ч, a прогнозная = 10 руб/МВт.ч, тогда MAPE = (0.01 – 10)/0.01 = 999%, хотя в действительности мы не так уж сильно ошиблись, всего на 10 руб/МВт.ч. Для рядов, содержащих значения близкие к нулю, применяют следующую оценку ошибки прогноза.

MAE – средняя абсолютная ошибка

Для оценки ошибки прогнозирования цен РСВ и индикатора БР корректнее использовать MAE.

После того, как получены значения для MAPE и/или MAE, то в работах обычно пишут: «Прогнозирование временного ряда энергопотребления с часовым разрешение проводилось на интервале с 01.01.2001 до 31.12.2001 (общее количество отсчетов N ~ 8500). Для данного прогноза значение MAPE = 1.5%». При этом, просматривая статьи, можно сложить общее впечатление об ошибки прогнозирования энергопотребления, для которого MAPE обычно колеблется от 1 до 5%; или ошибки прогнозирования цен на электроэнергию, для которого MAPE колеблется от 5 до 15% в зависимости от периода и рынка. Получив значение MAPE для собственного прогноза, вы можете оценить, насколько здорово у вас получается прогнозировать.

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

MSE – среднеквадратичная ошибка

RMSE – квадратный корень из среднеквадратичной ошибки

ME – средняя ошибка

SD – стандартное отклонение

, где ME – есть средняя ошибка, определенная по формуле выше.

Связь точности и ошибки прогнозирования

Точность прогнозирования есть понятие прямо противоположное ошибке прогнозирования. Если ошибка прогнозирования велика, то точность мала и наоборот, если ошибка прогнозирования мала, то точность велика. По сути дела оценка ошибки прогноза MAPE есть обратная величина для точности прогнозирования — зависимость здесь простая.

Точность прогноза в % = 100% – MAPE

Величину точности оценивать не принято, говоря о прогнозировании всегда оценивают, то есть определяют значение именно ошибки прогноза, то есть величину MAPE и/или MAE. Однако нужно понимать, что если MAPE = 5%, то точность прогнозирования = 95%. Говоря о высокой точности, мы всегда говорим о низкой ошибки прогноза и в этой области недопонимания быть не должно. Вы практически не найдете материалов о прогнозировании, в которых приведены оценки именно точности прогноза, хотя с точки зрения здравого маркетинга корректней говорить именно о высокой точности. В рекламных статьях всегда будет написано о высокой точности.

При этом величина MAPE является количественной оценкой именно ошибки, и эта величина нам ясно говорит и о точности прогнозирования, исходя из приведенной выше простой формулы. Таким образом, оценивая ошибку, мы всегда оцениваем точность прогнозирования.

Источник

Доброго дня! Мы продолжаем наш цикл статей открытого курса по машинному обучению и сегодня поговорим о временных рядах.

Посмотрим на то, как с ними работать в Python, какие возможные методы и модели можно использовать для прогнозирования; что такое двойное и тройное экспоненциальное взвешивание; что делать, если стационарность — это не про вас; как построить SARIMA и не умереть; и как прогнозировать xgboost-ом. И всё это будем применять к примеру из суровой реальности.

UPD 01.2022: С февраля 2022 г. ML-курс ODS на русском возрождается под руководством Петра Ермакова couatl. Для русскоязычной аудитории это предпочтительный вариант (c этими статьями на Хабре – в подкрепление), англоговорящим рекомендуется mlcourse.ai в режиме самостоятельного прохождения.

Видеозапись лекции по мотивам этой статьи в рамках второго запуска открытого курса (сентябрь-ноябрь 2017).

План этой статьи:

Движемся, сглаживаем и оцениваем
- Rolling window estimations
- Экспоненциальное сглаживание, модель Хольта-Винтерса
- Кросс-валидация на временных рядах, подбор параметров
Эконометрический подход
- Стационарность, единичные корни
- Избавляемся от нестационарности и строим SARIMA
Линейные и не очень модели на временных рядах
- Извлечение признаков (Feature extraction)
- Линейная регрессия vs XGBoost
Домашнее задание
Полезные ресурсы

Введение

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

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

Движемся, сглаживаем и оцениваем

Небольшое определение временного ряда:

Временной ряд – это последовательность значений, описывающих протекающий во времени процесс, измеренных в последовательные моменты времени, обычно через равные промежутки

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

Импортируем нужные библиотеки. В основном нам понадобится модуль statsmodels, в котором реализованы многочисленные методы статистического моделирования, в том числе для временных рядов. Для поклонников R, пересевших на питон, он может показаться очень родным, так как поддерживает написание формулировок моделей в стиле ‘Wage ~ Age + Education’.

import sys
import warnings
warnings.filterwarnings('ignore')
from tqdm import tqdm

import pandas as pd
import numpy as np
from sklearn.metrics import mean_absolute_error, mean_squared_error

import statsmodels.formula.api as smf
import statsmodels.tsa.api as smt
import statsmodels.api as sm
import scipy.stats as scs
from scipy.optimize import minimize

import matplotlib.pyplot as plt

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

Код для отрисовки графика

from plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot
from plotly import graph_objs as go
init_notebook_mode(connected = True)

def plotly_df(df, title = ''):
    data = []

    for column in df.columns:
        trace = go.Scatter(
            x = df.index,
            y = df[column],
            mode = 'lines',
            name = column
        )
        data.append(trace)

    layout = dict(title = title)
    fig = dict(data = data, layout = layout)
    iplot(fig, show_link=False)

dataset = pd.read_csv('hour_online.csv', index_col=['Time'], parse_dates=['Time'])
plotly_df(dataset, title = "Online users")

Rolling window estimations

Начнем моделирование с наивного предположения — «завтра будет, как вчера», но вместо модели вида $hat{y}_{t} = y_{t-1}$ будем считать, что будущее значение переменной зависит от среднего $n$ её предыдущих значений, а значит, воспользуемся скользящей средней.

$hat{y}_{t} = frac{1}{k} displaystylesum^{k-1}_{n=0} y_{t-n}$

Реализуем эту же функцию в питоне и посмотрим на прогноз, построенный по последнему наблюдаемому дню (24 часа)

def moving_average(series, n):
    return np.average(series[-n:])

moving_average(dataset.Users, 24)

Out: 29858.333333333332

К сожалению, такой прогноз долгосрочным сделать не удастся — для получения предсказания на шаг вперед предыдущее значение должно быть фактически наблюдаемой величиной. Зато у скользящей средней есть другое применение — сглаживание исходного ряда для выявления трендов. В пандасе есть готовая реализация — DataFrame.rolling(window).mean(). Чем больше зададим ширину интервала — тем более сглаженным окажется тренд. В случае, если данные сильно зашумлены, что особенно часто встречается, например, в финансовых показателях, такая процедура может помочь с определением общих паттернов.

Для нашего ряда тренды и так вполне очевидны, но если сгладить по дням, становится лучше видна динамика онлайна по будням и выходным (выходные — время поиграть), а недельное сглаживание хорошо отражает общие изменения, связанные с резким ростом числа активных игроков в феврале и последующим снижением в марте.

Код для отрисовки графика

def plotMovingAverage(series, n):

    """
    series - dataframe with timeseries
    n - rolling window size 

    """

    rolling_mean = series.rolling(window=n).mean()

    # При желании, можно строить и доверительные интервалы для сглаженных значений
    #rolling_std =  series.rolling(window=n).std()
    #upper_bond = rolling_mean+1.96*rolling_std
    #lower_bond = rolling_mean-1.96*rolling_std

    plt.figure(figsize=(15,5))
    plt.title("Moving averagen window size = {}".format(n))
    plt.plot(rolling_mean, "g", label="Rolling mean trend")

    #plt.plot(upper_bond, "r--", label="Upper Bond / Lower Bond")
    #plt.plot(lower_bond, "r--")
    plt.plot(dataset[n:], label="Actual values")
    plt.legend(loc="upper left")
    plt.grid(True)

plotMovingAverage(dataset, 24) # сглаживаем по дням
plotMovingAverage(dataset, 24*7) # сглаживаем по неделям

Модификацией простой скользящей средней является взвешенная средняя, внутри которой наблюдениям придаются различные веса, в сумме дающие единицу, при этом обычно последним наблюдениям присваивается больший вес.

$hat{y}_{t} = displaystylesum^{k}_{n=1} omega_n y_{t+1-n}$

def weighted_average(series, weights):
    result = 0.0
    weights.reverse()
    for n in range(len(weights)):
        result += series[-n-1] * weights[n]
    return result

weighted_average(dataset.Users, [0.6, 0.2, 0.1, 0.07, 0.03])

Out: 35967.550000000003

Экспоненциальное сглаживание, модель Хольта-Винтерса

Простое экспоненциальное сглаживание

А теперь посмотрим, что произойдёт, если вместо взвешивания последних $n$ значений ряда мы начнем взвешивать все доступные наблюдения, при этом экспоненциально уменьшая веса по мере углубления в исторические данные. В этом нам поможет формула простого экспоненциального сглаживания:

$hat{y}_{t} = alpha cdot y_t + (1-alpha) cdot hat y_{t-1}$

Здесь модельное значение представляет собой средневзвешенную между текущим истинным и предыдущим модельным значениями. Вес $alpha$ называется сглаживающим фактором. Он определяет, как быстро мы будем «забывать» последнее доступное истинное наблюдение. Чем меньше $alpha$ , тем больше влияния оказывают предыдущие модельные значения, и тем сильнее сглаживается ряд.

Экспоненциальность скрывается в рекурсивности функции — каждый раз мы умножаем на предыдущее модельное значение, которое, в свою очередь, также содержало в себе , и так до самого начала.

def exponential_smoothing(series, alpha):
    result = [series[0]] # first value is same as series
    for n in range(1, len(series)):
        result.append(alpha * series[n] + (1 - alpha) * result[n-1])
    return result

Код для отрисовки графика

with plt.style.context('seaborn-white'):    
    plt.figure(figsize=(20, 8))
    for alpha in [0.3, 0.05]:
        plt.plot(exponential_smoothing(dataset.Users, alpha), label="Alpha {}".format(alpha))
    plt.plot(dataset.Users.values, "c", label = "Actual")
    plt.legend(loc="best")
    plt.axis('tight')
    plt.title("Exponential Smoothing")
    plt.grid(True)

Двойное экспоненциальное сглаживание

До сих пор мы могли получить от наших методов в лучшем случае прогноз лишь на одну точку вперёд (и ещё красиво сгладить ряд), это здорово, но недостаточно, поэтому переходим к расширению экспоненциального сглаживания, которое позволит строить прогноз сразу на две точки вперед (и тоже красиво сглаживать ряд).

В этом нам поможет разбиение ряда на две составляющие — уровень (level, intercept) $ell$ и тренд $b$ (trend, slope). Уровень, или ожидаемое значение ряда, мы предсказывали при помощи предыдущих методов, а теперь такое же экспоненциальное сглаживание применим к тренду, наивно или не очень полагая, что будущее направление изменения ряда зависит от взвешенных предыдущих изменений.

$ell_x = alpha y_x + (1-alpha)(ell_{x-1} + b_{x-1})\ b_x = beta(ell_x - ell_{x-1}) + (1-beta)b_{x-1}\ hat{y}_{x+1} = ell_x + b_x$

В результате получаем набор функций. Первая описывает уровень — он, как и прежде, зависит от текущего значения ряда, а второе слагаемое теперь разбивается на предыдущее значение уровня и тренда. Вторая отвечает за тренд — он зависит от изменения уровня на текущем шаге, и от предыдущего значения тренда. Здесь в роли веса в экспоненциальном сглаживании выступает коэффициент $beta$ . Наконец, итоговое предсказание представляет собой сумму модельных значений уровня и тренда.

def double_exponential_smoothing(series, alpha, beta):
    result = [series[0]]
    for n in range(1, len(series)+1):
        if n == 1:
            level, trend = series[0], series[1] - series[0]
        if n >= len(series): # прогнозируем
            value = result[-1]
        else:
            value = series[n]
        last_level, level = level, alpha*value + (1-alpha)*(level+trend)
        trend = beta*(level-last_level) + (1-beta)*trend
        result.append(level+trend)
    return result

Код для отрисовки графика

with plt.style.context('seaborn-white'):    
    plt.figure(figsize=(20, 8))
    for alpha in [0.9, 0.02]:
        for beta in [0.9, 0.02]:
            plt.plot(double_exponential_smoothing(dataset.Users, alpha, beta), label="Alpha {}, beta {}".format(alpha, beta))
    plt.plot(dataset.Users.values, label = "Actual")
    plt.legend(loc="best")
    plt.axis('tight')
    plt.title("Double Exponential Smoothing")
    plt.grid(True)

Теперь настраивать пришлось уже два параметра — $alpha$ и $beta$ . Первый отвечает за сглаживание ряда вокруг тренда, второй — за сглаживание самого тренда. Чем выше значения, тем больший вес будет отдаваться последним наблюдениям и тем менее сглаженным окажется модельный ряд. Комбинации параметров могут выдавать достаточно причудливые результаты, особенно если задавать их руками. А о не ручном подборе параметров расскажу чуть ниже, сразу после тройного экспоненциального сглаживания.

Тройное экспоненциальное сглаживание a.k.a. Holt-Winters

Итак, успешно добрались до следующего варианта экспоненциального сглаживания, на сей раз тройного.

Идея этого метода заключается в добавлении еще одной, третьей, компоненты — сезонности. Соответственно, метод применим только в случае, если ряд этой сезонностью не обделён, что в нашем случае верно. Сезонная компонента в модели будет объяснять повторяющиеся колебания вокруг уровня и тренда, а характеризоваться она будет длиной сезона — периодом, после которого начинаются повторения колебаний. Для каждого наблюдения в сезоне формируется своя компонента, например, если длина сезона составляет 7 (например, недельная сезонность), то получим 7 сезонных компонент, по штуке на каждый из дней недели.

Получаем новую систему:

$ell_x = alpha(y_x - s_{x-L}) + (1-alpha)(ell_{x-1} + b_{x-1})\ b_x = beta(ell_x - ell_{x-1}) + (1-beta)b_{x-1}\ s_x = gamma(y_x - ell_x) + (1-gamma)s_{x-L}\ hat{y}_{x+m} = ell_x + mb_x + s_{x-L+1+(m-1)modL}$

Уровень теперь зависит от текущего значения ряда за вычетом соответствующей сезонной компоненты, тренд остаётся без изменений, а сезонная компонента зависит от текущего значения ряда за вычетом уровня и от предыдущего значения компоненты. При этом компоненты сглаживаются через все доступные сезоны, например, если это компонента, отвечающая за понедельник, то и усредняться она будет только с другими понедельниками. Подробнее про работу усреднений и оценку начальных значений тренда и сезонных компонент можно почитать здесь. Теперь, имея сезонную компоненту, мы можем предсказывать уже не на один, и даже не на два, а на произвольные $m$ шагов вперёд, что не может не радовать.

Ниже приведен код для построения модели тройного экспоненциального сглаживания, также известного по фамилиям её создателей — Чарльза Хольта и его студента Питера Винтерса.
Дополнительно в модель включен метод Брутлага для построения доверительных интервалов:

$hat y_{max_x}=ell_{x−1}+b_{x−1}+s_{x−T}+m⋅d_{t−T}\ hat y_{min_x}=ell_{x−1}+b_{x−1}+s_{x−T}-m⋅d_{t−T}\ d_t=gamma∣y_t−hat y_t∣+(1−gamma)d_{t−T},$

где $T$ — длина сезона, $d$ — предсказанное отклонение, а остальные параметры берутся из тройного сглаживани. Подробнее о методе и о его применении к поиску аномалий во временных рядах можно прочесть здесь

Код для модели Хольта-Винтерса

class HoltWinters:

    """
    Модель Хольта-Винтерса с методом Брутлага для детектирования аномалий
    https://fedcsis.org/proceedings/2012/pliks/118.pdf

    # series - исходный временной ряд
    # slen - длина сезона
    # alpha, beta, gamma - коэффициенты модели Хольта-Винтерса
    # n_preds - горизонт предсказаний
    # scaling_factor - задаёт ширину доверительного интервала по Брутлагу (обычно принимает значения от 2 до 3)

    """

    def __init__(self, series, slen, alpha, beta, gamma, n_preds, scaling_factor=1.96):
        self.series = series
        self.slen = slen
        self.alpha = alpha
        self.beta = beta
        self.gamma = gamma
        self.n_preds = n_preds
        self.scaling_factor = scaling_factor

    def initial_trend(self):
        sum = 0.0
        for i in range(self.slen):
            sum += float(self.series[i+self.slen] - self.series[i]) / self.slen
        return sum / self.slen  

    def initial_seasonal_components(self):
        seasonals = {}
        season_averages = []
        n_seasons = int(len(self.series)/self.slen)
        # вычисляем сезонные средние
        for j in range(n_seasons):
            season_averages.append(sum(self.series[self.slen*j:self.slen*j+self.slen])/float(self.slen))
        # вычисляем начальные значения
        for i in range(self.slen):
            sum_of_vals_over_avg = 0.0
            for j in range(n_seasons):
                sum_of_vals_over_avg += self.series[self.slen*j+i]-season_averages[j]
            seasonals[i] = sum_of_vals_over_avg/n_seasons
        return seasonals   

    def triple_exponential_smoothing(self):
        self.result = []
        self.Smooth = []
        self.Season = []
        self.Trend = []
        self.PredictedDeviation = []
        self.UpperBond = []
        self.LowerBond = []

        seasonals = self.initial_seasonal_components()

        for i in range(len(self.series)+self.n_preds):
            if i == 0: # инициализируем значения компонент
                smooth = self.series[0]
                trend = self.initial_trend()
                self.result.append(self.series[0])
                self.Smooth.append(smooth)
                self.Trend.append(trend)
                self.Season.append(seasonals[i%self.slen])

                self.PredictedDeviation.append(0)

                self.UpperBond.append(self.result[0] + 
                                      self.scaling_factor * 
                                      self.PredictedDeviation[0])

                self.LowerBond.append(self.result[0] - 
                                      self.scaling_factor * 
                                      self.PredictedDeviation[0])

                continue
            if i >= len(self.series): # прогнозируем
                m = i - len(self.series) + 1
                self.result.append((smooth + m*trend) + seasonals[i%self.slen])

                # во время прогноза с каждым шагом увеличиваем неопределенность
                self.PredictedDeviation.append(self.PredictedDeviation[-1]*1.01) 

            else:
                val = self.series[i]
                last_smooth, smooth = smooth, self.alpha*(val-seasonals[i%self.slen]) + (1-self.alpha)*(smooth+trend)
                trend = self.beta * (smooth-last_smooth) + (1-self.beta)*trend
                seasonals[i%self.slen] = self.gamma*(val-smooth) + (1-self.gamma)*seasonals[i%self.slen]
                self.result.append(smooth+trend+seasonals[i%self.slen])

                # Отклонение рассчитывается в соответствии с алгоритмом Брутлага
                self.PredictedDeviation.append(self.gamma * np.abs(self.series[i] - self.result[i]) 
                                               + (1-self.gamma)*self.PredictedDeviation[-1])

            self.UpperBond.append(self.result[-1] + 
                                  self.scaling_factor * 
                                  self.PredictedDeviation[-1])

            self.LowerBond.append(self.result[-1] - 
                                  self.scaling_factor * 
                                  self.PredictedDeviation[-1])

            self.Smooth.append(smooth)
            self.Trend.append(trend)
            self.Season.append(seasonals[i % self.slen])

Кросс-валидация на временных рядах, подбор параметров

Перед тем, как построить модель, поговорим, наконец, о не ручной оценке параметров для моделей.

Ничего необычного здесь нет, по-прежнему сначала необходимо выбрать подходящуюю для данной задачи функцию потерь: RMSE, MAE, MAPE и др., которая будет следить за качеством подгонки модели под исходные данные. Затем будем оценивать на кросс-валидации значение функции потерь при данных параметрах модели, искать градиент, менять в соответствии с ним параметры и бодро опускаться в сторону глобального минимума ошибки.

Небольшая загвоздка возникает только в кросс-валидации. Проблема состоит в том, что временной ряд имеет, как ни парадоксально, временную структуру, и случайно перемешивать в фолдах значения всего ряда без сохранения этой структуры нельзя, иначе в процессе потеряются все взаимосвязи наблюдений друг с другом. Поэтому придется использовать чуть более хитрый способ для оптимизации параметров, официального названия которому я так и не нашел, но на сайте CrossValidated, где можно найти ответы на всё, кроме главного вопроса Жизни, Вселенной и Всего Остального, предлагают название «cross-validation on a rolling basis», что не дословно можно перевести как кросс-валидация на скользящем окне.

Суть достаточно проста — начинаем обучать модель на небольшом отрезке временного ряда, от начала до некоторого $t$ , делаем прогноз на $t+n$ шагов вперед и считаем ошибку. Далее расширяем обучающую выборку до $t+n$ значения и прогнозируем с $t+n$ до $t+2*n$ , так продолжаем двигать тестовый отрезок ряда до тех пор, пока не упрёмся в последнее доступное наблюдение. В итоге получим столько фолдов, сколько $n$ уместится в промежуток между изначальным обучающим отрезком и всей длиной ряда.

Код для кросс-валидации на временном ряду

from sklearn.model_selection import TimeSeriesSplit

def timeseriesCVscore(x):
    # вектор ошибок
    errors = []

    values = data.values
    alpha, beta, gamma = x

    # задаём число фолдов для кросс-валидации
    tscv = TimeSeriesSplit(n_splits=3) 

    # идем по фолдам, на каждом обучаем модель, строим прогноз на отложенной выборке и считаем ошибку
    for train, test in tscv.split(values):

        model = HoltWinters(series=values[train], slen = 24*7, alpha=alpha, beta=beta, gamma=gamma, n_preds=len(test))
        model.triple_exponential_smoothing()

        predictions = model.result[-len(test):]
        actual = values[test]
        error = mean_squared_error(predictions, actual)
        errors.append(error)

    # Возвращаем средний квадрат ошибки по вектору ошибок 
    return np.mean(np.array(errors))

Значение длины сезона 24*7 возникло не случайно — в исходном ряде отчетливо видна дневная сезонность, (отсюда 24), и недельная — по будням ниже, на выходных — выше, (отсюда 7), суммарно сезонных компонент получится 24*7.

В модели Хольта-Винтерса, как и в остальных моделях экспоненциального сглаживания, есть ограничение на величину сглаживающих параметров — каждый из них может принимать значения от 0 до 1, поэтому для минимизации функции потерь нужно выбирать алгоритм, поддерживающий ограничения на параметры, в данном случае — Truncated Newton conjugate gradient.

%%time
data = dataset.Users[:-500] # отложим часть данных для тестирования

# инициализируем значения параметров
x = [0, 0, 0] 

# Минимизируем функцию потерь с ограничениями на параметры
opt = minimize(timeseriesCVscore, x0=x, method="TNC", bounds = ((0, 1), (0, 1), (0, 1)))

# Из оптимизатора берем оптимальное значение параметров
alpha_final, beta_final, gamma_final = opt.x
print(alpha_final, beta_final, gamma_final)

Out: (0.0066342670643441681, 0.0, 0.046765204289672901)

Передадим полученные оптимальные значения коэффициентов $alpha$ , $beta$ и $gamma$ и построим прогноз на 5 дней вперёд (128 часов)

# Передаем оптимальные значения модели, 
data = dataset.Users
model = HoltWinters(data[:-128], slen = 24*7, alpha = alpha_final, beta = beta_final, gamma = gamma_final, n_preds = 128, scaling_factor = 2.56)
model.triple_exponential_smoothing()

Код для отрисовки графика

def plotHoltWinters():
    Anomalies = np.array([np.NaN]*len(data))
    Anomalies[data.values<model.LowerBond] = data.values[data.values<model.LowerBond]
    plt.figure(figsize=(25, 10))
    plt.plot(model.result, label = "Model")
    plt.plot(model.UpperBond, "r--", alpha=0.5, label = "Up/Low confidence")
    plt.plot(model.LowerBond, "r--", alpha=0.5)
    plt.fill_between(x=range(0,len(model.result)), y1=model.UpperBond, y2=model.LowerBond, alpha=0.5, color = "grey")
    plt.plot(data.values, label = "Actual")
    plt.plot(Anomalies, "o", markersize=10, label = "Anomalies")
    plt.axvspan(len(data)-128, len(data), alpha=0.5, color='lightgrey')
    plt.grid(True)
    plt.axis('tight')
    plt.legend(loc="best", fontsize=13);

plotHoltWinters()

Судя по графику, модель неплохо описала исходный временной ряд, уловив недельную и дневную сезонность, и даже смогла поймать аномальные снижения, вышедшие за пределы доверительных интервалов. Если посмотреть на смоделированное отклонение, хорошо видно, что модель достаточно резко регирует на значительные изменения в структуре ряда, но при этом быстро возвращает дисперсию к обычным значениям, «забывая» прошлое. Такая особенность позволяет неплохо и без значительных затрат на подготовку-обучение модели настроить систему по детектированию аномалий даже в достаточно шумных рядах.

Эконометрический подход

Стационарность, единичные корни

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

Временной ряд справа не является стационарным, так как его матожидание со временем растёт

Здесь не повезло с дисперсией — разброс значений ряда существенно варьируется в зависимости от периода

Наконец, на последнем графике видно, что значения ряда внезапно становятся ближе друг ко другу, образуя некоторый кластер, а в результате получаем непостоянство ковариаций

Почему стационарность так важна? По стационарному ряду просто строить прогноз, так как мы полагаем, что его будущие статистические характеристики не будут отличаться от наблюдаемых текущих. Большинство моделей временных рядов так или иначе моделируют и предсказывают эти характеристики (например, матожидание или дисперсию), поэтому в случае нестационарности исходного ряда предсказания окажутся неверными. К сожалению, большинство временных рядов, с которыми приходится сталкиваться за пределыми учебных материалов, стационарными не являются, но с этим можно (и нужно) бороться.

Чтобы бороться с нестационарностью, нужно узнать её в лицо, потому посмотрим, как её детектировать. Для этого обратимся к белому шуму и случайному блужданию, чтобы выяснить как попасть из одного в другое бесплатно и без смс.

График белого шума:

white_noise = np.random.normal(size=1000)
with plt.style.context('bmh'):  
    plt.figure(figsize=(15, 5))
    plt.plot(white_noise)

Итак, процесс, порожденный стандартным нормальным распределением, стационарен, колеблется вокруг нуля с отклонением в 1. Теперь на основании него сгенерируем новый процесс, в котором каждое последующее значение будет зависеть от предыдущего: $x_t = rho x_{t-1} + e_t$

Код для отрисовки графиков

def plotProcess(n_samples=1000, rho=0):
    x = w = np.random.normal(size=n_samples)
    for t in range(n_samples):
        x[t] = rho * x[t-1] + w[t]

    with plt.style.context('bmh'):  
        plt.figure(figsize=(10, 3))
        plt.plot(x)
        plt.title("Rho {}n Dickey-Fuller p-value: {}".format(rho, round(sm.tsa.stattools.adfuller(x)[1], 3)))

for rho in [0, 0.6, 0.9, 1]:
    plotProcess(rho=rho)

На первом графике получился точно такой же стационарный белый шум, который строился раньше. На втором значение $rho$ увеличилось до 0.6, в результате чего на графике стали появляться более широкие циклы, но в целом стационарным он быть пока не перестал. Третий график всё сильнее отклоняется от нулевого среднего значения, но всё ещё колеблется вокруг него. Наконец, значение $rho$ равное единице дало процесс случайного блуждания — ряд не стационарен.

Происходит это из-за того, что при достижении критической единицы, ряд $x_t = rho x_{t-1} + e_t$ перестаёт возвращаться к своему среднему значению. Если вычесть из левой и правой части $x_{t-1}$ , то получим $x_t - x_{t-1} = (rho - 1) x_{t-1} + e_t$ , где выражение слева — первые разности. Если $rho=1$ , то первые разности дадут стационарный белый шум $e_t$ . Этот факт лёг в основу теста Дики-Фуллера на стационарность ряда (наличие единичного корня). Если из нестационарного ряда первыми разностями удаётся получить стационарный, то он называется интегрированным первого порядка. Нулевая гипотеза теста — ряд не стационарен, отвергалась на первых трех графиках, и принялась на последнем. Стоит сказать, что не всегда для получения стационарного ряда хватает первых разностей, так как процесс может быть интегрированным с более высоким порядком (иметь несколько единичных корней), для проверки таких случаев используют расширенный тест Дики-Фуллера, проверяющий сразу несколько лагов.

Бороться с нестационарностью можно множеством способов — разностями различного порядка, выделением тренда и сезонности, сглаживаниями и преобразованиями, например, Бокса-Кокса или логарифмированием.

Избавляемся от нестационарности и строим SARIMA

Попробуем теперь построить ARIMA модель для онлайна игроков, пройдя все круги ада стадии приведения ряда к стационарному виду. Про саму модель уже не раз писали на хабре — Построение модели SARIMA с помощью Python+R, Анализ временных рядов с помощью python, поэтому подробно останавливаться на ней не буду.

Код для отрисовки графиков

def tsplot(y, lags=None, figsize=(12, 7), style='bmh'):
    if not isinstance(y, pd.Series):
        y = pd.Series(y)
    with plt.style.context(style):    
        fig = plt.figure(figsize=figsize)
        layout = (2, 2)
        ts_ax = plt.subplot2grid(layout, (0, 0), colspan=2)
        acf_ax = plt.subplot2grid(layout, (1, 0))
        pacf_ax = plt.subplot2grid(layout, (1, 1))

        y.plot(ax=ts_ax)
        ts_ax.set_title('Time Series Analysis Plots')
        smt.graphics.plot_acf(y, lags=lags, ax=acf_ax, alpha=0.5)
        smt.graphics.plot_pacf(y, lags=lags, ax=pacf_ax, alpha=0.5)

        print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(y)[1])

        plt.tight_layout()
    return 

tsplot(dataset.Users, lags=30)

Out: Критерий Дики-Фуллера: p=0.190189

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

def invboxcox(y,lmbda):
    # обрабтное преобразование Бокса-Кокса
    if lmbda == 0:
        return(np.exp(y))
    else:
        return(np.exp(np.log(lmbda*y+1)/lmbda))

data = dataset.copy()
data['Users_box'], lmbda = scs.boxcox(data.Users+1) # прибавляем единицу, так как в исходном ряде есть нули
tsplot(data.Users_box, lags=30)
print("Оптимальный параметр преобразования Бокса-Кокса: %f" % lmbda)

Out: Критерий Дики-Фуллера: p=0.079760
     Оптимальный параметр преобразования Бокса-Кокса: 0.587270

Уже лучше, однако критерий Дики-Фуллера по-прежнему не отвергает гипотезу о нестационарности ряда. А автокорреляционная функция явно намекает на сезонность в получившемся ряде. Возьмём сезонные разности:

data['Users_box_season'] = data.Users_box - data.Users_box.shift(24*7)
tsplot(data.Users_box_season[24*7:], lags=30)

Out: Критерий Дики-Фуллера: p=0.002571

Критерий Дики-Фуллера теперь отвергает нулевую гипотезу о нестационарности, но автокорреляционная функция всё ещё выглядит нехорошо из-за большого числа значимых лагов. Так как на графике частной автокорреляционной функции значим лишь один лаг, стоит взять еще первые разности, чтобы привести наконец ряд к стационарному виду.

data['Users_box_season_diff'] = data.Users_box_season - data.Users_box_season.shift(1)
tsplot(data.Users_box_season_diff[24*7+1:], lags=30)

Out: Критерий Дики-Фуллера: p=0.000000

Наконец, получили стационарный ряд, по автокорреляционной и частной автокорреляционной функции прикинем параметры для SARIMA модели, на забыв, что предварительно уже сделали первые и сезонные разности.

Начальные приближения Q = 1, P = 4, q = 3, p = 4

ps = range(0, 5)
d=1
qs = range(0, 4)
Ps = range(0, 5)
D=1
Qs = range(0, 1)

from itertools import product

parameters = product(ps, qs, Ps, Qs)
parameters_list = list(parameters)
len(parameters_list)

Out: 100

Код для подбора параметров перебором

%%time
results = []
best_aic = float("inf")

for param in tqdm(parameters_list):
    #try except нужен, потому что на некоторых наборах параметров модель не обучается
    try:
        model=sm.tsa.statespace.SARIMAX(data.Users_box, order=(param[0], d, param[1]), 
                                        seasonal_order=(param[2], D, param[3], 24*7)).fit(disp=-1)
    #выводим параметры, на которых модель не обучается и переходим к следующему набору
    except ValueError:
        print('wrong parameters:', param)
        continue
    aic = model.aic
    #сохраняем лучшую модель, aic, параметры
    if aic < best_aic:
        best_model = model
        best_aic = aic
        best_param = param
    results.append([param, model.aic])

warnings.filterwarnings('default')

result_table = pd.DataFrame(results)
result_table.columns = ['parameters', 'aic']
print(result_table.sort_values(by = 'aic', ascending=True).head())

Лучшие параметры загоняем в модель:

%%time
best_model = sm.tsa.statespace.SARIMAX(data.Users_box, order=(4, d, 3), 
                                        seasonal_order=(4, D, 1, 24)).fit(disp=-1)
print(best_model.summary())

                                 Statespace Model Results                                 
==========================================================================================
Dep. Variable:                          Users_box   No. Observations:                 2625
Model:             SARIMAX(4, 1, 3)x(4, 1, 1, 24)   Log Likelihood              -12547.157
Date:                            Sun, 23 Apr 2017   AIC                          25120.315
Time:                                    02:06:39   BIC                          25196.662
Sample:                                         0   HQIC                         25147.964
                                           - 2625                                         
Covariance Type:                              opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
ar.L1          0.6794      0.108      6.310      0.000       0.468       0.890
ar.L2         -0.0810      0.181     -0.448      0.654      -0.435       0.273
ar.L3          0.3255      0.137      2.371      0.018       0.056       0.595
ar.L4         -0.2154      0.028     -7.693      0.000      -0.270      -0.161
ma.L1         -0.5086      0.106     -4.784      0.000      -0.717      -0.300
ma.L2         -0.0673      0.170     -0.395      0.693      -0.401       0.267
ma.L3         -0.3490      0.117     -2.976      0.003      -0.579      -0.119
ar.S.L24       0.1023      0.012      8.377      0.000       0.078       0.126
ar.S.L48      -0.0686      0.021     -3.219      0.001      -0.110      -0.027
ar.S.L72       0.1971      0.009     21.573      0.000       0.179       0.215
ar.S.L96      -0.1217      0.013     -9.279      0.000      -0.147      -0.096
ma.S.L24      -0.9983      0.045    -22.085      0.000      -1.087      -0.910
sigma2       873.4159     36.206     24.124      0.000     802.454     944.378
===================================================================================
Ljung-Box (Q):                      130.47   Jarque-Bera (JB):           1194707.99
Prob(Q):                              0.00   Prob(JB):                         0.00
Heteroskedasticity (H):               1.40   Skew:                             2.65
Prob(H) (two-sided):                  0.00   Kurtosis:                       107.88
===================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

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

tsplot(best_model.resid[24:], lags=30)

Out: Критерий Дики-Фуллера: p=0.000000

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

Код для построения прогноза и отрисовки графика

data["arima_model"] = invboxcox(best_model.fittedvalues, lmbda)
forecast = invboxcox(best_model.predict(start = data.shape[0], end = data.shape[0]+100), lmbda)
forecast = data.arima_model.append(forecast).values[-500:]
actual = data.Users.values[-400:]
plt.figure(figsize=(15, 7))
plt.plot(forecast, color='r', label="model")
plt.title("SARIMA modeln Mean absolute error {} users".format(round(mean_absolute_error(data.dropna().Users, data.dropna().arima_model))))
plt.plot(actual, label="actual")
plt.legend()
plt.axvspan(len(actual), len(forecast), alpha=0.5, color='lightgrey')
plt.grid(True)

В финале получаем достаточно адекватный прогноз, в среднем модель ошибалась на 1.3 K пользователей, что очень и очень неплохо, однако суммарные затраты на подготовку данных, приведение к стационарности, определение и перебор параметров могут такой точности и не стоить.

Линейные и не очень модели на временных рядах

Снова небольшое лирическое отступление. Часто на работе приходится строить модели, руководствуясь одним основополагающим принципом – быстро, качественно, недорого. Поэтому часть моделей могут банально не подойти для «продакшн-решений», так как либо требуют слишком больших затрат по подготовке данных (например, SARIMA), либо сложно настраиваются (хороший пример – SARIMA), либо требуют частого переобучения на новых данных (опять SARIMA), поэтому зачастую гораздо проще бывает выделить несколько признаков из имеющегося временного ряда и построить по ним обычную линейную регрессию или навесить решаюший лес. Дешево и сердито.

Возможно, этот подход не является значительно подкрепленным теорией, нарушает различные предпосылки, например, условия Гаусса-Маркова, особенно пункт про некоррелированность ошибок, однако на практике нередко выручает и достаточно активно используется в соревнованиях по машинному обучению.

Помимо стандартных признаков вроде лагов целевой переменной, много информации содержат в себе дата и время. Про извлечение признаков из них уже здорово описано в одной из предыдущих статей курса.

Добавлю только про еще один вариант кодирования категориальных признаков – кодирование средним. Если не хочется раздувать датасет множеством дамми-переменных, которые могут привести к потере информации о расстоянии, а в вещественном виде возникают противоречивые результаты а-ля «0 часов < 23 часа», то можно закодировать переменную чуть более интерпретируемыми значениями. Естественный вариант – закодировать средним значением целевой переменной. В нашем случае каждый день недели или час дня можно закодировать сооветствующим средним числом игроков, находившихся в этот день недели или час онлайн. При этом важно следить за тем, чтобы расчет среднего значения производился только в рамках тестового датасета (или в рамках текущего наблюдаемого фолда при кросс-валидации), иначе можно ненароком привнести в модель информацию о будущем.

def code_mean(data, cat_feature, real_feature):
    """
    Возвращает словарь, где ключами являются уникальные категории признака cat_feature, 
    а значениями - средние по real_feature
    """
    return dict(data.groupby(cat_feature)[real_feature].mean())

Создадим новый датафрейм и добавим в него час, день недели и выходной в качестве категориальных переменных. Для этого переводим имеющийся в датафрейме индекс в формат datetime, и извлекаем из него hour и weekday.

data = pd.DataFrame(dataset)
data.columns = ["y"]

data.index = data.index.to_datetime()
data["hour"] = data.index.hour
data["weekday"] = data.index.weekday
data['is_weekend'] = data.weekday.isin([5,6])*1
data.head()

Out:

Посмотрим на средние по дням недели

code_mean(data, 'weekday', "y")

Out:
{0: 38730.143229166664,
 1: 38632.828125,
 2: 38128.518229166664,
 3: 39519.035135135135,
 4: 41505.152777777781,
 5: 43717.708333333336,
 6: 43392.143603133161}

Помимо перечисленных преобразований для увеличения числа признаков используют и множество других метрик, например, максимальное/минимальное значение, наблюдавшееся в скользящем по ряду окне, медианы, число пиков, взвешенные дисперсии и многое другое. Автоматически этим занимается уже упоминавшаяся в курсе библиотека библиотека tsfresh.

Для удобства все преобразования можно записать в одну функцию, которая сразу же будет возвращать разбитые на трейн и тест датасеты и целевые переменные.

Функция для создания переменных

def prepareData(data, lag_start=5, lag_end=20, test_size=0.15):

    data = pd.DataFrame(data.copy())
    data.columns = ["y"]

    # считаем индекс в датафрейме, после которого начинается тестовыый отрезок
    test_index = int(len(data)*(1-test_size))

    # добавляем лаги исходного ряда в качестве признаков
    for i in range(lag_start, lag_end):
        data["lag_{}".format(i)] = data.y.shift(i)

    data.index = data.index.to_datetime()
    data["hour"] = data.index.hour
    data["weekday"] = data.index.weekday
    data['is_weekend'] = data.weekday.isin([5,6])*1

    # считаем средние только по тренировочной части, чтобы избежать лика
    data['weekday_average'] = map(code_mean(data[:test_index], 'weekday', "y").get, data.weekday)
    data["hour_average"] = map(code_mean(data[:test_index], 'hour', "y").get, data.hour)

    # выкидываем закодированные средними признаки 
    data.drop(["hour", "weekday"], axis=1, inplace=True)

    data = data.dropna()
    data = data.reset_index(drop=True)

    # разбиваем весь датасет на тренировочную и тестовую выборку
    X_train = data.loc[:test_index].drop(["y"], axis=1)
    y_train = data.loc[:test_index]["y"]
    X_test = data.loc[test_index:].drop(["y"], axis=1)
    y_test = data.loc[test_index:]["y"]

    return X_train, X_test, y_train, y_test

Линейная регрессия vs XGBoost

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

Построение линейной регрессии

from sklearn.linear_model import LinearRegression

X_train, X_test, y_train, y_test = prepareData(dataset.Users, test_size=0.3, lag_start=12, lag_end=48)
lr = LinearRegression()
lr.fit(X_train, y_train)
prediction = lr.predict(X_test)
plt.figure(figsize=(15, 7))
plt.plot(prediction, "r", label="prediction")
plt.plot(y_test.values, label="actual")
plt.legend(loc="best")
plt.title("Linear regressionn Mean absolute error {} users".format(round(mean_absolute_error(prediction, y_test))))
plt.grid(True);

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

Также можно провести оценку модели на кросс-валидации, тому же принципу, что был использован ранее. Для этого воспользуемся функцией (с небольшими модификациями), предложенной в посте Pythonic Cross Validation on Time Series

Код для кросс-валидации

def performTimeSeriesCV(X_train, y_train, number_folds, model, metrics):
    print('Size train set: {}'.format(X_train.shape))

    k = int(np.floor(float(X_train.shape[0]) / number_folds))
    print('Size of each fold: {}'.format(k))

    errors = np.zeros(number_folds-1)

    # loop from the first 2 folds to the total number of folds    
    for i in range(2, number_folds + 1):
        print('')
        split = float(i-1)/i
        print('Splitting the first ' + str(i) + ' chunks at ' + str(i-1) + '/' + str(i) )

        X = X_train[:(k*i)]
        y = y_train[:(k*i)]
        print('Size of train + test: {}'.format(X.shape)) # the size of the dataframe is going to be k*i

        index = int(np.floor(X.shape[0] * split))

        # folds used to train the model        
        X_trainFolds = X[:index]        
        y_trainFolds = y[:index]

        # fold used to test the model
        X_testFold = X[(index + 1):]
        y_testFold = y[(index + 1):]

        model.fit(X_trainFolds, y_trainFolds)
        errors[i-2] = metrics(model.predict(X_testFold), y_testFold)

    # the function returns the mean of the errors on the n-1 folds    
    return errors.mean()

%%time
performTimeSeriesCV(X_train, y_train, 5, lr, mean_absolute_error)

Size train set: (1838, 39)
Size of each fold: 367

Splitting the first 2 chunks at 1/2
Size of train + test: (734, 39)

Splitting the first 3 chunks at 2/3
Size of train + test: (1101, 39)

Splitting the first 4 chunks at 3/4
Size of train + test: (1468, 39)

Splitting the first 5 chunks at 4/5
Size of train + test: (1835, 39)
CPU times: user 59.5 ms, sys: 7.02 ms, total: 66.5 ms
Wall time: 18.9 ms

Out: 4613.17893150896

На 5 фолдах получили среднюю абсолютную ошибку в 4.6 K пользователей, достаточно близко к оценке качества, полученной на тестовом датасете.

Почему бы теперь не попробовать XGBoost…

Код для построения прогноза с XGBoost

import xgboost as xgb

def XGB_forecast(data, lag_start=5, lag_end=20, test_size=0.15, scale=1.96):

    # исходные данные
    X_train, X_test, y_train, y_test = prepareData(dataset.Users, lag_start, lag_end, test_size)
    dtrain = xgb.DMatrix(X_train, label=y_train)
    dtest = xgb.DMatrix(X_test)

    # задаём параметры
    params = {
        'objective': 'reg:linear',
        'booster':'gblinear'
    }
    trees = 1000

    # прогоняем на кросс-валидации с метрикой rmse
    cv = xgb.cv(params, dtrain, metrics = ('rmse'), verbose_eval=False, nfold=10, show_stdv=False, num_boost_round=trees)

    # обучаем xgboost с оптимальным числом деревьев, подобранным на кросс-валидации
    bst = xgb.train(params, dtrain, num_boost_round=cv['test-rmse-mean'].argmin())

    # можно построить кривые валидации
    #cv.plot(y=['test-mae-mean', 'train-mae-mean'])

    # запоминаем ошибку на кросс-валидации
    deviation = cv.loc[cv['test-rmse-mean'].argmin()]["test-rmse-mean"]

    # посмотрим, как модель вела себя на тренировочном отрезке ряда
    prediction_train = bst.predict(dtrain)
    plt.figure(figsize=(15, 5))
    plt.plot(prediction_train)
    plt.plot(y_train)
    plt.axis('tight')
    plt.grid(True)

    # и на тестовом
    prediction_test = bst.predict(dtest)
    lower = prediction_test-scale*deviation
    upper = prediction_test+scale*deviation

    Anomalies = np.array([np.NaN]*len(y_test))
    Anomalies[y_test<lower] = y_test[y_test<lower]

    plt.figure(figsize=(15, 5))
    plt.plot(prediction_test, label="prediction")
    plt.plot(lower, "r--", label="upper bond / lower bond")
    plt.plot(upper, "r--")
    plt.plot(list(y_test), label="y_test")
    plt.plot(Anomalies, "ro", markersize=10)
    plt.legend(loc="best")
    plt.axis('tight')
    plt.title("XGBoost Mean absolute error {} users".format(round(mean_absolute_error(prediction_test, y_test))))
    plt.grid(True)
    plt.legend()

XGB_forecast(dataset, test_size=0.2, lag_start=5, lag_end=30)

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

Заключение

Мы познакомились с разными методами и подходами к анализу и прогнозированию временных рядов. К сожалению, или к счастью, серебряной пули для решения такого рода задач пока не появилось. Методы, разработанные в 60-е годы прошлого века, (а некоторые и в начале 19-го), по-прежнему пользуются популярностью наравне с неразобранными в рамках данной статьи LSTM или RNN. Отчасти это связано с тем, что задача прогнозирования, как и любая другая задача, возникающая в процессе работы с данными — во многом творческая и уж точно исследовательская. Несмотря на обилие формальных метрик качества и способов оценки параметров, для каждого временного ряда часто приходится подбирать и пробовать что-то своё. Не последнюю роль играет и баланс между качеством и трудозатратами. Не раз уже упоминавшаяся здесь SARIMA-модель хотя и демонстрирует выдающиеся результаты при должной настройке, может потребовать не одного часа танцев с бубном манипуляций с рядом, в то время как простенькую линейную регрессию можно соорудить за 10 минут, получив при этом более-менее сравнимые результаты.

Домашнее задание №9

В демо-версии домашнего задания вы будете предсказывать просмотры wiki-страницы «Machine Learning». Веб-форма для ответов, там же найдете и решение.

Актуальные и обновляемые версии демо-заданий – на английском на сайте курса, вот первое задание. Также по подписке на Patreon («Bonus Assignments» tier) доступны расширенные домашние задания по каждой теме (только на англ.).

Полезные ресурсы

Open Machine Learning Course. Topic 9. Part 1. Time series analysis in Python
Видеозапись лекции по мотивам этой статьи
Open Machine Learning Course. Topic 9. Part 2. Predicting the future with Facebook Prophet
Онлайн-учебник курса по продвинутому статистическому прогнозированию университета Duke — разобраны всевозможные сглаживания, линейные модели и ARIMA модели
Статья Comparison of ARIMA and Random Forest time series models for prediction of avian influenza H5N1 outbreaks — одна из немногих, где активно защищается позиция случайного леса в задачах по прогнозированию временных рядов
Статья Time Series Analysis (TSA) in Python — Linear Models to GARCH семействе моделей ARIMA и их применении для моделирования финансовых показателей (Brian Christopher)

Материал статьи доступен в GitHub-репозитории курса
в виде тетрадок Jupyter.

Источник

Все курсы > Вводный курс > Занятие 20

Понятие временного ряда

Временной ряд (time series) — это данные, последовательно собранные в регулярные промежутки времени.

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

Основное отличие: перекрестные данные предполагают независимость наблюдений, во временных рядах будущее зависит от прошлого.

Работа с временными рядами предполагает два аспекта:

Анализ временного ряда (time series analysis), т.е. понимание его структуры и закономерностей; и
Моделирование и построение прогноза на будущее (time series forecasting)

Договоримся о терминах:

Во-первых, определим нотацию периодов. Временем t обозначим настоящее, t−1, t−2,… прошлое, t+1, t+2,… будущее.
Во-вторых, введем важное понятие временного лага (lag), т.е. запаздывания по сравнению с заданным периодом.

Датасеты

Мы будем использовать два популярных набора данных, а именно (1) ежемесячные данные о количестве пассажиров, перевезенных одной американской авиакомпанией с 1949 по 1960 годы и (2) ежедневные данные о родившихся в Калифорнии в 1959 году девочках.

Подгружать внешние данные в ноутбук Google Colab мы уже умеем. Можем переходить непосредственно к работе с кодом.

Откроем ноутбук к этому занятию⧉

Анализ временных рядов

# импортируем необходимые библиотеки

import pandas as pd

import numpy as np

import matplotlib.pyplot as plt

Импорт данных и работа в библиотеке Pandas

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

# импортируем файл с данными о пассажирах

passengers = pd.read_csv(«/content/passengers.csv»)

passengers.head()

Сделаем дату индексом.

# превратим дату в индекс и сделаем изменение постоянным

passengers.set_index(‘Month’, inplace = True)

passengers.head()

Питон воспринимает дату как число. Это не очень удобно, если мы хотим делать срезы и в целом изменять данные во времени. Дату можно преобразовать в специальный объект datetime.

# превратим дату (наш индекс) в объект datetime

passengers.index = pd.to_datetime(passengers.index)

# посмотрим на первые пять дат и на тип данных

passengers.index[:5]

DatetimeIndex([‘1949-01-01’, ‘1949-02-01’, ‘1949-03-01’, ‘1949-04-01’,

‘1949-05-01’],

dtype=’datetime64[ns]’, name=’Month’, freq=None)

Все эти операции также можно проделать в одну строчку.

passengers = pd.read_csv(‘/content/passengers.csv’, index_col = ‘Month’, parse_dates = True)

Теперь мы можем делать срезы за определенный период, например, с августа 1949 по март 1950 года.

passengers[‘1949-08’:‘1950-03’]

Обратите внимание, что в отличие от других перечней в Питоне (например, списков), во временном ряде мы получаем данные и по начальной, и по конечной дате среза (в частности, март 1950 года вошел в наш интервал).

Изменение шага временного ряда, сдвиг и скользящее среднее

Отдельно хотелось бы поговорить про возможность обобщения и изменения данных. Помимо прочего, мы можем изменить шаг (resample) нашего временного ряда, и посмотреть средние показатели перевозок, например, за год.

passengers.resample(rule = ‘AS’).mean().head()

Кроме того, мы можем сдвинуть (shift) наши данные на n периодов вперед или назад.

# произведем сдвиг на два периода (в данном случае месяца) вперед

passengers.shift(2, axis = 0).head()

Что логично, после сдвига первые два значения определяются как пропущенные (NaN или Not a number).

Мы также можем рассчитать скользящее среднее (moving average, rolling average) за n предыдущих периодов. Вначале посмотрим, что это такое.

Теперь давайте рассчитаем его для наших данных. Период, за который рассчитывается скользящее среднее, также называется окном (window).

# рассчитаем скользящее среднее для трех месяцев

passengers.rolling(window = 3).mean().head()

Опять же, так как в данном случае мы использовали три месяца для расчета скользящего среднего, этот показатель недоступен для первых двух значений. В целом, скользящее среднее сглаживает временные показатели (мы это увидим на графике в следующем разделе).

Построение графиков

Для того чтобы построить график временного ряда мы можем воспользоваться инструментами, которые уже содержатся в библиотеке Pandas. Например, простым методом .plot().

График можно усложнить.

# изменим размер графика, уберем легенду и добавим подписи

ax = passengers.plot(figsize = (12,6), legend = None)

ax.set(title = ‘Перевозки пассажиров с 1949 по 1960 год’, xlabel = ‘Месяцы’, ylabel = ‘Количество пассажиров’)

Визуализацию можно также построить с помощью библиотеки Matplotlib. Давайте выведем на одном графике перевозки пассажиров и скользящее среднее.

# зададим размер графика

plt.figure(figsize = (15,8))

# поочередно зададим кривые (перевозки и скользящее среднее) с подписями и цветом

plt.plot(passengers, label = ‘Перевозки пассажиров по месяцам’, color = ‘steelblue’)

plt.plot(passengers.rolling(window = 12).mean(), label = ‘Скользящее среднее за 12 месяцев’, color = ‘orange’)

# добавим легенду, ее положение на графике и размер шрифта

plt.legend(title = », loc = ‘upper left’, fontsize = 14)

# добавим подписи к осям и заголовки

plt.xlabel(‘Месяцы’, fontsize = 14)

plt.ylabel(‘Количество пассажиров’, fontsize = 14)

plt.title(‘Перевозки пассажиров с 1949 по 1960 год’, fontsize = 16)

# выведем обе кривые на одном графике

plt.show()

Как вы видите, скользящее среднее сильно сглаживает показатели. Также обратите внимание, что так как в данном случае мы взяли окно равное двенадцати месяцам, то первое значение скользящего среднего мы получили только за декабрь 1949 года (самое начало желтой кривой на графике).

В целом не стоит недооценивать важность визуальной оценки ряда на графике. Многие особенности можно выявить именно так.

Разложение временного ряда на компоненты

Выявление компонентов временного ряда (time series decomposition) предполагает его разложение на тренд, сезонность и случайные колебания. Дадим несколько неформальных определений.

Тренд — долгосрочное изменение уровня ряда
Сезонность предполагает циклические изменения уровня ряда с постоянным периодом
Случайные колебания — непрогнозируемое случайное изменение ряда

В Питоне в модуле statsmodels есть функция seasonal_decompose(). Воспользуемся ей для визуализации компонентов ряда.

Перед этим импортируем второй датасет для последующего сравнения.

# для наглядности импортируем второй датасет

# сразу превратим дату в индекс и преобразуем ее в объект datetime

births = pd.read_csv(‘/content/births.csv’, index_col = ‘Date’, parse_dates = True)

births.head(3)

Теперь давайте разложим наш временной ряд по авиаперевозкам на компоненты.

# импортируем функцию seasonal_decompose из statsmodels

from statsmodels.tsa.seasonal import seasonal_decompose

# задаем размер графика

from pylab import rcParams

rcParams[‘figure.figsize’] = 11, 9

# применяем функцию к данным о перевозках

decompose = seasonal_decompose(passengers)

decompose.plot()

plt.show()

Сделаем то же самое с данными о рождаемости.

decompose = seasonal_decompose(births)

decompose.plot()

plt.show()

Как мы видим, графики совершенно разные. В следующем разделе мы изучим это различие более подробно.

Стационарность

Стационарность (stationarity) временного ряда как раз означает, что такие компоненты как тренд и сезонность отсутствуют. Говоря более точно, среднее значение и дисперсия не меняются со смещением во времени.

Понимание того, стационарные ли у нас данные или нестационарные важно для последующего моделирования.

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

Для более точной оценки стационарности можно применить тест Дики-Фуллера (Dickey-Fuller test). О том, что такое статистический вывод мы с вами уже говорили.

В данном случае гипотезы звучат следующим образом.

Нулевая гипотеза предполагает, что процесс нестационарный
Альтернативная гипотеза соответственно говорит об обратном

Применим этот тест к обоим датасетам. Используем пороговое значение, равное 0,05 (5%).

# импортируем необходимую функцию

from statsmodels.tsa.stattools import adfuller

# передадим ей столбец с данными о перевозках и поместим результат в adf_test

adf_test = adfuller(passengers[‘#Passengers’])

# выведем p-value

print(‘p-value = ‘ + str(adf_test[1]))

p-value = 0.9918802434376409

Как мы видим, вероятность (p-value) для данных о перевозках существенно выше 0,05. Мы не можем отвергнуть нулевую гипотезу. Процесс нестанионарный. Проведем тест для второго набора данных.

# теперь посмотрим на данные о рождаемости

adf_test = adfuller(births[‘Births’])

# выведем p-value

print(‘p-value = ‘ +str(adf_test[1]))

p-value = 5.243412990149865e-05

Результат существенно меньше 5%. Временной ряд стационарен.

Надо сказать, что наша визуальная оценка полностью совпала с математическими вычислениями.

Автокорреляция

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

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

Допустим, у нас есть временной ряд и этот же ряд, взятый с лагом 1, 2 и 3.

Мы можем посчитать автокорреляцию ряда с лагом 1.

# для начала возьмем искусственные данные

data = np.array([16, 21, 15, 24, 18, 17, 20])

# для сдвига на одно значение достаточно взять этот ряд, начиная со второго элемента

lag_1 = data[1:]

# посчитаем корреляцию для лага 1 (у исходных данных мы убрали последний элемент)

# так как мы получим корреляционную матрицу, возьмем первую строку и второй столбец [0, 1]

np.round(np.corrcoef(data[:—1], lag_1)[0,1], 2)

Визуально это будет выглядеть следующим образом (вспомните график обычной корреляции).

# построим точечную диаграмму

plt.scatter(data[:—1], lag_1)

# добавим подписи

plt.xlabel(‘timeseries’, fontsize = 16)

plt.ylabel(‘lag 1’, fontsize = 16)

plt.title(‘Автокорреляция с лагом 1’, fontsize = 18)

Аналогично мы можем посчитать корреляцию для лагов 2 и 3 и на самом деле любого другого лага. Такие измерения автокорреляции удобно вычислить и изобразить с помощью графика автокорреляционной функции (autocorrelation function, ACF).

# импортируем автокорреляционную функцию (ACF)

from statsmodels.graphics.tsaplots import plot_acf

# применим функцию к нашему набору данных

plot_acf(data, alpha = None)

plt.show()

В частности, мы видим, что автокорреляция ряда с самим собой (первый столбец) равна 1, что логично. Второй столбец (то есть лаг 1) как раз примерно равен – 0,71.

Разные значения, полученные через np.corrcoef() и plot_acf(), объясняются небольшим различием в заложенных в этих функциях формулах.

Теперь построим график ACF для наших данных о перевозках.

# импортируем функцию для построения автокорреляционной функции (ACF)

from statsmodels.graphics.tsaplots import plot_acf

# применим ее к данным о пассажирах

plot_acf(passengers)

plt.show()

Автокорреляция позволяет выявлять тренд и сезонность, а также используется при подборе параметров моделей. В частности, мы видим, что лаг 12 сильнее коррелирует с исходным рядом, чем соседние лаги 10 и 11. То же самое можно сказать и про лаг 24. Такая автокорреляция позволяет предположить наличие (ежегодных) сезонных колебаний.

То, что корреляция постоянно положительная говорит о наличии тренда. Все это согласуется с тем, что мы узнали о данных, когда раскладывали их на компоненты.

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

Сравним полученный выше график с графиком автокорреляционной функции данных о рождаемости.

# построим аналогичный график для данных о рождаемости

plot_acf(births)

plt.show()

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

Моделирование и построение прогноза

На сегодняшнем занятии мы познакомимся с двумя типами моделей: экспоненциальное сглаживание и модель ARMA (и ее более продвинутые версии, ARIMA, SARIMA и SARIMAX).

Экспоненциальное сглаживание

Вновь обратимся к скользящему среднему (см. выше). В этой модели (1) всем предыдущим наблюдениям задавался одинаковый вес и (2) количество таких наблюдений было ограничено (мы называли это размером окна).

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

В модели экспоненциального сглаживания (exponential smoothing) или экспоненциального скользящего среднего мы как раз (1) берем все предыдущие значения и (2) задаем каждому из наблюдений определенный вес и (экспоненциально) уменьшаем этот вес по мере углубления в прошлое.

Приведем формулу.

$$ hat{y}_{t+1} = alpha cdot y_t + (1-alpha) cdot hat{y}_{t} $$

где ŷ_t+1 — это прогнозное значение, y_t — истинное значение в текущий период, ŷ_t — прогнозное значение в текущий период.

Как мы видим, прогнозное значение зависит как от истинного, так и от прогнозного значений. Важность этих значений определяется параметром альфа, который варьируется от 0 до 1. Чем альфа больше, тем больший вес у истинного наблюдения.

Формула рекурсивна, т.е. каждый раз мы умножаем (1 – α) на очередное прогнозное значение и так до конца временного ряда.

Практика

На Питоне эту модель не сложно прописать руками.

alpha = 0.2

# первое значение совпадает со значением временного ряда

exp_smoothing = [births[‘Births’][0]]

# в цикле for последовательно применяем формулу ко всем элементам ряда

for i in range(1, len(births[‘Births’])):

exp_smoothing.append(alpha * births[‘Births’][i] + (1 — alpha) * exp_smoothing[i — 1])

# выведем прогнозное значение для 366-го дня (1 января 1960 года)

exp_smoothing[—1]

Для временного ряда, состоящего из 365 наблюдений (весь 1959 год), мы получим 365 прогнозных значений (вплоть до 1 января 1960 года включительно).

len(births), len(exp_smoothing)

Теперь добавим эти данные в исходный датафрейм births.

# добавим кривую сглаживаия в качестве столбца в датафрейм

births[‘Exp_smoothing’] = exp_smoothing

births.tail(3)

Единственный нюанс, так как фактические значения описывают период с 1 января по 31 декабря 1959 года, а прогнозные со 2 января 1959 года по 1 января 1960, мы не можем просто их соединить. Второй столбец нужно сдвинуть на один день вперед.

# для этого импортируем класс timedelta

from datetime import timedelta

# возьмём последний индекс (31 декабря 1959 года)

last_date = births.iloc[[—1]].index

# # «прибавим» один день

last_date = last_date + timedelta(days = 1)

last_date

# добавим его в датафрейм

births = births.append(pd.DataFrame(index = last_date))

# значения за этот день останутся пустыми

births.tail()

Сдвинем второй столбец.

births[‘Exp_smoothing’] = births[‘Exp_smoothing’].shift(1)

И посмотрим на начало и конец датафрейма.

# как и должно быть первое прогнозное значение совпадает с предыдущим фактическим

births.head()

# и у нас есть прогноз на один день вперед

births.tail()

Проверим правильность расчетов, найдя последнее прогнозное значение по формуле выше.

$$ hat{y}_{366} = 0,2 times y_{355} + 0,8 times hat{y}_{355} = 0,2 times 50 + 0,8 times 45.756492 = 46.6051936 $$

Результаты совпадают. Теперь выведем фактические и прогнозные значения на графике.

# зададим размер

plt.figure(figsize = (15,8))

# выведем данные о рождаемости и кривую экспоненциального сглаживания

plt.plot(births[‘Births’], label = ‘Данные о рождаемости’, color = ‘steelblue’)

plt.plot(births[‘Exp_smoothing’], label = ‘Экспоненциальное сглаживание’, color = ‘orange’)

# добавим легенду, ее положение на графике и размер шрифта

plt.legend(title = », loc = ‘upper left’, fontsize = 14)

# добавим подписи к осям и заголовки

plt.ylabel(‘Количество родившихся’, fontsize = 14)

plt.xlabel(‘Месяцы’, fontsize = 14)

plt.title(‘Рождаемость в 1959 году. Прогноз на 1 января 1960 года’, fontsize = 16)

plt.show()

Очень советую в качестве упражнения попробовать построить график с несколькими значениями альфа от 0 до 1.

Модель экспоненциального сглаживания можно усложнить и тогда она будет улавливать тренд и сезонность. Кроме того, усложненные модели способны предсказывать более одного значения (обратите внимание, здесь мы смогли сделать прогноз лишь на один день вперёд).

Работа над ошибками. На видео в формуле экспоненциального сглаживания есть ошибка.

alpha = 0.2

# первое значение совпадает со значением временного ряда

exp_smoothing = [births[‘Births’][0]]

# в цикле for последовательно применяем формулу ко всем элементам ряда

for i in range(0, len(births[‘Births’])):

exp_smoothing.append(alpha * births[‘Births’][i] + (1 — alpha) * exp_smoothing[i — 1])

Должно быть.

alpha = 0.2

# первое значение совпадает со значением временного ряда

exp_smoothing = [births[‘Births’][0]]

# в цикле for последовательно применяем формулу ко всем элементам ряда

for i in range(1, len(births[‘Births’])):

exp_smoothing.append(alpha * births[‘Births’][i] + (1 — alpha) * exp_smoothing[i — 1])

Как вы видите расчет прогнозного значения в цикле следовало начинать с наблюдения
[births[‘Births’][1], а не
[births[‘Births’][0]. Из-за этого последующие результаты оказались некорректными. В лекции и ноутбуке код исправлен.

Теперь поговорим про модели семейства ARMA.

Модель ARMA

Модель ARMA состоит из двух компонентов.

компоненты модели ARMA: авторегрессия и скользящее среднее

Авторегрессия (autoregressive model, AR) — это регрессия ряда на собственные значения в прошлом. Другими словами, наши признаки в модели обычной регрессии мы заменяем значениями той же переменной, но за предыдущие периоды.

Когда мы прогнозируем значение в период t с помощью данных за предыдущий период (AR(1)), уравнение будет выглядеть следующим образом.

$$ y_t = c + varphi cdot y_{t-1} $$

где c — это константа, $varphi$ — вес модели, y_t–1 — значение в период t – 1.

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

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

Обратите внимание, что скользящее среднее временного ряда, которое мы рассмотрели выше, и модель скользящего среднего — это разные понятия.

Если использовать ошибку только предыдущего наблюдения, то уравнение будет выглядеть следующим образом.

$$ y_t = mu + varphi cdot varepsilon_{t-1} $$

где μ — это среднее значение временного ряда, φ — вес модели, ε_t–1 — ошибка в период t – 1.

Такую модель принято называть моделью скользящего среднего с параметром q = 1 или MA(1). Разумеется, параметр q может принимать и другие значения (MA(q)).

Модель ARMA с параметрами (или как еще говорят порядками, orders) p и q или ARMA(p, q) позволяет описать любой стационарный временной ряд.

ARMA предполагает, что в данных отсутствует тренд и сезонность (данные стационарны). Если данные нестационарны, нужно использовать более сложные версии этих моделей:

ARIMA, здесь добавляется компонент Integrated (I), который отвечает за удаление тренда (сам процесс называется дифференцированием); и
SARIMA, эта модель учитывает сезонность (Seasonality, S)
SARIMAX включает еще и внешние или экзогенные факторы (eXogenous factors, отсюда и буква X в названии), которые напрямую не учитываются моделью, но влияют на нее.

Параметров у модели SARIMAX больше. Их полная версия выглядит как SARIMAX(p, d, q) x (P, D, Q, s). В данном случае, помимо известных параметров p и q, у нас появляется параметр d, отвечающий за тренд, а также набор параметров (P, D, Q, s), отвечающих за сезонность.

Теперь давайте воспользуемся моделью SARIMAX для прогнозирования авиаперевозок.

Практика

В первую очередь нужно разбить данные на обучающую и тестовую выборки. Как мы помним, у нас есть данные с января 1949 года по декабрь 1960 года.

# обучающая выборка будет включать данные до декабря 1959 года включительно

train = passengers[:‘1959-12’]

# тестовая выборка начнется с января 1960 года (по сути, один год)

test = passengers[‘1960-01’:]

Посмотрим на разделение на графике.

plt.plot(train, color = «black»)

plt.plot(test, color = «red»)

# заголовок и подписи к осям

plt.title(‘Разделение данных о перевозках на обучающую и тестовую выборки’)

plt.ylabel(‘Количество пассажиров’)

plt.xlabel(‘Месяцы’)

# добавим сетку

plt.grid()

plt.show()

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

# принудительно отключим предупреждения системы

import warnings

warnings.simplefilter(action = ‘ignore’, category = Warning)

# обучим модель с соответствующими параметрами, SARIMAX(3, 0, 0)x(0, 1, 0, 12)

# импортируем класс модели

from statsmodels.tsa.statespace.sarimax import SARIMAX

# создадим объект этой модели

model = SARIMAX(train,

order = (3, 0, 0),

seasonal_order = (0, 1, 0, 12))

# применим метод fit

result = model.fit()

Теперь мы готовы делать прогноз. Вначале сделаем тестовый прогноз, соответствующий периоду тестовой выборки (1960 год), для того, чтобы оценить качество работы модели.

# тестовый прогнозный период начнется с конца обучающего периода

start = len(train)

# и закончится в конце тестового

end = len(train) + len(test) — 1

# применим метод predict

predictions = result.predict(start, end)

predictions

1960-01-01 422.703385

1960-02-01 404.947178

1960-03-01 466.293258

1960-04-01 454.781295

1960-05-01 476.848627

1960-06-01 527.162825

1960-07-01 601.449809

1960-08-01 610.821690

1960-09-01 513.229987

1960-10-01 455.692618

1960-11-01 409.200046

1960-12-01 450.754160

Freq: MS, Name: predicted_mean, dtype: float64

Построим соответствующий график.

# выведем три кривые (обучающая, тестовая выборка и тестовый прогноз)

plt.plot(train, color = «black»)

plt.plot(test, color = «red»)

plt.plot(predictions, color = «green»)

# заголовок и подписи к осям

plt.title(«Обучающая выборка, тестовая выборка и тестовый прогноз»)

plt.ylabel(‘Количество пассажиров’)

plt.xlabel(‘Месяцы’)

# добавим сетку

plt.grid()

plt.show()

В целом модель хорошо описывает временной ряд. Мы также можем использовать знакомые нам метрики среднеквадратической ошибки (MSE) и корня среднеквадратической ошибки (RMSE) для оценки качества.

# импортируем метрику

from sklearn.metrics import mean_squared_error

# рассчитаем MSE

print(mean_squared_error(test, predictions))

# и RMSE

print(np.sqrt(mean_squared_error(test, predictions)))

317.39564617023217

17.815601201481588

Теперь можно делать прогноз на будущее. Возьмём горизонт равный трем годам (1961, 1962 и 1963 год). Всего должно получиться 36 прогнозных значений.

# прогнозный период с конца имеющихся данных

start = len(passengers)

# и закончится 36 месяцев спустя

end = (len(passengers) — 1) + 3 * 12

# теперь построим прогноз на три года вперед

forecast = result.predict(start, end)

# посмотрим на весь 1963 год

forecast[—12:]

1963-01-01 518.603433

1963-02-01 497.908985

1963-03-01 556.406781

1963-04-01 542.133831

1963-05-01 561.524759

1963-06-01 609.244560

1963-07-01 681.016634

1963-08-01 687.950661

1963-09-01 587.995797

1963-10-01 528.167673

1963-11-01 479.454532

1963-12-01 518.856113

Freq: MS, Name: predicted_mean, dtype: float64

Посмотрим на прогнозные значения на графике.

# выведем две кривые (фактические данные и прогноз на будущее)

plt.plot(passengers, color = ‘black’)

plt.plot(forecast, color = ‘blue’)

# заголовок и подписи к осям

plt.title(‘Фактические данные и прогноз на будущее’)

plt.ylabel(‘Количество пассажиров’)

plt.xlabel(‘Месяцы’)

# добавим сетку

plt.grid()

plt.show()

Подведем итог

Занятие было насыщенным. Мы узнали, (1) чем временные ряды отличаются от структурных данных, (2) научились анализировать временные ряды с помощью Питона, а также (3) строить графики.

При прогнозировании временных рядов мы (4) использовали модель экспоненциального сглаживания и (5) познакомились с семейством моделей ARMA.

Вопросы для закрепления

Чем структурные/перекрестные данные отличаются от временных рядов?

Посмотреть правильный ответ

Ответ: во временных рядах значение наблюдаемой величины зависит от времени, в структурных данных — от других величин или переменных.

Что определяет параметр альфа в модели экспоненциального сглаживания?

Посмотреть правильный ответ

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

Из каких компонентов состоит модель ARMA?

Посмотреть правильный ответ

Ответ: модель ARMA состоит из компонента авторегрессии (AR), т.е. регрессии на собственные значения в прошлом, и модели скользящего среднего (MA), учитывающей отклонение от среднего значения или ошибку.

В рамках вводного курса осталась одна тема, основы нейронных сетей.

Ответы на вопросы

Вопрос. Вы не могли бы более детально рассказать, как подбираются параметры в модели SARIMAX. В частности, как были подобраны параметры (3, 0, 0)x(0, 1, 0, 12)?

Ответ. Параметры SARIMAX можно подобрать с помощью функции auto_arima. Выбор наиболее удачных параметров осуществляется на основе определенного критерия.

По умолчанию используется, так называемый, Akaike Information Criterion или AIC, который рассчитывается для каждой комбинации параметров. Чем ниже AIC модели, тем лучше.

Я привел код auto_arima в конце ноутбука⧉.

Вопрос. У меня уточнение: скажите, там, где мы считаем экспоненциальное сглаживание,

# в цикле for последовательно применяем формулу ко всем элементам ряда

for i in range(0, len(births[‘Births’])):

exp_smoothing.append(alpha * births[‘Births’][i] + (1 — alpha) * exp_smoothing[i — 1])

там, может быть, нужно в цикле
for i in range(0, len(births[‘Births’])): поставить начальное значение счетчика = 1:
for i in range(1, len(births[‘Births’]))?

Я проверял выкладки в Экселе — и, когда цикл стартует с нуля, получается ошибка. Возможно, это потому, что в начале списка получается задублированное значение первого элемента (мы сначала вручную вносим в список
[births[‘Births’][0]], а потом это же делаем на первой итерации, когда значение счетчика = 0)

Ответ. Да, действительно, первое прогнозное значение мы вносим вручную (оно совпадает с первым фактическим значением), поэтому должно быть

for i in range(1, len(births[‘Births’])):

Внес изменения в ноутбук и текст лекции.

Источник

Перевод

Ссылка на автора

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

Есть много разных показателей производительности на выбор. Может быть непонятно, какую меру использовать и как интерпретировать результаты.

В этом руководстве вы узнаете показатели производительности для оценки прогнозов временных рядов с помощью Python.

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

После завершения этого урока вы узнаете:

Основные показатели выполнения прогноза, включая остаточную ошибку прогноза и смещение прогноза.
Вычисления ошибок прогноза временного ряда, которые имеют те же единицы, что и ожидаемые результаты, такие как средняя абсолютная ошибка.
Широко используются вычисления ошибок, которые наказывают большие ошибки, такие как среднеквадратическая ошибка и среднеквадратичная ошибка.

Давайте начнем.

Ошибка прогноза (или остаточная ошибка прогноза)

ошибка прогноза рассчитывается как ожидаемое значение минус прогнозируемое значение.

Это называется остаточной ошибкой прогноза.

forecast_error = expected_value - predicted_value

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

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

expected = [0.0, 0.5, 0.0, 0.5, 0.0]
predictions = [0.2, 0.4, 0.1, 0.6, 0.2]
forecast_errors = [expected[i]-predictions[i] for i in range(len(expected))]
print('Forecast Errors: %s' % forecast_errors)

При выполнении примера вычисляется ошибка прогноза для каждого из 5 прогнозов. Список ошибок прогноза затем печатается.

Forecast Errors: [-0.2, 0.09999999999999998, -0.1, -0.09999999999999998, -0.2]

Единицы ошибки прогноза совпадают с единицами прогноза. Ошибка прогноза, равная нулю, означает отсутствие ошибки или совершенный навык для этого прогноза.

Средняя ошибка прогноза (или ошибка прогноза)

Средняя ошибка прогноза рассчитывается как среднее значение ошибки прогноза.

mean_forecast_error = mean(forecast_error)

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

Среднее значение ошибки прогноза, отличное от нуля, указывает на склонность модели к превышению прогноза (положительная ошибка) или занижению прогноза (отрицательная ошибка). Таким образом, средняя ошибка прогноза также называется прогноз смещения,

Ошибка прогноза может быть рассчитана непосредственно как среднее значение прогноза. В приведенном ниже примере показано, как среднее значение ошибок прогноза может быть рассчитано вручную.

expected = [0.0, 0.5, 0.0, 0.5, 0.0]
predictions = [0.2, 0.4, 0.1, 0.6, 0.2]
forecast_errors = [expected[i]-predictions[i] for i in range(len(expected))]
bias = sum(forecast_errors) * 1.0/len(expected)
print('Bias: %f' % bias)

При выполнении примера выводится средняя ошибка прогноза, также известная как смещение прогноза.

Bias: -0.100000

Единицы смещения прогноза совпадают с единицами прогнозов. Прогнозируемое смещение нуля или очень маленькое число около нуля показывает несмещенную модель.

Средняя абсолютная ошибка

средняя абсолютная ошибка или MAE, рассчитывается как среднее значение ошибок прогноза, где все значения прогноза вынуждены быть положительными.

Заставить ценности быть положительными называется сделать их абсолютными. Это обозначено абсолютной функциейабс ()или математически показано как два символа канала вокруг значения:| Значение |,

mean_absolute_error = mean( abs(forecast_error) )

кудаабс ()делает ценности позитивными,forecast_errorодна или последовательность ошибок прогноза, иимею в виду()рассчитывает среднее значение.

Мы можем использовать mean_absolute_error () функция из библиотеки scikit-learn для вычисления средней абсолютной ошибки для списка прогнозов. Пример ниже демонстрирует эту функцию.

from sklearn.metrics import mean_absolute_error
expected = [0.0, 0.5, 0.0, 0.5, 0.0]
predictions = [0.2, 0.4, 0.1, 0.6, 0.2]
mae = mean_absolute_error(expected, predictions)
print('MAE: %f' % mae)

При выполнении примера вычисляется и выводится средняя абсолютная ошибка для списка из 5 ожидаемых и прогнозируемых значений.

MAE: 0.140000

Эти значения ошибок приведены в исходных единицах прогнозируемых значений. Средняя абсолютная ошибка, равная нулю, означает отсутствие ошибки.

Средняя квадратическая ошибка

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

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

mean_squared_error = mean(forecast_error^2)

Мы можем использовать mean_squared_error () функция из scikit-learn для вычисления среднеквадратичной ошибки для списка прогнозов. Пример ниже демонстрирует эту функцию.

from sklearn.metrics import mean_squared_error
expected = [0.0, 0.5, 0.0, 0.5, 0.0]
predictions = [0.2, 0.4, 0.1, 0.6, 0.2]
mse = mean_squared_error(expected, predictions)
print('MSE: %f' % mse)

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

MSE: 0.022000

Значения ошибок приведены в квадратах от предсказанных значений. Среднеквадратичная ошибка, равная нулю, указывает на совершенное умение или на отсутствие ошибки.

Среднеквадратическая ошибка

Средняя квадратичная ошибка, описанная выше, выражается в квадратах единиц прогнозов.

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

rmse = sqrt(mean_squared_error)

Это можно рассчитать с помощьюSQRT ()математическая функция среднего квадрата ошибки, рассчитанная с использованиемmean_squared_error ()функция scikit-learn.

from sklearn.metrics import mean_squared_error
from math import sqrt
expected = [0.0, 0.5, 0.0, 0.5, 0.0]
predictions = [0.2, 0.4, 0.1, 0.6, 0.2]
mse = mean_squared_error(expected, predictions)
rmse = sqrt(mse)
print('RMSE: %f' % rmse)

При выполнении примера вычисляется среднеквадратичная ошибка.

RMSE: 0.148324

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

Дальнейшее чтение

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

Раздел 3.3 Измерение прогнозирующей точности, Практическое прогнозирование временных рядов с помощью R: практическое руководство,
Раздел 2.5 Оценка точности прогноза, Прогнозирование: принципы и практика
scikit-Learn Metrics API
Раздел 3.3.4. Метрики регрессии, scikit-learn API Guide

Резюме

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

В частности, вы узнали:

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

Есть ли у вас какие-либо вопросы о показателях эффективности прогнозирования временных рядов или об этом руководстве?
Задайте свои вопросы в комментариях ниже, и я сделаю все возможное, чтобы ответить.

Источник

1	2	3	4	5	6	7	8
37	7	18	4	87	7	18	3
38	8	19	4	88	8	19	4
39	9	20	4	89	9	20	5
40	10	21	4	90	10	21	6
41	1	12	5	91	1	12	7
42	2	13	5	92	2	13	1
43	3	14	5	93	3	14	2
44	4	15	5	94	4	15	3
45	5	16	5	95	5	16	4
46	6	17	5	96	6	17	5
47	7	18	5	97	7	18	6
48	8	19	5	98	8	19	7
49	9	20	5	99	9	20	1
50	10	21	5	0	10	21	2

1	2	3	4	5	6	7	8
37	7	18	4	87	7	18	3
38	8	19	4	88	8	19	4
39	9	20	4	89	9	20	5
40	10	21	4	90	10	21	6
41	1	12	5	91	1	12	7
42	2	13	5	92	2	13	1
43	3	14	5	93	3	14	2
44	4	15	5	94	4	15	3
45	5	16	5	95	5	16	4
46	6	17	5	96	6	17	5
47	7	18	5	97	7	18	6
48	8	19	5	98	8	19	7
49	9	20	5	99	9	20	1
50	10	21	5	0	10	21	2

Временные ряды в эконометрических исследованиях

Что такое ошибка прогноза?

Ошибка MAD — среднее абсолютное отклонение

Рассчитаем прогноз и оценим следующие модели:

Теперь рассчитаем прогноз с помощью «Логарифмического тренда с сезонностью».

Оценим эффективность использования в рамках года одной модели относительно другой.

Точных вам прогнозов!

Присоединяйтесь к нам!

Зарегистрируйтесь и скачайте решения

Статья полезная? Поделитесь с друзьями

Оценка ошибки прогнозирования временного ряда

MAPE – средняя абсолютная ошибка в процентах

MAE – средняя абсолютная ошибка

MSE – среднеквадратичная ошибка

RMSE – квадратный корень из среднеквадратичной ошибки

ME – средняя ошибка

SD – стандартное отклонение

Связь точности и ошибки прогнозирования

Введение

Движемся, сглаживаем и оцениваем

Rolling window estimations

Экспоненциальное сглаживание, модель Хольта-Винтерса

Простое экспоненциальное сглаживание

Двойное экспоненциальное сглаживание

Тройное экспоненциальное сглаживание a.k.a. Holt-Winters

Кросс-валидация на временных рядах, подбор параметров

Эконометрический подход

Стационарность, единичные корни

Избавляемся от нестационарности и строим SARIMA

Линейные и не очень модели на временных рядах

Линейная регрессия vs XGBoost

Заключение

Домашнее задание №9

Полезные ресурсы

Понятие временного ряда

Датасеты

Анализ временных рядов

Импорт данных и работа в библиотеке Pandas

Изменение шага временного ряда, сдвиг и скользящее среднее

Построение графиков

Разложение временного ряда на компоненты

Стационарность

Автокорреляция

Моделирование и построение прогноза

Экспоненциальное сглаживание

Модель ARMA

Подведем итог

Вопросы для закрепления

Ответы на вопросы

Ошибка прогноза (или остаточная ошибка прогноза)

Средняя ошибка прогноза (или ошибка прогноза)

Средняя абсолютная ошибка

Средняя квадратическая ошибка

Среднеквадратическая ошибка

Дальнейшее чтение

Резюме

Читайте также:

1	2	3	4	5	6	7	8
37	7	18	4	87	7	18	3
38	8	19	4	88	8	19	4
39	9	20	4	89	9	20	5
40	10	21	4	90	10	21	6
41	1	12	5	91	1	12	7
42	2	13	5	92	2	13	1
43	3	14	5	93	3	14	2
44	4	15	5	94	4	15	3
45	5	16	5	95	5	16	4
46	6	17	5	96	6	17	5
47	7	18	5	97	7	18	6
48	8	19	5	98	8	19	7
49	9	20	5	99	9	20	1
50	10	21	5	0	10	21	2