Фильтр ХП: плюсы, минусы, прокачка

Julia
Автор

Е. Тымченко

Дата публикации

23 сентября 2024 г.

Если добавить линию с общей тенденцией, становится понятнее.

О чём статья?

В анализе временных рядов мы часто выделяем тренды. Так мы решаем множество задач, которые я бы разбил их на 2 группы:

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

С первой задачей, неплохо справляется фильтр Ходрика-Прескотта, но иногда от инструмента требуют то, чего он дать не способен. Фильтром Ходрика-Прескотта измеряют разрыв выпуска.

shame

Отличная статья об этом есть у мегабазированного J.D. Hamilton. Несмотря на байтовое название Никогда не используйте фильтр ХП, Хэмилтон, последовательно описывает ограничения фильтра ХП, с которыми должен быть знаком практик. В статье по ссылке есть некоторые опечатки, в частности в формуле 1, где постулируется оптимизационная задача для фильтра ХП. Но это не так важно.

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

Достоинства фильтра ХП:

  • Очень прост в использовании.
  • Очень прост в реализации.
  • Действительно хорошо улавливает тенденции.

Недостатки фильтра ХП:

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

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

Классический двухсторонний фильтр ХП:

Классический фильтр ХП призван решить задачу подбора тренда \(\hat{y}\) к исходному ряду \(y\) так, чтобы минимизировать:

\[\min_{\hat{y}} \sum_{t = 1}^T (y_t - \hat{y}_t) + \lambda \times \sum_{t = 2}^{T-1}\left[ \hat{y}_{t+1} - \hat{y}_t) - (\hat{y}_t - \hat{y}_{t-1}) \right]\]

Он называется двухсторонним, потому что для момента \(t\) он использует данный слева \(t-1\) и данные справа \(t+1\). Есть также и односторонняя версия, которая использует только уже реализовавшиеся данные \(y\) для расчёта \(\hat{y}\).

Идея первого слагаемого в том, чтобы \(\hat{y}_t\) был как можно ближе к исходному ряду. Идея второго слагаемого - обеспечить гладкость тренда, чтобы он не слишком сильно менялся от наблюдения к наблюдению. Чем больше \(\lambda\) - тем больше веса при оптимизации мы вкладываем в гладкость. Оптимальные параметры \(\lambda\), которые часто даёт литература:

  • Для месячных данных: 14400 (или 129600),
  • Для квартальных данных: 1600,
  • Для годовых данных: 100 или (6.25).

Они примерно соответствуют идее, что тренд может значительно менять направление не чаще, чем раз в 8 лет. В зависимости от рядов и задач, значение \(\lambda\) может отличаться от приводимых. С этим связана значительная часть критики, которой подвергается фильтр ХП.

HP фильтр и нелинейная численная оптимизация в Julia

Чтобы запрограммировать фильтр ХП мы используем Julia, потому что это быстрый язык, на котором весело программировать. В отличие от всех библиотек, что я находил, мы не будем выводить closed-form решение, а оптимизируем целевую (нелинейную) функцию численно. Решение такой задачи в Julia занимает меньше миллисекунды, а алгоритм выглядит лаконично.

Библиотеки

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

using CSV, DataFrames, JuMP, Ipopt, Plots, Dates, LaTeXStrings, HypothesisTests

Данные

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

  • ВВП в постоянных ценах, тренд которого всегда положительный
  • Инфляционных ожиданиях из опросов ФОМ, тренд которых со временем может меняться.
data = CSV.read(download("https://raw.githubusercontent.com/ETymch/Econometrics_2023/refs/heads/main/Datasets/data_hp.csv"), DataFrame);
data.output = data.output ./ 10^3; # Переводим миллиарды в триллионы

Реализация двухстороннего фильтра ХП:

function hpfilter(y, λ) # y = исходный ряд,  lambda = параметр сглаживания
    # Модель
    model = Model(
        optimizer_with_attributes(
            Ipopt.Optimizer)
        ); # Создаём модель для оптимизации
    set_silent(model) # выключаем отображение результатов (можете убрать, тогда при оптимизации будет показана статистика оптимизации).

    # Размерности
    T = length(y) # длина ряда
    
    # Переменные

    @variable(model, yhat[1:T]) # Параметра для оценки - тренд ХП.

    # Целевая функция, как в формуле, которую я привёл:
    
    @NLobjective(model, Min,
    sum((y[t] - yhat[t])^2 for t in 1:T) + λ * sum(((yhat[t+1] - yhat[t]) - (yhat[t] - yhat[t-1]))^2 for t in 2:T-1)
    )

    optimize!(model) # Оптимизация
    
    # Результаты - в словарь
    
        yh = JuMP.value.(yhat);
        results = Dict(
            "yhat" => yh)
        return results
end
hpfilter (generic function with 1 method)

Поскольку мы используем алгоритмы численной оптимизации, решение будет чуть-чуть (на уровне 7 знака после точки) отличаться от того, которое мы получим из готовых библиотек. Последние используют готовое (closed form) решение оптимизационной задачки, т.е. они выписывают градиент и всё такое. На мой взгляд, попытка решения такой задачки, выписывание матрицы с градиентом вызывает у людей грусть. Вместо этого, пока одни берут производные, численные методы позволяют достаточно быстро запрограммировать всё что угодно.

Давайте применим алгоритм и сделаем визуализацию:

begin
    y_hp = hpfilter(data.output, 100); # наша функция
    yhat = y_hp["yhat"] # тренд
    
    trends = plot(data.date, data.output, label = "Исходный ряд") # График c исходными данными
    plot!(trends, data.date, yhat, label = "HP-тренд") # В Julia функции с ! изменяют уже созданные объекты. В данном случае мы добавляем на график слой с трендом.
    cycle = plot(data.date, data.output .- yhat, label = "Цикл") # график с циклом
    plot!(cycle, data.date, fill(0.0, nrow(data)), color = :black, linestyle = :dash, label = :none) # выделим дополнительно пересечение с 0.
    plot(trends, cycle, layout = (2,1), dpi = 300) # Выкладка. В R и ggplot2 для этого есть библиотека patchwork.
end

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit https://github.com/coin-or/Ipopt
******************************************************************************

Неоднородность фильтра ХП

Проблема фильтра ХП ещё и в том, что принцип, по которому рассчитывается значение тренда для начала, серединки и конца временного ряда - разный. В этом просто убедиться, если взять производную целевой функции для разных параметров \(\hat{y}_t\):

  • Для \(t = 1\):

\[y_1 = (1 + \lambda) \hat{y}_1 - 2 \lambda \hat{y}_2 + \lambda \hat{y}_3\]

  • Для \(t = 2\):

\[y_2 = -2 \lambda \hat{y}_1 + (1 - 5\lambda) \hat{y}_2 - 4 \lambda \hat{y}_3 + \lambda \hat{y}_4\]

  • Для \(3\leq t \leq T-2\):

\[y_t = \lambda \hat{y}_{t-2} - 4 \lambda \hat{y}_{t-1} + (1 + 6 \lambda) \hat{y}_t - 4 \lambda \hat{y}_{t+1} + \lambda \hat{y}_{t+2}\]

  • Для \(t = T-1\):

\[y_{T-1} = \lambda \hat{y}_{T-3} - 4 \lambda \hat{y}_{T-2} + (1 + 5\lambda) \hat{y}_{T-1} - 2\lambda \hat{y}_T\]

  • Для \(t = T\):

\[y_T = \lambda \hat{y}_{T-2} - 2\lambda \hat{y}_{T-1} + 2\lambda \hat{y}_T\]

Поэтому, если вам показывают актуальные, только что вышедшие статистические данные, рисуют при помощи двухстороннего фильтра ХП тренд и анализируют результаты, то это не вполне верно. Например, в библиотеке mFilter в R реализована только двухсторонняя версия фильтра ХП. Его односторонняя версия не столь популярна, хотя технически эту проблему она решает.

Односторонний фильтр ХП:

Можно также закодить односторонний фильтр. Односторонний фильтр решает проблему двухстороннего. Теперь производные для всех \(\hat{y}_t\) выглядят одинаково.

\[\min_{\hat{y}} \sum_{t = 1}^{T-1} (y_t - \hat{y}_t) + \lambda \times \sum_{t = 2}^{T-2}\left[ \hat{y}_{t+1} - \hat{y}_{t}) - (\hat{y}_{t} - \hat{y}_{t-1}) \right]\]

Алгоритм:

begin
    function hpfilter1(y, λ)
    # Модель
    model = Model(
        optimizer_with_attributes(
            Ipopt.Optimizer)
        );
    set_silent(model)

    # Размерности
    T = length(y)
    
    # Переменные

  @variable(model, yhat[1:T])

    @NLobjective(model, Min,
    sum((y[t] - yhat[t])^2 for t in 1:T-1) + λ * sum(((yhat[t+1] - yhat[t]) - (yhat[t] - yhat[t-1]))^2 for t in 2:T-2)
    )

    optimize!(model)
        yh = JuMP.value.(yhat);
        results = Dict(
            "yhat" => yh)
        return results
    end
end;

hp2

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

hp2

Прокаченный фильтр ХП

boost

В последние годы по-прежнему выходит много литературы на тему фильтра ХП. В работе P. Phillips и Zhentao Shi немного прокачали фильтр, добавили дополнительные условия в постановку задачки. У прокаченной версии фильтра ХП есть github repo с функцией BoostedHP, переведённой на все популярные языки: Python, R, Julia, Matlab и др. Давайте запрограммируем и его тоже.

Алгоритм тут очень простой:

  1. Берём временной ряд. Применяем к нему фильтр ХП.
  2. Проверяем циклическую компоненту на стационарность. Авторы используют тест Дики-Фуллера.
  3. Если она стационарна, -> мы достигли успеха.
  4. Если она не стационарна, применяем к тренду ХП фильтр ХП, получаем новое значение цикла \(\hat{c}_2\).
  5. Повторяем 2-4 пока циклическая компонента не станет стационарной.

Вот и вся идея. В Julia можно это лаконично запрограммировать:

function BHPFilter(y, λ, sig, Max_Iter)

T = length(y)
order = Int(floor((T-1)^(1/3)))
    
trend = hpfilter(y, λ)["yhat"]
cycle = y .- trend

# Тест на наличие тренда    
p_val = pvalue(ADFTest(cycle, :trend, order)) # функция ADFtest выполняет тест на наличие тренда, потом сразу достаём pvalue.
i = 1   
while (p_val > sig) & (i <= Max_Iter)
    trend = y .- cycle # Обновление тренда
  trend = hpfilter(trend, λ)["yhat"] # Применяем к нему фильтр ХП
    cycle = y .- trend # обновлённый цикл
    p_val = pvalue(ADFTest(cycle, :trend, order)) # Вычисляем p-value в тесте Д-Ф.
    i = i + 1 # Следующая итерация
end
    return trend
end
BHPFilter (generic function with 1 method)

Результаты можно сравнить с двухсторонним фильтром ХП:

bhp

Если ваша задача: по-умному остационарить временной ряд (т.е. действительно вычислить стохастический тренд), то прокаченный фильтр ХП - хороший вариант.

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

tu