Python复现Gertler-Karadi 2015 External Instrument VAR模型时脉冲响应发散

问题描述 #

我正尝试在Python中复现Gertler-Karadi 2015的外部工具VAR模型，但脉冲响应计算出现数值爆炸（发散）的异常结果，调试后仍未解决，不确定是否需要对数据进行缩放或归一化。

脉冲响应计算代码 #

# Define variables (based on your data)
num_impulses = 40  # number of periods
# based on the given matrices
col_names = ['IP', 'HICP', 'de1y', 'eurostoxx50', 'bbb_spread']  # example column names

# Initialize the impulse response storage
irs = np.zeros([num_lags + num_impulses, num_vars])
# Set initial impulse response (based on the third column of A0)
irs[num_lags, :] = A0[:, 2]  # Correctly initializing with the third column of A0

# Iterate to calculate the impulse responses
for jj in range(1, num_impulses):
    # Extract lagged values for the current impulse response calculation
    lvars = irs[np.arange(start=num_lags + jj - 1, stop=jj - 1, step=-1), :].T
    # Compute the next row of the impulse response
    irs[num_lags + jj, :] = lvars.flatten().T@ Bhat  # Ensure correct dimensions and transpose Bhat

# Slice the relevant impulse response values
irs = irs[num_lags:num_lags + num_impulses, :]

# Convert to DataFrame
irf_chol = pd.DataFrame(irs, columns=col_names)
print(irf_chol.head())

运行异常结果 #

IP      HICP        de1y  eurostoxx50  bbb_spread
0   0.000000  0.000000    1.000000     0.412065   -0.355296
1  -0.025553 -0.009461   -0.560410    -1.062703    6.280475
2   0.692767 -0.093960  -20.098512     4.914833   16.835246
3   0.342356  0.788316   91.130221    56.952628 -616.248000
4 -47.819026  0.379773  755.920970  -438.207420  647.235308

排查与解决方向 #

• 数据缩放/归一化：VAR模型对变量量纲高度敏感，你的变量包含工业生产指数、通胀率、利率、股票指数和利差，单位差异极大，系数估计会被量纲大的变量主导，直接导致脉冲响应发散。建议先对所有变量做标准化处理（将每个变量转换为均值0、方差1的形式），再重新估计模型的A0和Bhat参数，之后再计算脉冲响应。
• 模型稳定性验证：VAR模型稳定的核心是特征方程所有根的模长小于1。可以通过以下代码检查：

  import numpy as np
  # 构造VAR的伴随矩阵
  k = num_vars
  p = num_lags
  companion_matrix = np.zeros((k*p, k*p))
  companion_matrix[:k, :] = Bhat.T
  if p > 1:
      companion_matrix[k:, :-k] = np.eye(k*(p-1))
  # 计算特征值并检查模长
  eigenvalues = np.linalg.eigvals(companion_matrix)
  if any(np.abs(eigenvalues) >= 1):
      print("VAR模型不稳定，存在模长≥1的特征值")
  else:
      print("VAR模型稳定")
  

  若模型不稳定，需调整滞后阶数num_lags，或检查数据是否存在单位根（可通过ADF检验），考虑对非平稳变量做差分处理。
• 参数与迭代逻辑检查：
  • 确认Bhat的维度为(num_vars*num_lags, num_vars)，确保lvars.flatten().T @ Bhat的矩阵乘法维度匹配，避免因维度错误导致计算异常。
  • 验证A0[:,2]是否对应正确的结构冲击（Gertler-Karadi中通常是货币政策冲击，即短期利率冲击），若初始冲击方向或幅度错误，也会导致后续响应发散。
  • 检查迭代时滞后项提取逻辑：np.arange(start=num_lags+jj-1, stop=jj-1, step=-1)是否准确提取了前num_lags期的脉冲响应值，避免漏取或多取滞后项。

内容的提问来源于stack exchange，提问作者Mandar Priya Phatak