大数据宏观预测 Big Data Macroeconomic Forecasting

# -*- coding: utf-8 -*-
"""
大数据宏观预测：滚动窗口下比较 Lasso / Ridge / RF / XGBoost
作者: 经济学研究工作站
参考: Medeiros et al. (2021, JBES); 郑挺国等 (2025, 管理科学学报)
"""
import numpy as np
import pandas as pd
from sklearn.linear_model import LassoCV, RidgeCV, ElasticNetCV
from sklearn.ensemble import RandomForestRegressor
from sklearn.preprocessing import StandardScaler
from scipy import stats
import warnings
warnings.filterwarnings("ignore")

np.random.seed(20260911)

# ================================================================
# 1. 模拟"数据丰富"宏观数据：T=480 个月, p=50 个预测子
#    真实 DGP: y 由 5 个变量 + 1 个非线性交互项驱动,
#              其余 45 个预测子是噪声 (模拟 FRED-MD 式高维数据)
# ================================================================
T, p = 480, 50
X = np.random.randn(T, p) * 0.5

# y 的真实生成: AR(1) + 5 个线性项 + 1 个非线性门限项
y = np.zeros(T)
for t in range(1, T):
    lin = 0.7 * X[t, 0] + 0.5 * X[t, 1] - 0.4 * X[t, 2] \
          + 0.3 * X[t, 3] + 0.2 * X[t, 4]
    # 非线性: 只有当 X[t,5] > 0 时 X[t,6] 才起作用 (门限效应)
    nonlinear = 1.2 * X[t, 6] * (X[t, 5] > 0)
    y[t] = 0.8 * y[t-1] + lin + nonlinear + np.random.randn() * 0.5

# ================================================================
# 2. 构造 h=1 步预测的特征矩阵 (用 t 时刻信息预测 t+1)
#    标准做法: y_{t+1} ~ X_t  (不能用 y_{t+1} 的同期 X, 那是数据泄露)
# ================================================================
X_feat = X[:-1]            # 用 t 时刻的 X
y_target = y[1:]           # 预测 y_{t+1}

# ================================================================
# 3. 滚动窗口样本外预测
#    初始训练窗口 240 期, 每次滚动 1 期, 预测下一期
# ================================================================
WARMUP = 240
oos_index = np.arange(WARMUP, len(y_target))

def rolling_predict(model_factory, Xf, yt, warmup):
    """通用滚动窗口预测函数"""
    preds = np.zeros(len(yt) - warmup)
    for i, t in enumerate(range(warmup, len(yt))):
        X_tr, X_te = Xf[:t], Xf[t:t+1]
        y_tr = yt[:t]
        # 必须在训练窗口内 fit scaler, 不能用全样本均值!
        scaler = StandardScaler().fit(X_tr)
        X_tr_s = scaler.transform(X_tr)
        X_te_s = scaler.transform(X_te)
        model = model_factory()
        model.fit(X_tr_s, y_tr)
        preds[i] = model.predict(X_te_s)[0]
    return preds

# --- 模型工厂 (每次滚动重新拟合, 避免未来信息) ---
models = {
    "AR(1)基准":   lambda: LassoCV(cv=5, alphas=[0.01]),  # 简化: 仅截距+滞后项近似
    "Ridge":      lambda: RidgeCV(alphas=np.logspace(-2, 2, 10)),
    "Lasso":      lambda: LassoCV(cv=5, n_alphas=50, max_iter=5000),
    "ElasticNet": lambda: ElasticNetCV(cv=5, l1_ratio=[0.2,0.5,0.8,0.95],
                                       n_alphas=30, max_iter=5000),
    "RandomForest": lambda: RandomForestRegressor(
            n_estimators=500, max_depth=8, min_samples_leaf=5,
            max_features="sqrt", random_state=42, n_jobs=-1),
}

# --- 为了让 AR(1) 基准真正是 AR(1), 单独构造 (只用 y 的滞后项) ---
y_lag1 = np.roll(y_target, 1); y_lag1[0] = y[0]
X_ar = y_lag1.reshape(-1, 1)

results = {}
for name, factory in models.items():
    print(f"Rolling-window OOS: {name} ...")
    if name == "AR(1)基准":
        results[name] = rolling_predict(factory, X_ar, y_target, WARMUP)
    else:
        results[name] = rolling_predict(factory, X_feat, y_target, WARMUP)

# ================================================================
# 4. 计算 RMSE / MAE / 相对 RMSE (vs AR(1))
# ================================================================
y_oos = y_target[WARMUP:]
rmse = {k: np.sqrt(np.mean((v - y_oos)**2)) for k, v in results.items()}
mae  = {k: np.mean(np.abs(v - y_oos))       for k, v in results.items()}

base_rmse = rmse["AR(1)基准"]
print("\n{:15s} {:>10s} {:>10s} {:>10s}".format("模型", "RMSE", "MAE", "相对RMSE"))
for k in results:
    print(f"{k:15s} {rmse[k]:10.4f} {mae[k]:10.4f} {rmse[k]/base_rmse:10.4f}")

# ================================================================
# 5. Diebold-Mariano 检验 (以 AR(1) 为基准)
#    H0: 模型 k 与 AR(1) 的预测精度相同
# ================================================================
def dm_test(e_a, e_b, h=1):
    """Diebold-Mariano 检验, 平方损失, h 步预测需 HAC 调整"""
    d = e_a**2 - e_b**2                     # 损失差
    T = len(d)
    d_bar = d.mean()
    # Newey-West HAC 长期方差 (h 步预测有 h-1 阶自相关)
    gamma0 = np.var(d, ddof=1)
    lrvar = gamma0
    for lag in range(1, h):
        gamma_lag = np.cov(d[lag:], d[:-lag], ddof=1)[0,1]
        lrvar += 2 * (1 - lag/h) * gamma_lag
    dm_stat = d_bar / np.sqrt(lrvar / T)
    p_value = 2 * (1 - stats.norm.cdf(abs(dm_stat)))
    return dm_stat, p_value

print("\nDiebold-Mariano 检验 (基准 = AR(1), 负 DM 表示模型更优):")
e_ar = results["AR(1)基准"] - y_oos
for k in ["Ridge", "Lasso", "ElasticNet", "RandomForest"]:
    e_k = results[k] - y_oos
    dm, p = dm_test(e_k, e_ar, h=1)
    sig = "***" if p < 0.01 else ("**" if p < 0.05 else ("*" if p < 0.1 else ""))
    print(f"  {k:15s} DM = {dm:7.3f}  p = {p:6.4f} {sig}")

# ================================================================
# 6. 变量重要性 (用最后一个训练窗口的 RF 输出)
# ================================================================
scaler = StandardScaler().fit(X_feat[:WARMUP])
rf_final = RandomForestRegressor(n_estimators=500, max_depth=8,
                                 min_samples_leaf=5, random_state=42)
rf_final.fit(scaler.transform(X_feat[:WARMUP]), y_target[:WARMUP])
imp = pd.Series(rf_final.feature_importances_,
                index=[f"x{i}" for i in range(p)]).sort_values(ascending=False)
print("\nRF 变量重要性 Top 10:")
print(imp.head(10).round(4))
# 期望看到 x0..x6 排前几名 (因为 DGP 就是这 7 个变量驱动的)

*==============================================================*
* Stata 大数据宏观预测: Elastic Net + 滚动窗口
* 命令: elasticnet (Stata 16+ 内置, 无需 ssc install)
* 参考: Medeiros et al. (2021 JBES)
*==============================================================*
clear all
set more off
set seed 20260911

* --- 1. 构造模拟宏观数据: T=480, p=50 ---
set obs 480
gen t = _n
forvalues j=1/50 {
    gen x`j' = rnormal()*0.5
}
gen y = .
replace y = 0 in 1
forvalues t=2/480 {
    local lin = 0.7*x1[`t'] + 0.5*x2[`t'] - 0.4*x3[`t'] ///
              + 0.3*x4[`t'] + 0.2*x5[`t']
    local nonlin = 1.2*x6[`t']*(x7[`t']>0)
    replace y = 0.8*y[`t'-1] + `lin' + `nonlin' + rnormal()*0.5 in `t'
}

* --- 2. Elastic Net 主命令 (在全样本上做 CV, 仅作演示) ---
* 语法:
*   elasticnet  y  x1-x50, l1options(...) cv(5)
*   l1options(alpha(0.5)) 表示 Elastic Net; alpha(1)=Lasso, alpha(0)=Ridge
*   可选: selection(cv) 用 5 折 CV 选 lambda
elasticnet y x1-x50, l1options(alpha(0.5)) cv(5) selection(cv)

* 查看 CV 结果: 哪条 lambda 最优
elasticnet, cvplot
* 列出被选中的变量 (系数非零)
elasticnet, lassocoef display(coef)

* --- 3. 滚动窗口样本外预测 ---
* 初始训练: t=1..240, 测试: t=241..480
gen yhat_en = .
forvalues t=241/479 {
    quietly {
        * 在 [1, t] 上训练
        elasticnet y x1-x50 if t <= `t', ///
            l1options(alpha(0.5)) cv(5) selection(cv)
        * 预测 t+1
        predict yhat_temp if t==`t', xb
        replace yhat_en = yhat_temp[_n+1] if t==`t'
        drop yhat_temp
    }
    display "Rolling OOS: t = " `t'
}

* --- 4. 计算 RMSE / MAE (样本外) ---
gen err = y - yhat_en if t >= 242 & !missing(yhat_en)
quietly summarize err
gen rmse_en = sqrt(mean(err^2))
gen mae_en  = mean(abs(err))
display "Elastic Net OOS RMSE = " rmse_en
display "Elastic Net OOS MAE  = " mae_en

* --- 5. 对照: 简单 AR(1) 基准 ---
gen y_lag = L.y
reg y y_lag if t <= 240
predict yhat_ar if t >= 241, xb
gen err_ar = y - yhat_ar if t >= 241 & !missing(yhat_en)
quietly summarize err_ar
gen rmse_ar = sqrt(mean(err_ar^2))
display "AR(1)    OOS RMSE = " rmse_ar
display "相对 RMSE (EN/AR) = " rmse_en / rmse_ar

* --- 6. Python 集成: 在 Stata 里直接调 sklearn 做 RF/XGBoost ---
* Stata 16+ 支持 python 命令, 可直接调用 Python 环境
python:
import pandas as pd, numpy as np
from sklearn.ensemble import RandomForestRegressor
from sklearn.preprocessing import StandardScaler
# sfi (Stata Function Interface) 把当前 Stata 数据读进 pandas
from sfi import Data
df = pd.DataFrame({v: Data.get(v) for v in Data.getVarNames()})
X = df[[f"x{i}" for i in range(1,51)]].values
y = df["y"].values
# 滚动窗口 RF
WARMUP = 240
preds = np.full(len(y), np.nan)
for t in range(WARMUP, len(y)-1):
    sc = StandardScaler().fit(X[:t])
    rf = RandomForestRegressor(n_estimators=300, max_depth=8,
                               random_state=42, n_jobs=-1)
    rf.fit(sc.transform(X[:t]), y[:t])
    preds[t+1] = rf.predict(sc.transform(X[t:t+1]))[0]
# 把结果写回 Stata
Data.addVarDouble("yhat_rf")
Data.store("yhat_rf", None, preds.tolist())
end

* 计算 RF 相对 RMSE
gen err_rf = y - yhat_rf if t >= 242 & !missing(yhat_rf)
quietly summarize err_rf
display "Random Forest OOS RMSE = " sqrt(mean(err_rf^2))