DengAI: Прогнозирование распространения заболеваний - Прогнозирование STL / ARIMA / Box-Jenkins

Использование метода прогнозирования STL с моделью ARIMA, параметризованной с помощью метода Бокса-Дженкинса.

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

Для тех, кто не знаком с проблемой прогнозирования, этот конкурс посвящен прогнозированию лихорадки денге в двух городах, или, по словам самих DrivenData:

Ваша цель - предсказать метку total_cases для каждого (city, year, weekofyear) в наборе тестов. Есть два города, Сан-Хуан и Икитос, с данными испытаний для каждого города, охватывающими 5 и 3 года соответственно. Вы сделаете одну заявку, содержащую прогнозы для обоих городов.

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

В центре нашего подхода к прогнозированию находится вышеупомянутый так называемый метод STL. STL расшифровывается как S easonal - T метод разложения с использованием L OESS. Этот метод разбивает временной ряд на три компонента, а именно на его тренд, сезонность и остатки. LOESS означает сглаживание локально оцененной диаграммы рассеяния и извлекает гладкие оценки трех вышеупомянутых компонентов.

Структурные изменения

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

Здесь важно подчеркнуть более высокую важность новых данных по сравнению со старыми данными в задачах временных рядов. Обычно в проектах по науке о данных мы всегда заинтересованы в сборе большего количества данных, так как это в целом повысит надежность и точность нашей модели. Хотя, когда дело доходит до данных временных рядов, мы должны взвесить важность новых данных в зависимости от того, когда именно данные были отобраны. Если мы, например, заинтересованы в прогнозировании курса акций определенной компании, финансовые отчеты компании за последние пару лет, несомненно, более важны, чем результаты компании 100 лет назад. Почему это так? Потому что вполне возможно, что базовый процесс генерации данных со временем изменился. Это означало бы, что более старые данные не дают нам никакой информации о будущих данных и, следовательно, не только неверны, но также могут отрицательно сказаться на производительности нашей модели.

На трех диаграммах ниже показаны различия между первой и второй половиной временного ряда.

Помимо индивидуального визуального суждения, есть также несколько тестов для формальной проверки того, отличается ли распределение данных, а именно Тест Чоу, Двухвыборочный тест Колмогрова Смирнова и Дисперсионный анализ (ANOVA). . Для последних двух мы находим результаты теста между первой и второй половиной данных ниже.

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

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

Winsorizing

Winsorizing представляет собой действенную альтернативу простому удалению данных. Этот метод используется, когда вместо исключения потенциальных выбросов мы хотели бы сохранить их, но изменить их величину. Это делается путем указания, что потенциальный выброс должен быть скорректирован. Например, выигрышное значение X% означает, что все значения, которые выше, чем (100-X) процентиль, устанавливаются на значение (100-X) процентиля. Более подробное объяснение процесса можно найти здесь.

На диаграмме ниже показано влияние выигрыша наших данных до уровня 2,5% (значение, выбранное для выравнивания выбросов в первой половине данных). Если посмотреть на масштаб оси ординат, эффект меры становится очевидным, величина выбросов уменьшилась.

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

Обнаружение сезонности

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

На верхнем графике ниже мы видим фактическую целевую переменную, а именно количество случаев лихорадки денге в Сан-Хуане. График ниже показывает обсуждаемую автокорреляционную функцию этого временного ряда. График ACF показывает, что изначально (до лага 400) мы находим годичный образец. Это означает, что амплитуда кратна 52 (52 недели равны году). После отставания 400 мы обнаруживаем другую закономерность. То, что больше напоминает двухлетнюю сезонность. Учитывая вышеупомянутую более высокую важность более свежих данных, мы продолжим сезонность в 104 недели (или 2 года).

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

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

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

Учитывая, что деасонализированные временные ряды не страдают какой-либо нестационарностью, нам не нужно беспокоиться о количестве интегрирований в модели ARIMA, поскольку в них нет необходимости. Более важный вопрос заключается в том, сколько необходимо лагов авторегрессии (AR) и скользящих средних (MA). На практике длительность задержки процессов AR и MA обозначается буквами p и q соответственно.

Одним из хорошо известных подходов к определению оптимальной длины лага процессов AR и MA является так называемый метод Бокса-Дженкинса. Этот метод состоит из трех этапов:

  1. После обеспечения стационарности используйте графики функции автокорреляции (ACF) и функции частичной автокорреляции (PACF), чтобы сделать первоначальное предположение о длине запаздывания для процесса AR и MA.
  2. Используйте параметры, измеренные на первом шаге, и укажите модель
  3. Наконец, проверьте остатки подобранной модели на серийную корреляцию (независимость друг от друга) и на стационарность. Первое может быть достигнуто с помощью теста Юнга-Бокса. Если указанная модель не проходит эти тесты, вернитесь к шагу 2 и используйте немного другой параметр для p и / или q.

Следуя трем шагам, описанным выше, мы начинаем с построения автокорреляционной функции и частичной автокорреляционной функции. Это можно увидеть на графике ниже.

Частичная автоматическая корреляция - это степень корреляции между переменной и ее запаздыванием, которая не объясняется корреляциями на всех младших -пакетах. Это делается путем регрессии временного ряда со всеми лагами n-1 и, следовательно, их учета при оценке n-й частичной корреляции. Это контрастирует с функцией автокорреляции, которая не учитывает другие задержки.

Способ определения подходящей длины лага из двух приведенных выше графиков представлен в следующей таблице:

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

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

Важно подчеркнуть, что метод Бокса-Дженкинса основан на принципе экономии. Цитата оксфордского эконометриста Кевина Шеппарда:

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

Это означает, что если мы найдем несколько моделей, которые удовлетворяют нашим условиям в отношении последовательной корреляции, мы выберем модель с наименьшим количеством лагов. Поэтому мы проверяем результаты теста последовательной корреляции Ljung-Box не только для длин лагов 6 и 20 (что уже слишком много), но и для всех меньших длин. P-значения тестов Ljung-Box показаны на изображении ниже.

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

Проверка и прогнозирование по выборке

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

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

Весь код