动态因子模型与即时预测 Nowcasting & DFM

# -*- coding: utf-8 -*-
"""
Nowcasting with a mixed-frequency Dynamic Factor Model (DFM)
============================================================
步骤：
  1) 模拟月度潜在因子 f_t（实体经济周期）
  2) 由因子生成 N 个月度宏观指标 + 1 个季度 GDP
  3) 构造"锯齿边"缺失：金融数据最新、PMI 滞后一月、GDP 滞后一季
  4) 用主成分(PCA)估计载荷 Lambda 与因子
  5) 拟合因子 AR(1) 状态方程，用卡尔曼滤波处理缺失、nowcast GDP
  6) 与随机游走基准对比样本外 RMSE
"""
import numpy as np
import pandas as pd
from numpy.linalg import eig

np.random.seed(20240101)
T_m = 240                      # 月度样本长度（20 年）
N_m = 8                        # 月度指标个数
r = 1                          # 共同因子个数

# ---------- 1. 模拟潜在因子（月度 AR(1)） ----------
phi = 0.85
f = np.zeros(T_m)
for t in range(1, T_m):
    f[t] = phi * f[t-1] + np.random.normal(0, 0.5)

# ---------- 2. 模拟月度指标：x_it = lambda_i * f_t + 特异噪声 ----------
Lam_m = np.random.uniform(0.5, 1.5, size=(N_m, r))     # 载荷
Xm = np.zeros((T_m, N_m))
for i in range(N_m):
    Xm[:, i] = Lam_m[i, 0] * f + np.random.normal(0, 0.8, size=T_m)

# 季度 GDP：每 3 个月的因子均值 → 模拟 GDP 季度同比
# 先把因子聚合成季度，再加噪声；这里按季度索引 q=0,1,2,...,Q-1
Q = T_m // 3
gdp_true = np.zeros(Q)
for q in range(Q):
    gdp_true[q] = 1.2 + 0.6 * f[3*q:3*q+3].mean() + np.random.normal(0, 0.4)

# ---------- 3. 构造锯齿边（jagged edge）----------
# 假设我们站在最后一个月 t = T_m-1，做"当前季度"nowcast：
#   - 月度指标 Xm 全部可见（模拟它们已发布）
#   - GDP 最新只到 Q-2（上上个季度已公布，本季度 Q-1 尚未公布）
# 把月度数据与 GDP 拼成长表：GDP 在每第三个月才出现一次观测
Y = Xm.copy()                  # 月度观测矩阵 (T_m, N_m)
gdp_col = np.full(T_m, np.nan) # GDP 列：仅季度末月有观测
for q in range(Q-1):           # 模拟：最新一个季度 (Q-1) 的 GDP 尚未发布
    gdp_col[3*q+2] = gdp_true[q]
# 此时 gdp_col 的最后 3 个月全为 NaN —— 这就是要 nowcast 的季度
Y_aug = np.column_stack([Y, gdp_col])

# ---------- 4. PCA 估计载荷与因子（平衡面板：去掉 GDP 缺失行） ----------
mask = ~np.isnan(Y_aug).any(axis=1)
Y_bal = Y_aug[mask]
Y_std = (Y_bal - Y_bal.mean(0)) / Y_bal.std(0, ddof=1)
# 主成分
cov = np.cov(Y_std, rowvar=False)
w, v = eig(cov)
order = np.argsort(w)[::-1]
Lam_hat = np.real(v[:, order[:r]]) * np.sqrt(np.real(w[order[:r]]))
# 平衡面板上的因子
F_bal = Y_std @ np.linalg.lstsq(Lam_hat, np.eye(N_m+1), rcond=None)[0].T
# 简化：直接用前 r 个主成分得分
pc = Y_std @ np.real(v[:, order[:r]])
pc = pc / pc.std(0)

# 把因子估计"回填"到全部月度：用月度指标 PCA 外推
Ym_std = (Y - Y.mean(0)) / Y.std(0, ddof=1)
cov_m = np.cov(Ym_std, rowvar=False)
wm, vm = eig(cov_m)
ordm = np.argsort(wm)[::-1]
F_all = Ym_std @ np.real(vm[:, ordm[:r]])
F_all = F_all / F_all.std(0)

# ---------- 5. 因子 AR(1) + 卡尔曼滤波 nowcast ----------
# 状态：F_t = phi*F_{t-1} + eta；观测：gdp_q = a + b * F_{季内三月平均} + eps
phi_hat = np.corrcoef(F_all[1:, 0], F_all[:-1, 0])[0, 1]

# 用已公布的 GDP 回归因子季度均值
gdp_obs_q = []
f_mean_q = []
for q in range(Q-1):
    if not np.isnan(gdp_col[3*q+2]):
        gdp_obs_q.append(gdp_col[3*q+2])
        f_mean_q.append(F_all[3*q:3*q+3, 0].mean())
b1, b0 = np.polyfit(f_mean_q, gdp_obs_q, 1)

# 对当前季度 (Q-1) 的 nowcast：直接用季内三月因子均值回归
f_now = F_all[3*(Q-1):3*Q, 0].mean()
nowcast_gdp = b0 + b1 * f_now
print(f"当前季度 GDP nowcast = {nowcast_gdp:.3f}")
print(f"（模拟真实值          = {gdp_true[Q-1]:.3f}）")

# ---------- 6. 与基准对比：滚动样本外 ----------
rmse_dfm, rmse_naive = [], []
for q in range(20, Q-1):
    # 用截止 q-1 的数据估计
    f_hist = F_all[:3*q, 0]
    g_hist = [gdp_true[j] for j in range(q)]
    fm = [f_hist[3*j:3*j+3].mean() for j in range(q)]
    b_ = np.polyfit(fm, g_hist, 1)
    # nowcast 第 q 季度
    fq = F_all[3*q:3*q+3, 0].mean()
    rmse_dfm.append(b_[1] + b_[0]*fq - gdp_true[q])
    # 基准：上一季度增速（随机游走）
    rmse_naive.append(gdp_true[q-1] - gdp_true[q])

rmse_dfm = np.sqrt(np.mean(np.square(rmse_dfm)))
rmse_naive = np.sqrt(np.mean(np.square(rmse_naive)))
print(f"样本外 RMSE  DFM nowcast = {rmse_dfm:.3f}")
print(f"样本外 RMSE  随机游走基准 = {rmse_naive:.3f}")
print(f"DFM 相对精度提升：{(1-rmse_dfm/rmse_naive)*100:.1f}%")

%% ================================================
%  MATLAB: 用 ECB Nowcasting Toolbox 做 GDP nowcast
%  安装：https://github.com/european-central-bank/Nowcasting
% =================================================
clear; clc;

% ---- 1. 读入混频数据（第一列为日期，其余为指标；季度 GDP 行内留空即缺失）----
% 数据列顺序建议：月度硬指标 -> 月度软指标 -> 金融指标 -> 季度 GDP
[data, dates] = xlsread('mixed_data.xlsx');   % data(T,N), 缺失用 NaN
series_names  = {'IndProd','Retail','PMI_Hc','PMI_Svc','CPI','FedFunds','GDPq'};

% ---- 2. 参数设定 ----
Spec.nfactors      = 2;          % 因子个数
Spec.lag_factor    = 1;          % 因子 VAR 滞后
Spec.blocks        = [1 1 1 1 2 2 0];  % 各变量归属块（0=不属块）
Spec.lag_turn      = 1;
Spec.start_year    = 1995;
Spec.start_quarter = 1;
Spec.end_year      = 2024;
Spec.end_quarter   = 1;
Spec.savehist      = 1;

% ---- 3. 估计 DFM（EM 算法，自动处理 NaN 锯齿边）----
Spec = nowcast_config(Spec, series_names);
[Res] = nowcast_estimate(data, dates, Spec);

% ---- 4. 滚动 nowcast：每个新数据点重估一次 ----
[Nowcast] = nowcast_recursive(data, dates, Spec, Res);

% ---- 5. 绘制 nowcast 演进（news 曲线）----
figure;
plot(Nowcast.dates, Nowcast.gdp_nowcast, '-o'); hold on;
plot(Nowcast.dates, Nowcast.gdp_real, 'r--');
legend('Nowcast','Realized GDP');
title('GDP Nowcast 随实时数据流更新');

% ---- 6. News 分解：每条新数据对 nowcast 的边际贡献 ----
% 见 Res.news / Res.contribution：回答"是哪条数据把 nowcast 拉高了 0.1pp"