分位数回归与 Growth-at-Risk 尾部风险测度

* ============================================================
*  Stata: Quantile Regression & Growth-at-Risk 完整复现
* ============================================================
clear all
set seed 20240101
set obs 200

* ---------- 1. 模拟数据：金融条件 FCI -> GDP 增长 ----------
* FCI：正值=金融条件宽松，负值=收紧
gen t = _n
tsset t
gen fci = 0.8*l.fci + rnormal()            // FCI 自回归
replace fci = fci[_n-1] if _n==1

* GDP 增长：均值受 FCI 影响，条件方差也受 FCI 影响（脆弱增长机制）
gen sigma = exp(-0.8*fci)                   // FCI 越差 sigma 越大
gen gdp_growth = 0.5 + 0.6*fci + sigma*rnormal()
gen l_gdp = l.gdp_growth

* ---------- 2. 单分位点 QR：tau=0.05 ----------
qreg gdp_growth fci l_gdp, quantile(0.05) vce(bootstrap, reps=500)
est store q05
* 关注 fci 的系数：tau=0.05 时应显著为负且绝对值大

* ---------- 3. 多分位点联立 ----------
foreach tau in 0.05 0.25 0.50 0.75 0.95 {
    qreg gdp_growth fci l_gdp, quantile(`tau')
    est store q`=substr("`tau'",3,2)'
}
esttab q05 q25 q50 q75 q95, keep(fci l_gdp) se ///
    mtitles("p5" "p25" "p50" "p75" "p95") ///
    star(* 0.1 ** 0.05 *** 0.01)

* ---------- 4. Bootstrap 分位数回归（更可靠的标准误）----------
* bsqreg 用 bootstrap 估计 QR 的协方差矩阵
bsqreg gdp_growth fci l_gdp, quantile(0.05) reps(500)

* ---------- 5. 工具变量分位数回归 ivqreg ----------
* 若担心 FCI 内生（央行逆周期调整），可用工具变量
* ssc install ivqreg
* ivqreg gdp_growth l_gdp (fci = instrument), quantile(0.05)

* ---------- 6. 计算 GaR 时间序列（滚动估计）----------
gen gar_p5 = .
gen gar_p50 = .
gen gar_p95 = .
forvalues t = 60/200 {
    * 用前 t-1 期数据估计
    quietly qreg gdp_growth fci l_gdp in 1/`=`t'-1', quantile(0.05)
    predict resid_p5, xb
    qui replace gar_p5 = resid_p5[`t'] in `t'
    quietly qreg gdp_growth fci l_gdp in 1/`=`t'-1', quantile(0.50)
    predict resid_p50, xb
    qui replace gar_p50 = resid_p50[`t'] in `t'
    quietly qreg gdp_growth fci l_gdp in 1/`=`t'-1', quantile(0.95)
    predict resid_p95, xb
    qui replace gar_p95 = resid_p95[`t'] in `t'
    drop resid_p5 resid_p50 resid_p95
}

* ---------- 7. 画 GaR 时间序列 + 系数分位图 ----------
twoway (rarea gar_p5 gar_p95 t, fcolor(gs12) lcolor(none)) ///
       (line gar_p50 t, lcolor(navy) lwidth(med)) ///
       (line gdp_growth t, lcolor(red) lwidth(med)), ///
       title("GDP 增长分布与 Growth-at-Risk") ///
       legend(label(1 "5%-95% 区间") label(2 "中位数") label(3 "实际增长"))

* 系数沿分位变化图：手动循环估计并记录 fci 系数
gen tau_list = .
gen beta_fci = .
local i = 1
forvalues tau = 5(5)95 {
    quietly qreg gdp_growth fci l_gdp, quantile(`=`tau'/100')
    qui replace tau_list = `tau'/100 in `i'
    qui replace beta_fci = _b[fci] in `i'
    local i = `i' + 1
}
twoway connected beta_fci tau_list, ///
    title("FCI 系数沿 GDP 增长分布分位的变化") ///
    yline(0, lcolor(red) lpattern(dash))

# -*- coding: utf-8 -*-
"""
Growth-at-Risk (GaR) with Quantile Regression
=============================================
模拟金融条件 FCI 与 GDP 增长，做分位数回归，
构造 GaR (5% 分位) 时间序列，并与 OLS 均值预测对比。
"""
import numpy as np
import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt

np.random.seed(20240101)
T = 200

# ---------- 1. 模拟数据 ----------
# FCI：AR(1)，正值=宽松，负值=收紧
fci = np.zeros(T)
fci[0] = 0.0
for t in range(1, T):
    fci[t] = 0.8 * fci[t-1] + np.random.normal(0, 0.5)

# GDP 增长：均值随 FCI 上升；条件方差在 FCI 低（收紧）时放大
sigma = np.exp(-0.8 * fci)               # 金融条件越差，波动越大
gdp = np.zeros(T)
gdp[0] = 0.5
for t in range(1, T):
    gdp[t] = 0.5 + 0.6*fci[t-1] + 0.3*gdp[t-1] + sigma[t]*np.random.normal()

lag_gdp = np.r_[np.nan, gdp[:-1]]
df = pd.DataFrame({'gdp': gdp, 'fci': fci, 'lag_gdp': lag_gdp}).dropna()
X = sm.add_constant(df[['fci', 'lag_gdp']])
y = df['gdp']

# ---------- 2. OLS 均值回归（基准） ----------
ols = sm.OLS(y, X).fit()
print(ols.summary().tables[1])

# ---------- 3. 多分位数 QR ----------
quantiles = [0.05, 0.25, 0.50, 0.75, 0.95]
coefs = {}
for tau in quantiles:
    qr = sm.QuantReg(y, X).fit(q=tau)
    coefs[tau] = qr.params
    print(f"\n--- tau={tau:.2f} ---")
    print(qr.summary().tables[1])

# ---------- 4. GaR 时间序列：滚动估计 ----------
start = 60
gar5, gar50, gar95, actual = [], [], [], []
for t in range(start, len(df)):
    tr = df.iloc[:t]
    Xtr = sm.add_constant(tr[['fci','lag_gdp']])
    ytr = tr['gdp']
    Xte = sm.add_constant(df.iloc[[t]][['fci','lag_gdp']], has_constant='add')
    q5  = sm.QuantReg(ytr, Xtr).fit(q=0.05).predict(Xte).values[0]
    q50 = sm.QuantReg(ytr, Xtr).fit(q=0.50).predict(Xte).values[0]
    q95 = sm.QuantReg(ytr, Xtr).fit(q=0.95).predict(Xte).values[0]
    gar5.append(q5); gar50.append(q50); gar95.append(q95)
    actual.append(df['gdp'].iloc[t])

# ---------- 5. 样本外评估：GaR 5% 分位命中率 ----------
# 理论上 5% 的实际值应落在 GaR 之下
hit_rate = np.mean(np.array(actual) < np.array(gar5))
print(f"\nGaR(5%) 样本外命中率 = {hit_rate:.3f} （应接近 0.05）")

# ---------- 6. 可视化：FCI 系数沿分位变化 ----------
taus = np.linspace(0.05, 0.95, 19)
beta_fci = []
for tau in taus:
    beta_fci.append(sm.QuantReg(y, X).fit(q=tau).params['fci'])
plt.figure(figsize=(8,4.5))
plt.plot(taus, beta_fci, 'o-', label='QR: beta(FCI)')
plt.axhline(ols.params['fci'], color='red', ls='--', label='OLS mean')
plt.xlabel('Quantile'); plt.ylabel('Coefficient on FCI')
plt.title('FCI 对 GDP 增长的效应沿分位变化')
plt.legend(); plt.tight_layout()
# plt.savefig('gar_coef.png', dpi=150)
plt.show()

# ---------- 7. 输出：GaR 在 FCI 收紧时应显著下移 ----------
print("\n样本外 RMSE:")
print(f"  均值预测(OLS) = {np.sqrt(np.mean((np.array(actual)-np.array(gar50))**2)):.3f}")