Анализ временных рядов — методические заметки #

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

Роли переменных #

Шкала столбца времени определяется один раз для всего столбца, а дата без времени суток привязывается к локальной полуночи. Если выбирать шкалу для каждого значения отдельно, в одном векторе оказались бы небольшие порядковые номера вперемешку с миллисекундами эпохи, а интервал между значениями с разных шкал — не шаг. Поэтому одной ячейки, которая распознаётся только как дата, достаточно, чтобы весь столбец считался столбцом дат, и любая похожая на число ячейка в нём считается ошибкой разбора, а не годом или индексом. Date.parse интерпретирует строку вида YYYY-MM-DD как полночь UTC, а YYYY-MM-DD HH:MM — как местное время, так что в столбце с обоими видами было бы две точки привязки. Модуль дописывает к дате без времени T00:00:00, и все значения становятся местными, то есть оказываются в той шкале, в которой их выводят ось времени на графиках и таблицы результатов. Если распознанных строк меньше трёх, интервалы считать не из чего, и расчёт отклоняется с уведомлением.

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

Сезонный период — календарный шаг, ближайший к медианному интервалу в логарифмической шкале, а при интервале короче примерно семнадцати часов — число наблюдений в сутки. Медианный интервал устойчив к немногим пропускам и дубликатам, которые встречаются в реальных файлах. Расстояние до шага оценивается как отношение, потому что сами шаги отличаются друг от друга в разы (неделя длиннее суток в 7 раз, месяц длиннее недели в 4.3 раза), а разброс внутри одного шага мал: месяцы от 28 до 31 дня отличаются от шага в 30.44 дня меньше чем на 8%, так что все они без всякой таблицы порогов считаются месячным шагом. Если наблюдений больше одного в сутки, циклом служат сами сутки, а период равен единице, делённой на интервал в сутках, с округлением (24 для часовых данных, 3 для восьмичасовых). Граница проходит на $\sqrt{0.5}$ суток, примерно на семнадцати часах: это среднее геометрическое между одним и двумя наблюдениями в сутки. Числовой индекс о шаге ничего не говорит, поэтому без ручного задания частота равна 1. Годовой шаг тоже даёт 1, поскольку более короткого цикла для моделирования нет.

Проверка. Месяцы в таблице авиапассажиров отстоят друг от друга на 28–31 день. Медиана приводится к месячному шагу, в строке частоты написано Месячный шаг — частота = 12 (годовой цикл), а findfrequency подтверждает период 12 по самому ряду.

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

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

Предварительная проверка (Проверить данные) #

Пропуск — это интервал, более чем в 1.5 раза шире медианного интервала группы, а дубликат — нулевой интервал. Календарный разброс мал (самый длинный месяц в 31/28 = 1.11 раза длиннее самого короткого), а одно пропущенное наблюдение удваивает интервал, так что множитель 1.5 разделяет эти случаи с запасом в обе стороны. Интервал меньше медианного, то есть два близких по времени измерения, пропуском не считается, а дубликатом — только если он в точности равен нулю. Строки для всего набора считает код на JavaScript, таблицу по группам — код на R, по одному и тому же правилу.

Проверка. Интервалы в таблице авиапассажиров составляют от 28 до 31 дня при медиане 31, так что порог пропуска равен 46.5 дня. Ни один месяц его не превышает: 0 пропусков, 0 дубликатов, регулярный шаг — Да.

Перед подсчётом разностей и поиском периода пропуски интерполируются, а поиск выбросов идёт по исходному ряду. ndiffs, nsdiffs и findfrequency не принимают NA, а удаление пустых значений сдвинуло бы все последующие на шаг раньше и разрушило бы сезонное выравнивание, на которое опираются критерии. Поэтому ряд с пропусками интерполируется функцией forecast::na.interp (сезонно, если позволяет частота) только для этих трёх расчётов, а если интерполяция не удалась, берётся исходный ряд. ndiffs применяет свой критерий KPSS по умолчанию с фиксированным α = 0.05, а не с уровнем из настроек приложения, и допускает до двух разностей. nsdiffs применяет эвристику по силе сезонности (одна сезонная разность, если сила сезонности по STL больше 0.64), и только при частоте больше 1 и хотя бы двух полных циклах в ряде. Обеим функциям, как и findfrequency, нужно восемь непропущенных значений. findfrequency удаляет линейный тренд, подгоняет к остатку авторегрессионный спектр и сообщает период в его пике, если пик достаточно высок, иначе 1.

Проверка. Для таблицы авиапассажиров приложение рекомендует d = 1 и D = 1 (сила сезонности ряда 0.94 при пороге 0.64) и обнаруживает период 12.

Для оценки λ по Герреро нужны десять положительных значений, а оценивается λ при сезонном периоде ряда. forecast::BoxCox.lambda(method = "guerrero") получает ряд как объект ts с установленной частотой, так что критерий Герреро на основе коэффициента вариации сравнивает разброс и уровень целых сезонов, то есть ту изменчивость, которую λ должно стабилизировать, а при частоте 1 — последовательных блоков из двух значений. Именно так forecast описывает работу с сезонным рядом, и так же поступает auto.arima(lambda = "auto"). Пропуски для этой оценки заполняются функцией forecast::na.interp, потому что их удаление сдвинуло бы последующие значения в чужой сезон. Если после заполнения появилось нулевое или отрицательное значение, оценка считается неудачной. При менее чем десяти значениях критерию не по чему подгоняться, и в ячейке стоит Не вычислено. Нулевое или отрицательное значение выводит ряд за область определения степенного преобразования, а сдвигать ряд модуль не станет, поэтому в ячейке стоит Неприменимо. Оценка, завершившаяся ошибкой, даёт Оценивание не удалось. Ряд не длиннее двух сезонных периодов получает Не оценено: на такой длине BoxCox.lambda возвращает λ = 1, ничего не подгоняя, и вывести эту единицу как оценку значило бы выдать тождественное преобразование за оценку по Герреро. Те же условия действуют на автоматический шаг Бокса — Кокса в наборе преобразований (он пропускает такой ряд с указанием причины) и на каждое обучающее окно, где набор применяется заново при оценке прогноза на отложенной выборке и при скользящей проверке. Окно, слишком короткое для оценки, подгоняется без преобразования, и строка набора считает такие окна, поскольку в них оценивается не та модель, что подогнана ко всему ряду. Заданное вручную λ применяется при любой длине.

Проверка. Для таблицы авиапассажиров предварительная проверка даёт λ = −0.295 (если передать те же значения простым вектором без периода, получилось бы 0.111), а подгонка с включённым флажком Заменить выбросы — 0.071, поскольку набор оценивает λ после замены семи летних пиков.

Выбросы ищет forecast::tsoutliers при заданной частоте, по исходному ряду, начиная с десяти непропущенных значений. Если ряд сезонный и длиннее двух циклов, робастная подгонка STL удаляет сезонную компоненту, когда её сила достигает 0.6, тренд подгоняется суперсглаживателем, а выбросом считается остаток ниже $Q_1-3\text{IQR}$ или выше $Q_3+3\text{IQR}$ остатков. Проходов два, и между ними отмеченные точки интерполируются. При менее чем десяти значениях ничего не отмечается. Поиск идёт по исходному ряду, поскольку его задача — заранее показать, что заменила бы очистка в наборе преобразований, а список позиций обрезается после двенадцатой, чтобы на ряде со множеством всплесков строка не разрасталась до абзаца.

Проверка. В таблице авиапассажиров отмечено семь летних пиков: июль 1949, 1950 и 1951 годов, август 1958 и 1959 годов, июль и август 1960 года. Летний размах растёт вместе с уровнем, поэтому при аддитивной сезонности самые высокие пики остаются крайними остатками. Это настоящие сезонные максимумы, а не ошибки, и решать, включать ли флажок Заменить выбросы, должен сам исследователь.

Рекомендация строится сначала по сезонной разности, затем по обычной, а пометка о коротком ряде ставится при любой из них. Рекомендованная сезонная разность означает сезонную модель, что бы ни показывало число обычных разностей. Если сезонной разности нет, обычная разность означает нестационарный ряд. Неопределённое число обычных разностей значит, что ни один критерий не удалось выполнить, а две нулевые разности — что ряд стационарен. Пометка о коротком ряде (меньше 30 непропущенных значений, то есть около двух циклов месячных данных с запасом) намеренно не зависит от рекомендации: короткий сезонный ряд — как раз тот случай, когда SARIMA и Холт — Винтерс принимают шум за сигнал, поэтому обе строки выводятся вместе.

Набор преобразований #

Замена выбросов идёт первой, по исходному ряду при заданной частоте, а λ оценивается уже по очищенному ряду. forecast::tsclean — это tsoutliers с последующей интерполяцией отмеченных точек и всех пропущенных значений. Точка считается заменённой, если её исходное и очищенное значения конечны и различаются больше чем на $\sqrt{\epsilon }$. Заполненные пропуски считаются отдельно, а сбой оставляет ряд без изменений и отражается в строке набора, а не теряется. Порядок важен, потому что всплеск искажает λ Герреро: оценка по очищенному ряду даёт степень для того ряда, который модель действительно получит.

Проверка. Для таблицы авиапассажиров набор заменяет 7 точек и затем оценивает λ = 0.071, тогда как предварительная проверка по исходному ряду давала −0.295.

λ передаётся функции подгонки, а не применяется заранее, поэтому обратное преобразование несмещено по среднему. Там, где функция λ не принимает, ряд преобразуется вручную. ets, auto.arima, Arima, naive, rwf, snaive и stlf принимают lambda = вместе с biasadj = TRUE, так что их подогнанные значения, прогнозы и интервалы возвращаются на исходную шкалу с поправкой на смещение от обратного преобразования, а величины, основанные на правдоподобии (информационные критерии, параметры, остатки), остаются на преобразованной шкале. Об этом и сообщает примечание о шкале на карточке. ma() и thetaf λ не принимают, а периодограмме и поиску разрывов нужен уже преобразованный ряд, поэтому там BoxCox применяется вручную и затем обращается. Заданное вручную λ — одно число на весь расчёт, одинаковое для всех рядов и групп. Пустое поле означает оценку для каждого ряда. При менее чем десяти непропущенных значениях, при нулевом или отрицательном значении в ряде (сдвигать ряд модуль не станет) или при неудачной оценке преобразование пропускается, и причина выводится на карточке, чтобы непреобразованный ряд не подгонялся без предупреждения.

В наборе нет дифференцирования. ETS, тэта-метод и Холт — Винтерс моделируют сам уровень: заранее продифференцированный ряд дал бы им ряд приращений, уровень которого — шум, а их прогнозы пришлось бы интегрировать обратно вручную. В ARIMA порядок d выполняет дифференцирование внутри её собственного правдоподобия, поэтому он задаётся на её карточке, а в наборе остаются только преобразования формы, которые корректно сочетаются с любым методом.

Методы анализа #

Разведочный анализ #

АКФ и ЧАКФ вычисляются по интерполированному ряду с его заданной частотой функциями forecast::Acf и Pacf, с na.pass и полосой 1.96/√n. forecast::Acf расширяет диапазон лагов по умолчанию до двух сезонных периодов, если это длиннее обычных 10·log₁₀(n), так что сезонный пик и его первое повторение остаются на коррелограмме, тогда как stats::acf оборвал бы месячный ряд на лаге 21. na.action = na.pass оставляет возможный остаточный NA на месте, а не смыкает пробел, потому что в укороченном векторе сдвинулся бы каждый последующий лаг. Лаг 0 на графике не показывается, как и в ggAcf: по построению он равен 1 и только растягивал бы ось. Полоса — это ±1.96/√n нулевой гипотезы белого шума, где n — число использованных наблюдений. Нижняя граница — четыре непропущенных значения, при меньшем числе оценка бессмысленна.

Проверка. Для таблицы авиапассажиров полоса равна ±0.163. АКФ выводится до лага 24, то есть на два периода, а не до лага 21, как дало бы 10·log₁₀(144). На лаге 1 она равна 0.948, на лаге 12 — 0.760. Значение ЧАКФ −0.540 на лаге 13 — характерный след сезонной разности.

Декомпозиция — это STL с периодическим сезонным окном и робастной подгонкой или stats::decompose. Ей нужны частота больше 1 и два полных цикла, а мультипликативной форме ещё и положительный ряд. s.window = "periodic" фиксирует форму сезонности во всех циклах, благодаря чему сила сезонности отражает заданную сезонность и ничего больше, а robust = TRUE уменьшает вес выпадающего месяца, так что всплеск попадает в остаток, а не в тренд. Классические формы — стандартные функции R, где тренд — центрированное скользящее среднее длиной в один период. При менее чем двух циклах оценивать сезонную компоненту не по чему (stl в этом случае просто отказывается работать), а у нулевого или отрицательного значения нет мультипликативной доли, поэтому такие расчёты пропускаются с заранее предусмотренным примечанием, а не завершаются ошибкой внутри R. Ряд с пропусками раскладывается после заполнения na.interp, которое само использует STL, когда позволяет сезонность, и не даёт спектральных артефактов простой линейной интерполяции. Если это заполнение не удалось, запасной вариант — линейная интерполяция с продлением крайних значений, поскольку stl не принимает NA.

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

Проверка. В таблице авиапассажиров отклонения остатка STL приходятся на летние месяцы последних лет, размах которых перерастает фиксированную периодическую сезонность: отмечено 18 точек, в том числе июнь — август 1956–1960 годов.

Сила тренда и сила сезонности — это $1-\mathrm{Var}(R)/\mathrm{Var}(T+R)$ и $1-\mathrm{Var}(R)/\mathrm{Var}(S+R)$ по Хайндману, ограниченные снизу нулём, а при мультипликативной декомпозиции — по логарифмам компонент. Определения взяты из функции stl_features пакета tsfeatures, так что числа сопоставимы с литературой, которая ими пользуется. Ограничение нулём срабатывает, когда дисперсия остатка больше дисперсии суммы компоненты и остатка и формула иначе дала бы отрицательную силу. Мультипликативная декомпозиция аддитивна в логарифмах, а её остаток колеблется около 1 при любом уровне шума, так что аддитивная формула, применённая к исходным мультипликативным частям, даёт число, прижатое к 1 независимо от данных. Логарифмирование трёх компонент возвращает мере смысл.

Проверка. Для таблицы авиапассажиров STL даёт силу тренда 0.942 и силу сезонности 0.618. nsdiffs в предварительной проверке находит для того же ряда силу сезонности 0.94: у этой функции своя STL с другим сезонным окном, так что совпадения двух чисел ожидать не следует.

Критерии стационарности — это tseries::adf.test, kpss.test с выбранной нулевой гипотезой и pp.test в форме Z(t_alpha), по интерполированному ряду, начиная с восьми непропущенных значений. p за пределами таблицы выводится как граница. adf.test подгоняет дрейф и тренд и выбирает порядок лагов как ⌊(n − 1)^(1/3)⌋, pp.test проверяет ту же альтернативу с поправкой Ньюи — Уэста и берёт лаг усечения Бартлетта по длине ряда, а форма Z(t_alpha) — статистика t-отношения, единственный вариант PP в шкале ADF, отсюда и обозначение PP (Zτ) в таблице. kpss.test получает null = "Level" или "Trend" из списка. Все три функции берут p-значения из таблиц интерполяцией и, если статистика выходит за пределы таблицы, возвращают её крайнее значение, сообщая об усечении только предупреждением. Поэтому модуль перехватывает это предупреждение, и в ячейке выводится < 0.01 или > 0.10, а не крайнее значение таблицы, как если бы оно было точным. Ячейка с границей сохраняет свой класс значимости и больше не перерисовывается как число. Критерии единичного корня не принимают na.action, так что ряд должен быть полным: берётся интерполированная копия, а если сама интерполяция не удалась, критерии пропускаются, а не выполняются на укороченном векторе, сдвинутые индексы которого исказили бы любую статистику. Ниже восьми значений порядок лагов adf.test равен 1, а разностей остаются считанные единицы, поэтому эта нижняя граница общая для всех критериев.

Проверка. Для таблицы авиапассажиров ADF равен −7.32 при лаге 5, PP равен −5.05 при лаге 4, оба < 0.01, а KPSS в форме для уровня равен 2.74 при лаге 4, тоже < 0.01: все ячейки — границы, поскольку статистики лежат далеко за пределами таблиц.

KPSS оценивается на уровне для проверки допущений, критерии единичного корня — на основном уровне, а заключение переворачивает утверждение самой нулевой гипотезы. Нулевая гипотеза KPSS — «ряд стационарен», то есть та же гипотеза «с данными всё в порядке», что у любой проверки допущений, поэтому KPSS использует settings.assumptionAlpha, как критерии нормальности и сферичности. У ADF и PP нулевая гипотеза противоположная, поэтому они сохраняют основной уровень значимости, и так пара остаётся согласованной, когда пользователь разводит два уровня. Отвержение гипотезы стационарности — вывод о нестационарности, а отвержение гипотезы единичного корня — вывод о стационарности, поэтому столбец заключения вычисляется по коду нулевой гипотезы и собственному α строки, и в каждом заключении выводится то α, на котором оценена эта строка.

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

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

Критерий Льюнга — Бокса для исходного ряда — это Box.test с max(10, 2m) лагами, но не больше n/5. По правилу Хайндмана для сезонного ряда берётся 2m лагов, чтобы сезонное повторение попало в окно, а для остальных — 10. Ограничение пятой частью ряда сохраняет точность χ²-аппроксимации статистики на коротком ряде. Нулевая гипотеза — отсутствие автокорреляции, поэтому критерий использует уровень для проверки допущений.

Проверка. Для таблицы авиапассажиров лаг равен min(24, 28) = 24, Q(24) = 1606, p < 0.001.

Поиск сезонных периодов — это пошаговый g-критерий Фишера по периодограмме без тейпера и без дополнения нулями, с границей Бонферрони для p, подавлением ±2 бинов на каждом шаге, не более чем десятью шагами и без повторной поправки. spec.pgram запускается с taper = 0 и fast = FALSE, потому что тейпер и дополнение нулями нарушают нулевую гипотезу о независимых одинаково распределённых экспоненциальных ординатах, на которую опирается g Фишера: тейпер снижает число степеней свободы ординаты ниже 2, а дополнение нулями увеличивает число ординат и делает соседние коррелированными. demean и detrend включены, чтобы масса на нулевой частоте не выдавала себя за цикл. p оценивается границей $N{(1-g)}^{N-1}$ вместо точной знакопеременной суммы Фишера: та численно неустойчива на длинных рядах, возможных в модуле (N ≈ 2500), а граница если и ошибается, то в консервативную сторону. На каждом шаге проверяется наибольшая оставшаяся ордината, после чего её главный лепесток (по два бина с каждой стороны) опускается до медианы периодограммы: цикл, не совпадающий точно с частотой Фурье, растекается по соседним бинам, и после удаления одного бина остались бы плечи, которые следующий шаг нашёл бы снова как слегка смещённый пик. Процедура останавливается через десять шагов или когда p очередного шага достигает 0.5, поскольку оставшаяся периодограмма тогда похожа на белый шум. Сохраняются пики с p не выше уровня значимости и периодом от 2 наблюдений до половины длины ряда. g более поздней строки измеряется на рабочей периодограмме с подавленными пиками и потому больше доли той же строки в нетронутой сумме, которую выводит столбец Доля дисперсии. Граница уже контролирует групповую вероятность ошибки по всем ординатам, так что глобальная поправка p-значений по K пикам скорректировала бы одно и то же семейство дважды.

Проверка. Для таблицы авиапассажиров проверено 10 пиков, и при α = 0.05 значимы 7: 12 месяцев с g = 0.502 и p = 2 × 10⁻²⁰, затем 6, 4, 3 и 2.4 месяца и два слабых длинных периода, 8 и 9.6 месяца, с p = 0.029 и 0.044. g строки 6 месяцев на периодограмме с подавленными пиками равно 0.404, а её доля в нетронутой периодограмме — 15.3%.

Строка считается гармоникой, если её интервал бинов Фурье пересекается с интервалом более сильной и более длинной строки, делённым на целое число от двух и больше. Несинусоидальный цикл с периодом P даёт настоящие пики на P/2, P/3, …, и критерий справедливо считает их значимыми. Допуском для отношения служит спектральное разрешение, которое и так выводится как границы периода, поэтому проверяется пересечение интервалов, а не округление P/p. Гармоника приписывается самому сильному подходящему основному циклу, поскольку кандидаты перебираются в порядке убывания g. Строка остаётся на своём месте в таблице (это настоящий пик), но приглушается и снабжается пометкой, а не удаляется.

Проверка. Для таблицы авиапассажиров 6, 4, 3 и 2.4 месяца помечены как ×2, ×3, ×4 и ×5 от 12. Период 8 месяцев не помечен: 12/8 округляется до 2, но интервал 6-месячной строки [5.9, 6.1] не доходит до интервала 8-месячной [7.8, 8.2].

Сглаживание #

Скользящее среднее — это forecast::ma(order = k, centre = TRUE) с окном, изначально равным сезонному периоду и не меньше 2, а при преобразовании Бокса — Кокса ряд преобразуется вручную и возвращается на шкалу медианы. При чётном порядке ma применяет взвешенную форму 2×m, так что окно 12 — это классическая центрированная оценка сезонного тренда, а каждое значение кривой охватывает m + 1 точку. Начальная ширина окна равна действующей частоте, потому что сезон убирает окно шириной в один сезон. Если периода нет, берётся 7 (неделя суточных данных). Окно 1 — это сам ряд, а 0 не определено, отсюда и нижняя граница. ma не принимает lambda, поэтому ряд преобразуется BoxCox до фильтра и обращается InvBoxCox после него. У фильтра нет дисперсии остатков для формулы поправки на смещение, поэтому обращённая кривая — медиана, а не среднее, о чём и сообщает карточка. Если непропущенных значений меньше ширины окна, ни одного полного окна нет, и подгонка пропускается.

Проверка. Для таблицы авиапассажиров окно 12 (начальное значение при частоте 12) сглаживает 132 точки из 144, теряя по шесть на каждом конце.

SES, Холт и Холт — Винтерс — это forecast::ses, holt и hw с biasadj = TRUE по интерполированному ряду, начиная с 4, 5 и 2m непропущенных значений. Это модели ETS: ANN, AAN или AAdN, AAA / MAM с их затухающими формами. Функции ets нужно на несколько наблюдений больше, чем параметров (четыре для уровня, пять для уровня и тренда, два полных цикла для сезонности), иначе она отказывается подгонять модель или возвращает вырожденную. Эти методы не могут перешагнуть через пропуск: они либо отказываются работать, либо сводят каждый параметр сглаживания к нижней границе. Поэтому ряд с пропусками заполняется na.interp до подгонки, а заполнение учитывается в строке набора преобразований. Скользящее среднее прямо показывает окна с NA, и его ряд остаётся без изменений. hw сразу отвергает сезонный период больше 24 и мультипликативную сезонность на ряде после преобразования Бокса — Кокса, а мультипликативные индексы не определены на ряде, касающемся нуля, поэтому каждый из этих случаев пропускается с заранее предусмотренным примечанием ещё до обращения к R. Выводимая сигма равна $\sqrt{{\sigma }^2}$ подогнанной модели. Таблица параметров — это вектор model$par: первые четыре его значения содержат веса сглаживания, а всё, что после них, ets хранит как начальные состояния.

Проверка. Для таблицы авиапассажиров SES оценивает α = 0.9999, то есть по сути даёт наивный прогноз (RMSE 33.6), а мультипликативный Холт — Винтерс ETS(M,A,M) даёт α = 0.315, β = 0.007, γ = 0.598, σ = 0.041, RMSE 10.6 и MAPE 2.91%.

Диагностика остатков — это критерий Льюнга — Бокса из checkresiduals (2m или 10 лагов, но не больше n/5, с fitdf, равным числу степеней свободы модели), АКФ остатков с полосой на уровне доверия и критерий Шапиро — Уилка не более чем на 5000 остатках, отобранных с зерном. Box.test вычитает fitdf из числа лагов, чтобы получить степени свободы, и при df ≤ 0 не выдаёт ошибку, а сообщает это значение. Поэтому если число лагов не превышает forecast::modeldf (число параметров, оценённых подгонкой, у ARIMA — вместе с константой), критерий пропускается с примечанием, а не выводит бессмысленное df. NA в остатках доходят до критериев как есть, потому что смыкание пробелов сдвинуло бы оценку на каждом лаге. shapiro.test принимает не больше 5000 наблюдений, поэтому более длинный вектор остатков оценивается по выборке из 5000, взятой с settings.randomSeed, и об этой выборке сообщается рядом с выводом. Нижняя граница — три остатка. Нулевые гипотезы обоих критериев утверждают, что модель подходит (остатки — белый шум, остатки нормальны), поэтому оба используют уровень для проверки допущений.

Проверка. Для таблицы авиапассажиров остатки мультипликативного метода Холта — Винтерса дают Q(24) = 40.2 при 24 − 16 = 8 степенях свободы, p = 3 × 10⁻⁶, и W = 0.990, p = 0.43. Остатки автоматически подобранной ARIMA дают Q(24) = 41.8 при 23 степенях свободы, p = 0.010, и W = 0.978, p = 0.021.

Моделирование ARIMA и SARIMA #

Автоматический режим — это forecast::auto.arima с выбранным критерием, где stepwise и approximation по умолчанию включены, а при полном переборе выключены. Нужно не меньше десяти непропущенных значений. Пошаговый поиск с приближённым критерием — вариант Хайндмана и Хандакара по умолчанию, и на длинном ряде он занимает секунды. Точная оценка каждого порядка нередко находит модель с меньшим значением критерия, но на длинном или высокочастотном ряде занимает минуты, поэтому она включается отдельно, а не по умолчанию. Порядок считывается функцией arimaorder, которая сообщает P, D и Q только для сезонной модели, а константа — по именам коэффициентов (drift на дифференцированном ряде, intercept на недифференцированном), чтобы в заголовке, как и в выводе summary.Arima, стояло с дрейфом. При менее чем десяти значениях у auto.arima слишком мало разностей для выбора порядка, а сетка соседних порядков для ручной подгонки теряет смысл.

Проверка. Для таблицы авиапассажиров поиск по умолчанию возвращает ARIMA(0,1,1)(0,1,0)[12] с ma1 = −0.318, AICc 1020.7, RMSE 11.2 и MASE 0.258.

Ручной режим — это forecast::Arima с include.constant, сезонными порядками, обнуляемыми при частоте 1, и сеткой соседних порядков шириной 3 по каждой оси при закреплённых d и D. Arima без ошибки оценивает сезонный член при периоде 1: он смешивается с несезонным, а arimaorder сообщает о нём так, будто его нет. Поэтому при частоте 1 сезонная строка отключается на панели и обнуляется в расчёте. include.constant добавляет среднее к недифференцированному ряду и дрейф к однократно продифференцированному, как auto.arima, а при более высоком дифференцировании не добавляет ничего. Сетка меняет p, q, P и Q на шаг в обе стороны от выбранного порядка (при нулевом порядке — {0, 1, 2}, чтобы сравнение оставалось симметричным и показывало, помогает ли переход к 2) и закрепляет d и D, поскольку информационные критерии при разном дифференцировании несопоставимы. Это 9 кандидатов при несезонной частоте и 81 при сезонной, и каждый подгоняется полностью, поэтому сетка строится только в ручном режиме и показывает собственный индикатор хода. Границы (p, q ≤ 10; d ≤ 3; P, Q ≤ 5; D ≤ 2) — это min/max элементов управления. Пустой порядок или порядок вне диапазона останавливает запуск с уведомлением ещё до подгонки, а замена на 0 в сборщике параметров нужна только для нетронутой панели.

Проверка. Для таблицы авиапассажиров классическая модель авиаперевозок ARIMA(0,1,1)(0,1,1)[12] даёт ma1 = −0.309 и sma1 = −0.107 при AICc 1021.2. Из 81 соседней модели подгоняются 79, и лучшая из них — ARIMA(1,1,2)(1,1,2)[12] с AICc 1013.2.

Интервалы коэффициентов — это оценка ± z × SE на глобальном уровне доверия, а p — по нормальному распределению. Точность на обучающей выборке считает forecast::accuracy по обратно преобразованным подогнанным значениям, а при периоде больше 24 карточка выводит предупреждение. Arima сообщает var.coef, диагональ которой даёт стандартные ошибки, а асимптотическое распределение оценок нормально, поэтому правильное опорное распределение — z, а не t. accuracy() оценивает fitted(), которые при biasadj = TRUE возвращаются на исходной шкале, так что при преобразовании Бокса — Кокса четыре строки точности — единственные числа карточки на исходной шкале, а критерии, σ² и остатки остаются на преобразованной. Об этом и сообщает примечание о шкале. Примерно с сезонного периода 24 auto.arima незаметно сводится к одной сезонной разности, а ручной сезонный член перестаёт сходиться. Это собственный предел Хайндмана для SARIMA. Поэтому здесь карточка предупреждает и называет альтернативу: динамическую гармоническую регрессию и методы на основе STL.

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

Семь методов прогноза — это naive, rwf(drift = TRUE), snaive, forecast(ets(…)), forecast(auto.arima(…)), обёртка тэта-метода и stlf(robust = TRUE, s.window = "periodic"), каждый с lambda и biasadj = TRUE, если он их принимает, и с заранее проверяемым условием применимости. Набор охватывает базовые прогнозы, два семейства автоматических моделей, тэта-метод, победивший в соревновании M3, и путь через STL для длинных или неудобных сезонов. Интервалы предсказания запрашиваются на уровнях 80 и 95. stl нужны два полных периода, snaive один, и обоим нужна частота больше 1. Шаг АКФ в thetaf не выполняется на любом пропуске, который не удалось заменить заполнением. Во всех этих случаях R выдаёт лишь невнятную ошибку, поэтому модуль сначала проверяет условие, и в строке стоит Неприменимо с причиной, а подгонка, не удавшаяся на самих данных, даёт Не выполнено. thetaf не принимает lambda, поэтому её вход преобразуется по Боксу — Коксу вручную, а среднее и подогнанные значения обращаются с biasadj, причём дисперсия прогноза восстанавливается по полуширине 95-процентного интервала. Так прогноз тэта-метода, как и остальные на том же графике, несмещён по среднему. Кроме того, thetaf раскладывает ряд мультипликативно всякий раз, когда срабатывает её собственный критерий сезонности, а на ряде, пересекающем ноль, такое разложение не определено. Поэтому этот критерий воспроизводится в модуле, и неположительный сезонный ряд корректируется аддитивно.

Отложенная выборка — это последние min(h, ⌊0.2 n⌋) наблюдений ряда, в котором не меньше 30 непропущенных значений и обучающий отрезок не короче 10. Она оценивается по конечным фактическим значениям, а набор преобразований заново применяется к обучающему отрезку. Пятая часть ряда — обычный предел для проверочного окна, а горизонт показывает, на сколько вперёд прогноз действительно потребуется, поэтому берётся меньшее из двух. При менее чем 30 наблюдениях любое окно, которое стоило бы оценивать, оставляет для подгонки слишком мало данных, так что блок пропускается, а строка валидации называет условие. Окно задаётся по позициям, а условие считает непропущенные значения, поэтому ошибки вычисляются по конечным фактическим значениям внутри окна. Набор преобразований применяется к обучающему отрезку заново, потому что замена выбросов или λ, оценённое с учётом окна, передали бы сведения об окне в ту самую подгонку, которую оно должно проверять.

Проверка. Для таблицы авиапассажиров отложенная выборка составляет min(12, 28) = 12 месяцев, то есть 1960 год. Авто-ARIMA получает на ней RMSE 23.9 и MASE 0.61.

Скользящая перекрёстная проверка размещает точки отсчёта на целый горизонт друг от друга, использует min(запрошенное число, ⌊(n − 10)/h⌋) из них и оценивает только те, до которых дошли все подогнанные методы. Непересекающиеся окна не дают оценить одно наблюдение дважды, поэтому объединённая ошибка — равновзвешенное среднее, а не среднее с завышенным весом последнего отрезка. Поэтому же каждое окно требует ещё одного горизонта наблюдений, а модуль использует столько окон, сколько позволяет ряд, и сообщает a из n, если их меньше запрошенного. Если подгонка метода не удалась в какой-то точке отсчёта, без особого правила эта точка осталась бы в общем наборе для соперников, но не для него. Поэтому точка оценивается только тогда, когда прогноз для неё дали все участники состязания, а метод, не дошедший ни до одной точки, выбывает из перекрёстной проверки и не опустошает общий набор. Максимум в 20 окон ограничивает затраты (20 × 7 подгонок на ряд) и при запуске без интерфейса, и в элементе управления.

Проверка. Для таблицы авиапассажиров оцениваются все пять окон с шагом в год: авто-ARIMA — RMSE 27.0 и MASE 0.71; ETS — 30.5 и 0.81; тэта-метод — 33.3 и 0.90; сезонный наивный — 42.5 и 1.33; STL + ETS — 60.8; дрейф — 80.8; наивный — 90.4 и 2.30.

MASE делится на MAE сезонного наивного прогноза на том отрезке, где обучалась модель, MAPE не определена, если хоть одно оцениваемое фактическое значение не положительно, а ранжирует весь ряд одна метрика — RMSE перекрёстной проверки, иначе RMSE отложенной выборки. Масштаб Хайндмана и Кёлер — ошибка наивного прогноза на лаге m на обучающих данных (на лаге 1 при несезонной частоте), так что MASE меньше 1 означает, что метод лучше повторения прошлого сезона на данных, которые он тоже видел. При нулевом знаменателе MASE не определена. MAPE делит ошибку на фактическое значение, поэтому значение, равное нулю или близкое к нему (в пределах 10⁻⁶ от наибольшего фактического), делает её неопределённой, а не огромной. Перекрёстная проверка и отложенная выборка оценивают разные данные, поэтому корректно лишь ранжирование всего ряда по одной метрике: RMSE перекрёстной проверки, если оценено хотя бы одно окно, иначе RMSE отложенной выборки. Метод без этой метрики идёт после всех оценённых, а метод, который так и не удалось подогнать, — последним, что бы ни содержала метрика. Заголовок ранга называет метрику, чтобы порядок не принимали за результат перекрёстной проверки, когда он получен иначе.

Аномалии внутри выборки — это остатки, отстоящие от своей медианы больше чем на три MAD, относительно исходного ряда, начиная с пяти конечных остатков. Остатки считаются по исходным наблюдениям, а не по преобразованному ряду, поэтому точка, которую заменил tsclean, тоже может быть отмечена: отметка означает, что наблюдение необычно для модели, что бы ни сделал с ним набор преобразований. Правило 3·MAD — то же, что на карточке разведочного анализа, так что оба маркера определяют выброс одинаково.

Будущие метки времени идут по календарным месяцам, если интервал с допуском равен целому числу месяцев, и отсчитываются от опорной даты, а не накапливаются. При любом другом шаге ось продолжается линейно. Месячный интервал составляет от 28 до 31 дня, и экстраполяция с постоянным числом миллисекунд уже за год уходит от границ месяцев. Поэтому интервал, равный целому числу месяцев по 30.44 дня с допуском 3% (но не меньше двух дней, что поглощает разброс длины месяцев при месячном шаге и всё же отвергает четырёхнедельный шаг в 28 дней), проходится по календарным месяцам. Каждый шаг отсчитывается от опорной даты и ограничивается длиной целевого месяца, потому что при накоплении шагов опорная дата в конце месяца переполняется (31 января плюс месяц — это 3 марта), и период просто пропускается. Интервалы предсказания фиксированы на 80% и 95%: эти два уровня все методы forecast возвращают по умолчанию, и только их график и таблица могут показать без отдельного элемента управления уровнем.

Спектральный анализ / точки разрыва #

Периодограмма — это spec.pgram с 10-процентным тейпером и удалением тренда на истинной сетке Фурье, для графика сглаженная модифицированным ядром Даниэля шириной порядка √n. Пять сильнейших локальных максимумов с периодом до n/2 проверяются g-критерием Фишера по второй периодограмме, без тейпера. fast = FALSE сохраняет истинную сетку Фурье: дополнение нулями до длины с большим числом делителей создаёт бины ниже 1/n, то есть периоды длиннее самого ряда. Тейпер — стандартное косинусное окно на 10% с каждого края — уменьшает просачивание сильных пиков в соседние частоты. У несглаженной периодограммы около двух степеней свободы на ординату и выглядит шумным гребнем, поэтому рядом с ней рисуется сглаживание модифицированным ядром Даниэля шириной max(3, 2⌊√n/4⌋ + 1) (рядом, а не вместо, поскольку g Фишера определён на несглаженных ординатах), и его полоса указывается под графиком. Пики — локальные максимумы, чтобы соседние бины одного пика не заполняли список. Один-два наблюдённых цикла ещё не доказывают период, поэтому ранжированный список обрывается на n/2, а на графике остаётся весь спектр. g каждого пика берётся по второй периодограмме, без тейпера, потому что тейпер снижает число степеней свободы ординаты ниже 2 и нарушает нулевую гипотезу о независимых одинаково распределённых экспоненциальных ординатах. Здесь g — доля пика во всём спектре с той же границей Бонферрони для p, что и на карточке разведочного анализа, так что первые строки двух карточек совпадают, а следующие расходятся: на карточке разведочного анализа каждый более сильный пик подавлялся перед проверкой следующего. После удаления тренда постоянный ряд даёт тождественно нулевой спектр, в котором все бины равны, поэтому такой случай пропускается, а не ранжируется. Восемь наблюдённых значений — нижняя граница для спектра хоть с каким-то разрешением.

Проверка. Для таблицы авиапассажиров 72 ординаты сглаживаются с шириной 7 и полосой 0.0125 цикла на наблюдение. Пики приходятся на 12 месяцев (g = 0.502, p = 2 × 10⁻²⁰), 6 (g = 0.153, p = 0.0005), 4, 3 и 48 (p = 1), тогда как на карточке разведочного анализа строка 6 месяцев даёт g = 0.404 по периодограмме с подавленными пиками.

Точки разрыва — это strucchange::breakpoints(y ~ 1) или (y ~ t) с h, равным доле минимального сегмента, числом разрывов по BIC и интервалами confint, ограниченными пределами ряда. Средние сегментов берутся по непреобразованному ряду, и нужно не меньше двадцати наблюдённых значений. y ~ 1 — чистая модель сдвига среднего, которая представляет любой тренд лестницей из ступенек уровня, а y ~ t — модель кусочно-линейного тренда. breakpoints выбирает число разрывов по BIC в поиске динамическим программированием, а h ограничивает, сколько разрывов можно поставить: диапазон элемента управления от 5 до 45% оставляет место хотя бы для двух сегментов, поскольку поиск упирается в ограничение задолго до половины ряда. Получившаяся длина в наблюдениях выводится, так что видно, если значение было ограничено. confint завершается ошибкой, если разрыв не найден, отсюда отдельная ветка для случая без разрывов. Возвращаемый ею интервал может выходить за любой конец ряда и там обрезается. Среднее, взятое на шкале Бокса — Кокса и пропущенное через InvBoxCox, — медиана, а не среднее данных, которые рисует график, поэтому средние сегментов считаются по непреобразованному интерполированному ряду. При менее чем двадцати значениях сегменты при любом допустимом h слишком коротки, чтобы оценить уровень или прямую. strucchange зависит от пакетов zoo, sandwich и lattice, поэтому загружается, только если поиск включён, а сбой загрузки сводится к предупреждению в разделе точек разрыва, тогда как периодограмма всё равно строится.

Проверка. Для таблицы авиапассажиров модель сдвига среднего при h = 15% (минимальный сегмент — 21 наблюдение) ставит пять разрывов, после наблюдений 26, 50, 77, 98 и 123 (первый датирован в пределах [24, 28]), со средними сегментов 134, 188, 239, 318, 379 и 464. Модель «уровень + тренд» не ставит ни одного.

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

Импутация, пропуски и отброшенные строки #

Пустое значение ряда остаётся на месте и заполняется каждой карточкой по-своему, никогда не удаляясь вместе со строкой, а глобальная обработка пропусков действует ещё до модуля. ts() не использует столбец времени: время в создаваемом ею объекте задаётся позицией. Поэтому удаление неполной строки сдвигает все последующие наблюдения на период раньше и сбивает сезоны, а сохранённый NA оставляет сетку на месте, и каждая подгонка обрабатывает его по-своему. Arima и auto.arima проходят пропуск фильтром Калмана, методы прогноза с шагом по лагам (naive, rwf и snaive) через него перешагивают, stlf сам интерполирует внутри mstl, а подгонки и критерии, не принимающие NA (ses, holt, hw, критерии единичного корня и KPSS, Box.test, stl, spec.pgram, breakpoints), получают заполнение na.interp, которое следует за трендом и сезонностью ряда, и это заполнение учитывается на карточке. ets и thetaf на карточке прогнозирования тоже не принимают NA и получают такое же заполнение, тогда как остальные методы прогноза продолжают подгоняться прямо по ряду с пропусками. Заполнение выполняется заново на каждом обучающем отрезке (на всём ряде, на отрезке отложенной выборки и на отрезке каждого окна перекрёстной проверки), так что на отрезок не влияют значения после его конца, а точность оценивается только по наблюдённым значениям, и заполненная точка в оценку не входит. Модуль отбрасывает только строки, которым нет места на оси, то есть с нераспознанным временем или с пустой группой при группирующей переменной. Их число выводится в предварительной проверке и в предупреждении на каждой карточке. Глобальная обработка пропусков выполняется над набором данных раньше, чем его получает любой модуль, поэтому модуль не может её отключить и не делает вид, что может. Построчное удаление убирает неполные строки ещё до разбора оси времени, и они появляются в диагностике шага как пропуски на оси, а не как пропущенные значения. Импутация одним значением заполняет каждую пустую ячейку столбца одной константой, самым плоским заполнением, какое только может получить сезонный ряд, и карточкам уже нечего сообщать. Руководство описывает оба варианта и рекомендует попарное удаление, при котором пустые ячейки доходят до карточек как NA и применяется собственное заполнение каждой карточки.

Проверка. Если в таблице авиапассажиров сделать пустыми три месяца (май 1953, ноябрь 1956 и февраль 1959 года), карточка прогнозирования проводит состязание всех семи методов, ETS и тэта-метод подгоняются по заполненному ряду (заполнено 3 значения), и первое место получает авто-ARIMA с RMSE перекрёстной проверки 27.4, опережая ETS с 30.4 и тэта-метод с 33.3. При построчном удалении предварительная проверка показывает 141 наблюдение, 3 пропуска на оси и 0 пропущенных значений.

Воспроизводимость #

Ссылки на литературу соответствуют расчёту, а не возможностям модуля: forecast и его статья — всегда, tseries — только если выбран критерий стационарности, strucchange — только если выполнялся поиск разрывов, а каждый метод — из той ветви, которая его вычислила. Список литературы со всеми методами, которые модуль умеет применять, был бы бесполезен для раздела «Методы», поэтому каждый расчёт ссылается на то, что сделал. Предварительная проверка ссылается на KPSS, потому что ndiffs принимает решение по нему, на Бокса и Кокса и на Герреро, потому что оценивает λ, на суперсглаживатель, потому что его подгоняет tsoutliers, а на STL и меру силы — только при частоте больше 1. Разведочный анализ при любых параметрах ссылается на АКФ и ЧАКФ, критерий Льюнга — Бокса, периодограмму, g Фишера и границу Бонферрони, на которых стоит поиск сезонных периодов. На критерий стационарности он ссылается, только если критерий был отмечен, а на метод декомпозиции и правило 3·MAD — только если декомпозиция выполнялась. Карточки с подгонкой ссылаются на своё семейство моделей (выбранный метод сглаживания, ARIMA и автоматический поиск в автоматическом режиме, семь методов прогноза), а если карточка выводит три информационных критерия и MASE, то и на них. На диагностику они ссылаются, только если она была включена, а на скользящую перекрёстную проверку и правило аномалий — только если они выполнялись. Набор преобразований ссылается на суперсглаживатель и STL только при включённом флажке Заменить выбросы, на Бокса и Кокса только при включённом флажке Преобразование Бокса — Кокса, а на Герреро только тогда, когда λ оценивалось, а не было задано вручную. Единственная случайная выборка — подвыборка для критерия Шапиро — Уилка, когда остатков больше 5000, — обоснована в разделе о диагностике остатков. Больше нигде модуль случайных чисел не использует.