Частотный
vs Байесовский
подход

Сравнительный анализ методов оценки структурной неопределённости на данных реального ВНП США (1909–1970) из оригинального набора Нельсона-Плоссера.

Центральный вопрос

«Являются ли макроэкономические шоки перманентными, или экономика всегда возвращается к устойчивому тренду?»

Python 3.13 PyMC 5.27 statsmodels Kalman Filter

Тест Зивота-Эндрюса

Тест определяет эндогенную дату структурного излома, оценивая ADF-статистику для каждой потенциальной точки разрыва. Минимальная (наиболее отрицательная) \(t\)-статистика указывает на год наибольшего структурного изменения.

t-статистика
-4.887
p-value
0.085
Год излома
1938
Лаги (AIC)
1
$$ \Delta y_t = c + \alpha\, y_{t-1} + \beta\, t + \gamma\, DU_t(\lambda) + \sum_{j=1}^{k} d_j\, \Delta y_{t-j} + e_t $$

где \(DU_t(\lambda) = 1\) если \(t > T\lambda\) (фиктивная переменная излома), \(\lambda \in (0.15,\, 0.85)\) -- параметр локализации.

Последовательность \(t\)-статистик Зивота-Эндрюса

Каждая точка -- ADF \(t\)-стат. для кандидатной даты излома. Минимум при 1938.

Критические значения

CV 1% -5.576 Не отвергнут
CV 5% -5.073 Не отвергнут
CV 10% -4.827 Отвергнут

«Серая зона»

На уровне 5% единичный корень не отвергается: шоки перманентны (Нельсон-Плоссер правы). На уровне 10% -- отвергается: экономика возвращается к тренду (Перрон прав).

В частотной парадигме это тупик: бинарный выбор без градации уверенности.

ln(Real GNP) со структурным изломом 1938

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

Классический парадокс

ARIMA Model B с детерминированным трендом и фиктивной переменной излома при 1938. Модель формально предпочтена по AIC/BIC, однако коэффициент излома статистически незначим.

Model B: ARMA(2,0) + Trend + Break

$$ y_t = \underbrace{\beta_0}_{4.626} + \underbrace{\beta_1}_{0.033}\, t + \underbrace{\beta_2}_{-0.093}\, D_{1938} + \phi_1\, u_{t-1} + \phi_2\, u_{t-2} + \varepsilon_t $$

где \(u_t = y_t - \beta_0 - \beta_1 t - \beta_2 D_{1938}\), \(\phi_1 = 1.302\), \(\phi_2 = -0.453\), \(\sigma^2 = 0.003\).

AIC (Model B)
-170.12
BIC (Model B)
-157.35
Break coef.
-0.093
p-value (break)
0.289
CI width
0.215

Fitted vs Actual: ARIMA Model B

Проблема: ригидность доверительного интервала

Ширина 95% CI составляет 0.2152 и не изменяется на протяжении всего ряда. Модель ARIMA оценивает дисперсию остатков как единственную глобальную константу \(\hat{\sigma}^2 = 0.003\). Она «размазывает» ошибку равномерно: интервал одинаков для стабильных 1950-х и турбулентных 1930-х. Это фундаментальное ограничение параметрической спецификации.

Байесовское решение

Локальный линейный тренд (LLT) моделирует уровень и наклон как случайные блуждания. Маргинализация латентных состояний через фильтр Калмана сводит задачу сэмплирования к 5 гиперпараметрам дисперсии.

Формулировка пространства состояний (State-Space)

Уравнение наблюдения

$$ y_t = \mu_t + \varepsilon_t, \quad \varepsilon_t \sim \mathcal{N}(0,\, \sigma^2_{\mathrm{obs}}) $$

Уравнения перехода

$$ \mu_t = \mu_{t-1} + \nu_{t-1} + \xi_t, \quad \xi_t \sim \mathcal{N}(0,\, \sigma^2_{\mathrm{level}}) $$ $$ \nu_t = \nu_{t-1} + \zeta_t, \quad \zeta_t \sim \mathcal{N}(0,\, \sigma^2_{\mathrm{slope}}) $$

Матричная форма

$$ \underbrace{\begin{bmatrix} \mu_t \\ \nu_t \end{bmatrix}}_{\mathbf{x}_t} = \underbrace{\begin{bmatrix} 1 & 1 \\ 0 & 1 \end{bmatrix}}_{\mathbf{F}} \mathbf{x}_{t-1} + \mathbf{w}_t, \qquad y_t = \underbrace{\begin{bmatrix} 1 & 0 \end{bmatrix}}_{\mathbf{H}} \mathbf{x}_t + \varepsilon_t $$
Probabilistic Programming

Маргинализация Калмана

Прямая параметризация (сэмплирование инноваций \(\xi_t, \zeta_t\)) приводит к 666 дивергенциям из-за «воронки Нила» (Neal's funnel) в 127-мерном пространстве латентных состояний.

Решение: аналитически интегрировать латентные состояния через фильтр Калмана, сведя задачу к 5 гиперпараметрам. Апостериорное распределение тренда восстанавливается сглаживателем Рауха-Тунга-Штрибеля (RTS).

Дивергенции
0
max R-hat
< 1.01
Цепи
4
Draws
2000

Маргинальное правдоподобие (Kalman filter log-likelihood)

$$ \log p(\mathbf{y} \mid \boldsymbol{\theta}) = -\frac{n}{2}\log 2\pi - \frac{1}{2}\sum_{t=1}^{n}\!\left[\log S_t + \frac{v_t^2}{S_t}\right] $$

где \(v_t = y_t - \hat{\mu}_{t|t-1}\) -- инновация, \(S_t = P^{(11)}_{t|t-1} + \sigma^2_{\mathrm{obs}}\) -- дисперсия инновации.

Реализация: PyMC + PyTensor Scan

phase3_bayesian_llt.py
# Kalman filter step (vectorised via pytensor.scan) def kf_step(y_t, mu_p, nu_p, P11_p, P12_p, P22_p, s_obs2, s_lev2, s_slp2): v = y_t - mu_p # innovation S = P11_p + s_obs2 # innovation variance K1 = P11_p / S # Kalman gain (level) K2 = P12_p / S # Kalman gain (slope) mu_f = mu_p + K1 * v # filtered level nu_f = nu_p + K2 * v # filtered slope P11_f = P11_p - K1 * P11_p P12_f = P12_p - K2 * P11_p P22_f = P22_p - K2 * P12_p ll_t = -0.5 * (pt.log(2*np.pi) + pt.log(S) + v**2/S) # Predict next step mu_p_next = mu_f + nu_f P11_p_next = P11_f + 2*P12_f + P22_f + s_lev2 P12_p_next = P12_f + P22_f P22_p_next = P22_f + s_slp2 return (mu_p_next, nu_f, P11_p_next, P12_p_next, P22_p_next, ll_t, mu_f, nu_f) # Marginal log-likelihood as PyMC Potential log_lik = results[5].sum() pm.Potential("kalman_loglik", log_lik)

Апостериорная сводка гиперпараметров

Параметр Mean SD HDI 2.5% HDI 97.5% ESS bulk R-hat
\(\sigma_{\mathrm{obs}}\) 0.0090 0.0070 0.0000 0.0220 2350 1.000
\(\sigma_{\mathrm{level}}\) 0.0570 0.0140 0.0260 0.0810 1889 1.000
\(\sigma_{\mathrm{slope}}\) 0.0180 0.0170 0.0000 0.0530 1321 1.000
\(\mu_0\) (level init) 4.7540 0.4390 3.9160 5.6280 5919 1.000
\(\nu_0\) (slope init) 0.0300 0.0500 -0.0680 0.1240 6734 1.000

Слабо информативные априорные распределения: \(\sigma \sim \mathrm{HalfCauchy}\), \(\mu_0 \sim \mathcal{N}(y_1, 0.5)\), \(\nu_0 \sim \mathcal{N}(0.03, 0.05)\). Сэмплер: NUTS (nutpie, Rust), 4 цепи, target_accept = 0.99.

Апостериорный тренд \(\hat{\mu}_t\) vs Observed

Стохастический наклон \(\nu_t\): динамика роста

Наклон естественно замедляется в Депрессию (\(\bar{\nu} = 0.023\), ~2.3%/год) и ускоряется после 1938 (\(\bar{\nu} = 0.036\), ~3.6%/год) -- рост в 1.58 раз.

Avg slope pre-1938
0.023
~2.3% annual
Avg slope post-1938
0.036
~3.6% annual
Acceleration
1.58x
Break detected by
Data, not
researcher

Синтез: калибровка неопределённости

Кульминация исследования: прямое сравнение ширины интервалов неопределённости в критический период экономической турбулентности (1938–1945).

CI vs HDI: наложение интервалов неопределённости

Серая полоса -- фокусный период 1938-1945. Обратите внимание на постоянную ширину CI vs адаптивную ширину HDI.

Ширина интервалов по годам (1938-1945)

Год ARIMA CI Bayes HDI Ratio
19380.21520.04870.23
19390.21520.04410.21
19400.21520.04520.21
19410.21520.04340.20
19420.21520.04350.20
19430.21520.04480.21
19440.21520.04730.22
19450.21520.04540.21
Mean0.21520.04530.21
HDI / CI Ratio
0.21x

Байесовский HDI в ~5 раз уже частотного CI

Почему HDI уже -- и почему это корректно

ARIMA CI: цена ригидности

Фиксированный параметрический тренд не может отследить нелинейную динамику перехода от Депрессии к WWII. Большие остатки раздувают оценку \(\hat{\sigma}^2\), создавая широкую, равномерную полосу, которая одновременно слишком широка в спокойные периоды и недостаточно адаптивна в период шоков.

Bayesian HDI: честная калибровка

Стохастический тренд плотно отслеживает данные. Неопределённость HDI отражает апостериорное распределение гиперпараметров \((\sigma_{\mathrm{obs}}, \sigma_{\mathrm{level}}, \sigma_{\mathrm{slope}})\), маргинализованных через фильтр Калмана -- а не мисфит жёсткой структуры. Интервалы не просто уже, а честнее откалиброваны.

Итог исследования

Кто был прав?

«Байесовский HDI в 5 раз уже частотного интервала. Это не значит, что мы занижаем риск. Это значит, что гибкая модель поглощает динамику шока, которую ARIMA оставляет в остатках. Интервалы уже, потому что они более честно откалиброваны -- концентрируя неопределённость там, где латентный тренд действительно неоднозначен.»

Фаза 1

Zivot-Andrews выявил излом в 1938, но результат попал в «серую зону» (p = 0.085).

Фаза 2-3

ARIMA навязывает жёсткий излом (p = 0.289, незначим). Байесовский LLT находит его органически.

Фаза 4

HDI/CI = 0.21. Не занижение риска, а лучшая калибровка.

Python 3.13· PyMC 5.27· statsmodels 0.14· nutpie 0.16 (Rust NUTS)· Nelson-Plosser 1909-1970