计量经济学(九)——向量自回归VAR模型检验
向量自回归（VAR，Vector Autoregression）模型是一种广泛用于时间序列分析的统计工具，特别是在经济学和金融学领域中。VAR模型的关键优势在于其可以捕捉多个变量之间的相互依赖关系，而无需预设变量之间的因果顺序。这使得VAR模型在处理复杂动态系统时极具灵活性。VAR模型的基本结构是将多个时间序列变量构建成一个多元系统，假设每个变量不仅受其自身的历史值影响，还受到其他变量的滞后项影响。实际应用中，VAR模型可以用于政策冲击分析、预测经济指标、并通过冲击响应函数衡量变量间的短期和长期影响，广泛应用于经济预测、货币政策分析等领域。
一、向量自回归VAR模型
向量自回归模型（VAR, Vector Autoregression）是一种多变量时间序列模型，它假设每个变量不仅受自身过去的滞后影响，还受其他变量的滞后影响。VAR模型的核心思想是将多个变量的自回归模型结合起来，从而构建一个包含多个方程的系统。在这个系统中，每个方程都以一个变量为因变量，其他变量作为自变量出现在模型中。
1.1 单变量自回归模型
在讨论VAR模型之前，我们先回顾单变量自回归（AR，Autoregressive）模型。对于一个时间序列 \(y_t\)，AR模型可以表示为：
其中，$ \alpha $ 是常数项，$ \beta_i $ 是滞后系数，$ p $ 是滞后阶数，$ \epsilon_t $ 是误差项，服从均值为0、方差为 $ \sigma^2 $ 的白噪声过程。AR模型的基本思想是，当前值 $ y_t $ 可以由过去的 $ p $ 个值来解释。
1.2 向量自回归模型
在VAR模型中，多个时间序列变量相互影响。设我们有 $ K $ 个时间序列变量 $ y_{1t}, y_{2t}, \ldots, y_{Kt} $，这些变量构成一个向量 $ Y_t $，即：
VAR(p) 模型的基本形式如下：
其中：
- $ Y_t $ 是 $ K \times 1 $ 的向量，表示时间 $ t $ 上的 $ K $ 个时间序列变量。
- $ A_0 $ 是 $ K \times 1 $ 的常数项向量。
- $ A_1, A_2, \ldots, A_p $ 是 $ K \times K $ 的系数矩阵。
- $ \epsilon_t $ 是 $ K \times 1 $ 的误差向量，假设 $ \epsilon_t $ 服从均值为0、协方差矩阵为 $ \Sigma $ 的正态分布，即 $ \epsilon_t \sim N(0, \Sigma) $。
VAR模型的每一个变量都由自身的滞后项以及其他变量的滞后项来解释，这就使得VAR模型能够捕捉变量之间复杂的动态关系。
二、VAR模型的估计
VAR模型的参数估计通常通过最小二乘法（OLS, Ordinary Least Squares）来实现。由于VAR模型是一个线性系统，因此可以分别对每个方程使用OLS进行估计。假设我们有 $ T $ 期数据，估计步骤如下：
2.1 数据矩阵表示
首先，将VAR模型表示为矩阵形式。设 $ Y $ 是 $ T \times K $ 的观测矩阵，$ X $ 是 $ T \times (Kp + 1) $ 的自变量矩阵，模型可以表示为：
其中：
- $ Y $ 是包含所有时间序列观测值的矩阵，每一列对应一个变量，每一行为一个时间点的观测值。
- $ X $ 是包含所有滞后项和常数项的自变量矩阵。
- $ B $ 是 $ (Kp + 1) \times K $ 的系数矩阵，每列对应一个方程的系数。
- $ E $ 是 $ T \times K $ 的残差矩阵。
2.2 最小二乘估计
对每个方程进行OLS估计，即对 $ Y $ 中的每一列进行回归，得到系数矩阵 $ B $ 的估计值 $ \hat{B} $。最小二乘估计公式为：
估计得到的 $ \hat{B} $ 是每个方程中滞后项的系数，代表了每个变量对其他变量滞后值的影响。
2.3 VAR模型的阶数选择
在实际应用中，选择VAR模型的滞后阶数 $ p $ 是一个关键问题。如果滞后阶数过小，模型可能不能充分捕捉变量之间的动态关系；如果滞后阶数过大，则会导致模型过拟合并降低预测精度。为了确定最佳的滞后阶数，常用的信息准则有AIC（Akaike Information Criterion）和BIC（Bayesian Information Criterion）。
AIC准则的定义为：
其中，$ \hat{\Sigma} $ 是残差协方差矩阵的估计值，$ T $ 是样本的观测期数，$ K $ 是时间序列的维度。
BIC准则与AIC类似，但对模型的复杂度惩罚更为严格。BIC的定义为：
AIC和BIC准则的计算结果随着滞后阶数 $ p $ 的变化而变化，通常选择使AIC或BIC取最小值的滞后阶数作为模型的最佳滞后阶数。
2.4 VAR模型的平稳性分析
VAR模型的平稳性是一个重要的假设。只有在模型是平稳的情况下，VAR模型的估计结果才是可靠的。平稳性意味着时间序列的均值、方差和协方差在时间上不随时间变化。对于VAR模型来说，平稳性可以通过检查特征根来判断。
一个VAR(p) 模型是平稳的，当且仅当模型的特征根的模小于1。具体而言，首先计算VAR模型对应的矩阵 $ A_1, A_2, \ldots, A_p $ 的特征值，然后判断这些特征值的模是否小于1。如果所有特征根的模都小于1，则VAR模型是平稳的。
2.5 VAR模型的冲击响应函数
在VAR模型中，不同变量之间的动态关系通过各自的滞后项表现出来。当系统中的一个变量受到某种冲击时（例如，经济中的某项政策变化，或金融市场的一个突然波动），其他相关变量可能会受到影响。简单来说，冲击响应函数回答了这样一个问题：如果系统中的一个变量突然发生变化，那么这个变化将如何影响其他变量，并且这种影响会持续多久？
冲击响应函数是指在给定初始冲击的情况下，系统中其他变量对这一冲击的响应。设 $ \epsilon_t $ 是时间 $ t $ 的误差项，冲击响应函数描述了一个单位冲击对未来时期变量的影响。具体来说，假设在时间 $ t $ 上对变量 $ y_{1t} $ 施加一个单位冲击，则冲击响应函数 $ \psi_{ij}(h) $ 描述了第 $ j $ 个变量在 $ t+h $ 时刻对这个冲击的响应：
冲击响应函数通常通过对VAR模型进行迭代计算得到。具体步骤是：
- 估计VAR模型的系数矩阵 $ A_1, A_2, \ldots, A_p $。
- 计算误差项协方差矩阵 $ \Sigma $，并对其进行Cholesky分解。
- 使用迭代方法计算不同滞后期的冲击响应函数。
冲击响应函数的目的是在给定初始冲击后，研究在未来的每个时间点\(t+h\)上，系统中各个变量的响应。通常，我们用图形展示各个变量对特定冲击的响应情况，称为冲击响应图。冲击响应函数通过研究一个特定的外生冲击（即误差项\(\epsilon_t\)的变化）对所有内生变量（即\(y_t\)中的所有变量）的影响，揭示了系统内部的动态联系。影响路径：冲击如何影响系统中不同的变量？影响大小：某个冲击对某个变量的影响有多大？影响持续时间：该冲击的影响会持续多长时间？这些信息对于政策评估和金融市场的分析都至关重要。冲击响应函数可以帮助分析经济变量之间的动态传递机制。例如，在宏观经济学中，冲击响应函数常用于研究货币政策冲击如何影响GDP、通货膨胀率和利率等宏观经济变量。
2.6 VAR模型的预测
VAR模型可以用于多个时间序列变量的联合预测。预测过程基于已估计的VAR模型，通过逐步外推未来的滞后值来进行。
一步预测：一步预测是指在时间 $ t $ 预测时间 $ t+1 $ 的值。根据已估计的VAR模型的系数矩阵 $ A_1, A_2, \ldots, A_p $，我们可以计算出未来的预测值：
多步预测：多步预测是指在时间 $ t $ 预测时间 $ t+h $ 的值。多步预测通过将一步预测的结果作为滞后值继续迭代，直到达到目标期：
VAR模型的预测性能在经济学、金融市场等领域有广泛的应用。例如，VAR模型可以用于预测股票市场的价格走势、利率的变化，以及货币政策对经济的影响。
2.7 VAR模型的局限性与改进
虽然VAR模型在时间序列分析中具有广泛应用，但它也存在一些局限性。主要的局限性包括：
- 滞后阶数选择难题：如果滞后阶数选择过多，会导致模型复杂性增加，参数估计不稳定；如果滞后阶数选择过少，可能无法充分捕捉变量之间的动态关系。
- 平稳性要求：VAR模型要求时间序列是平稳的，但在实际应用中，很多经济时间序列往往是非平稳的。在这种情况下，需要通过差分等方式将非平稳序列转化为平稳序列。
- 解释性较弱：VAR模型的解释性相对较弱，它并不区分变量之间的因果关系，而只是捕捉变量之间的相互依赖关系。为了解决这一问题，研究者常常结合因果检验（如Granger因果检验）或结构VAR（SVAR）模型进行分析。
结构VAR模型
结构VAR（SVAR）模型是一种改进的VAR模型，它通过引入额外的约束条件来识别因果关系。SVAR模型的基本思想是，经济变量之间的冲击可能有着不同的来源，SVAR模型试图通过识别这些结构性冲击来进行更为精确的分析。
SVAR模型的形式为：
其中，$ A $ 和 $ B $ 是结构系数矩阵，$ \epsilon_t $ 是结构性冲击。SVAR模型的估计通常需要结合经济理论进行约束，通过对VAR模型中的误差项进行转化，得到更加符合实际经济关系的冲击响应分析结果。
VECM模型
对于非平稳的时间序列，向量误差修正模型（VECM，Vector Error Correction Model）是VAR模型的一个重要改进。VECM模型通过引入误差修正项，能够捕捉非平稳序列之间的长期均衡关系。VECM模型的基本形式为：
其中，$ \Delta Y_t $ 表示差分后的时间序列，$ \beta' Y_{t-1} $ 是误差修正项，表示变量之间的长期均衡关系。VECM模型特别适合于处理存在协整关系的非平稳时间序列。
三、VAR检验步骤
VAR模型有着较为复杂的分析流程步骤，具体可参考下图和表格所述。下述分析步骤和流程仅供参考使用，实际研究中可能仅需要其中一部分的分析即可。比如很多时候并不需要关注残差自相关检验和残差正态性检验，也或者有时对格兰杰因果检验关注度较少等，具体以研究者实际研究情况为准即可。
| 项 | 说明 | 项 | 说明 |
|---|---|---|---|
| 单位根检验 | 模型要求研究变量保持平稳性，可使用单位根检验检验平稳性性，也或者使用时序图亦可。 | 脉冲响应 | 其用于研究变量间的动态影响关系情况。 |
| 协整关系 | 协整关系研究系统的长期稳定性情况，以避免出现伪回归问题，通常有两种方法，分别是Johansen法和Engle-Granger法研究协整关系。 | 方差分解 | 其用于研究变量间的影响比例情况。 |
| 定阶 | 进行VAR模型时，需要针对滞后阶数进行选择，SPSSAU可自动识别最佳滞后阶数。 | 残差自相关检验 | 研究模型残差的自相关性情况。 |
| AR特征根检验 | AR特征根检验可用于分析VAR模型系数稳定有效性情况。 | 残差正态性检验 | 研究模型残差的正态性情况。 |
| 格兰杰因果 | 研究变量间的‘相互影响关系’情况。 | 预测 | 得到模型的预测值，用于研究使用。 |
VAR模型的构建涉及多个重要步骤，从单位根检验、定阶到格兰杰因果关系和脉冲响应函数，每一个步骤都对模型的准确性至关重要。
单位根检验：在建模之前，首先需要进行单位根检验（如ADF检验）以判断时间序列的平稳性。如果数据平稳（即无单位根），可以直接进入VAR模型的构建。如果数据具有单位根（不平稳），则需要对其进行差分处理，使其成为平稳序列。如果变量经过同阶单整后可以满足平稳条件，则继续进行协整关系分析，否则无法进行有效的VAR建模。
协整关系分析：当时间序列数据经过差分后，若发现变量之间存在长期均衡关系（协整），则可以结合协整检验，如Johansen检验，来确认这些变量之间的协整向量。若无单位根或数据极为平稳，则直接进行VAR模型；若存在协整关系，则应考虑使用VECM（向量误差修正模型）来处理。
滞后阶数选择：确定VAR模型的滞后阶数是关键的一步，可以通过信息准则如AIC、BIC、HQ等方法来确定。滞后阶数影响模型的动态特性，准则值越小表示模型的解释能力越好。SPSS等工具通常可以自动给出建议的滞后阶数。
AR特征根检验：模型构建后，必须对其稳定性进行检验，通常通过AR特征根检验来验证。如果特征根的模数都在单位圆内，说明模型是稳定的，参数也具有可靠性。
格兰杰因果检验：构建VAR模型后，可以进行格兰杰因果检验，来判断某个变量的过去值是否能够预测另一个变量的未来变化。它可以揭示变量间的动态相关性，帮助理解变量间的因果链。
脉冲响应函数：通过正交化脉冲响应分析，可以研究一个变量的冲击如何随着时间影响其他变量。这个分析对理解经济冲击的传播路径和力度十分重要。
方差分解：在VAR模型的后期分析中，方差分解用于了解每个变量的方差有多大比例是由自己和其他变量解释的。这有助于量化不同变量之间的相互影响。
残差检验：构建模型后，还需要对残差进行自相关和正态性检验，以评估模型的拟合效果。Portmanteau检验可以判断残差是否存在自相关，而Jarque-Bera检验可以检测残差的正态性。
预测分析：VAR模型的最终目标通常是进行预测。通过模型可以对未来的多个时间点进行预测，常见的预测期数为12期。
通过以上步骤，研究者能够构建出具有较强解释力的VAR模型，并根据模型的输出进行深入的经济或金融领域的分析与预测。
四、案例分析
这里将国内钢铁现货价格（普钢价格指数pg）、期货价格（螺纹钢主力合约收盘价lwg）和中国铁矿石指数(tks)纳入到动态的 VAR 系统中，以考察钢铁价格的传导，计算各变量的贡献度与影响效应。
4.1 数据的选取
所有变量选取自2015年2月25至2021年9月30日，样本容量为1649。数据源自前瞻数据库（https://x.qianzhan.com），有的直接就选用了指数，有的进行了定基化处理，为消除异方差数据取对数来进行分析。全部数据蓝奏云链接下载
| 时间序号 | TKS | PG | LWG | 时间序号 | TKS | PG | LWG |
|---|---|---|---|---|---|---|---|
| 1 | 228.95 | 94.64 | 2483 | 5 | 230.35 | 94.4 | 2472 |
| 2 | 231.08 | 94.67 | 2505 | 6 | 231.19 | 94.18 | 2573 |
| 3 | 230.2 | 94.67 | 2464 | 7 | 229.66 | 94.03 | 2565 |
| 4 | 205.54 | 94.62 | 2468 | 8 | 224.69 | 93.84 | 2523 |
| 时序图 | 自相关图 | 偏相关图 |
|---|---|---|
 |
 |
 |
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm
# 读取当前路径下的data1.csv文件
data = pd.read_csv('data1.csv')
print(data)
# 提取出pg列，并将其转换为时间序列
pg_series = pd.Series(data['PG'])
# 画出pg列的时间序列图
plt.figure(figsize=(10, 6))
plt.plot(pg_series, label='PG Time Series')
plt.title('PG Column Time Series')
plt.xlabel('Time')
plt.ylabel('PG')
plt.legend()
plt.grid(True)
plt.show()
# 单独绘制自相关图（ACF）
plt.figure(figsize=(10, 6))
sm.graphics.tsa.plot_acf(pg_series, lags=40)
plt.title('ACF of PG Column')
plt.grid(True)
plt.show()
# 单独绘制偏自相关图（PACF）
plt.figure(figsize=(10, 6))
sm.graphics.tsa.plot_pacf(pg_series, lags=40)
plt.title('PACF of PG Column')
plt.grid(True)
plt.show()
4.2 平稳性检验
鉴于时间序列分析和回归分析有许多前提假定，如序列的平稳性、正态性等，为规避出现伪回归现象，先对序列进行平稳性检验。一个平稳过程的数据图形特征为：数据围绕长期均值波动，偏离均值之后有复归均值的调整；方差variance有限且不随时间改变；其自相关函数随时间衰减。非平稳过程其数据图形特征为：不存在长期均值；方差具有时变性且趋于无穷；自相关不随时间衰减。
import pandas as pd
import numpy as np
import statsmodels.api as sm
from statsmodels.tsa.stattools import adfuller
# 读取当前路径下的data1.csv文件
data = pd.read_csv('data1.csv')
# 提取出 PG, LWG, TKS 列数据
series_pg = data['PG']
series_lwg = data['LWG']
series_tks = data['TKS']
# ADF检验函数
def adf_test(series, series_name):
print(f'ADF Test for {series_name}:')
result = adfuller(series, autolag='AIC')
print(f'ADF Statistic: {result[0]}')
print(f'p-value: {result[1]}')
print('Critical Values:')
for key, value in result[4].items():
print(f'\t{key}: {value}')
if result[1] < 0.05:
print(f"The series {series_name} is stationary.\n")
else:
print(f"The series {series_name} is not stationary. Differencing is required.\n")
print('-' * 50)
# 检验PG, LWG, TKS原序列的平稳性
adf_test(series_pg, 'PG')
adf_test(series_lwg, 'LWG')
adf_test(series_tks, 'TKS')
# 对不平稳的序列进行一阶差分
diff_pg = series_pg.diff().dropna()
diff_lwg = series_lwg.diff().dropna()
diff_tks = series_tks.diff().dropna()
# 检验差分后的平稳性
print("\n差分后的序列平稳性检验：")
adf_test(diff_pg, 'Differenced PG')
adf_test(diff_lwg, 'Differenced LWG')
adf_test(diff_tks, 'Differenced TKS')
# 生成类似Eviews输出表格的结果
def generate_adf_output(series, diff_series, series_name):
original_test = adfuller(series, autolag='AIC')
diff_test = adfuller(diff_series, autolag='AIC')
# 打印类似Eviews输出的结果表格
print(f"{'Test':<20}{'ADF Statistic':<20}{'p-value':<20}{'Critical Values'}")
print(f"{'-' * 80}")
print(f"{series_name:<20}{original_test[0]:<20.4f}{original_test[1]:<20.4f}{original_test[4]}")
print(f"Differenced {series_name:<12}{diff_test[0]:<20.4f}{diff_test[1]:<20.4f}{diff_test[4]}")
print(f"{'-' * 80}")
# 生成输出表格
generate_adf_output(series_pg, diff_pg, 'PG')
generate_adf_output(series_lwg, diff_lwg, 'LWG')
generate_adf_output(series_tks, diff_tks, 'TKS')
--------------------------------------------------
Test ADF Statistic p-value Critical Values
--------------------------------------------------------------------------------
PG -0.0773 0.9517 {'1%': -3.434363244707632, '5%': -2.8633126160514863, '10%': -2.567713701716287}
Differenced PG -10.0992 0.0000 {'1%': -3.434363244707632, '5%': -2.8633126160514863, '10%': -2.567713701716287}
--------------------------------------------------------------------------------
Test ADF Statistic p-value Critical Values
--------------------------------------------------------------------------------
LWG -0.2741 0.9291 {'1%': -3.4343387651403106, '5%': -2.8633018108486343, '10%': -2.5677079480625067}
Differenced LWG -17.1263 0.0000 {'1%': -3.4343387651403106, '5%': -2.8633018108486343, '10%': -2.5677079480625067}
--------------------------------------------------------------------------------
Test ADF Statistic p-value Critical Values
--------------------------------------------------------------------------------
TKS -1.4342 0.5658 {'1%': -3.434363244707632, '5%': -2.8633126160514863, '10%': -2.567713701716287}
Differenced TKS -8.5975 0.0000 {'1%': -3.434363244707632, '5%': -2.8633126160514863, '10%': -2.567713701716287}
--------------------------------------------------------------------------------
由上表可得, \(PG，LWG，TKS\) 对于任何显著性水平都接受原假设, 为非平稳时间序列；但其一阶差分序列均拒绝原假设，各变量均是一阶单整的时间序列。
4.3 VAR模型滞后项选取
import pandas as pd
import numpy as np
import statsmodels.api as sm
from statsmodels.tsa.api import VAR
# 读取当前路径下的data1.csv文件
data = pd.read_csv('data1.csv')
# 提取 PG, LWG, TKS 列数据
data_var = data[['PG', 'LWG', 'TKS']]
# 创建VAR模型
model = VAR(data_var)
# 使用信息准则确定最佳滞后阶数
lag_selection = model.select_order(maxlags=15) # 假设最多滞后15阶
print(lag_selection.summary())
# 打印各个信息准则的选择结果
print("AIC-selected lag: ", lag_selection.aic)
print("BIC-selected lag: ", lag_selection.bic)
print("HQIC-selected lag: ", lag_selection.hqic)
# 根据AIC准则选择最佳滞后阶数进行模型拟合
optimal_lag = lag_selection.aic
var_result = model.fit(optimal_lag)
# 打印VAR模型结果
print(var_result.summary())
VAR Order Selection (* highlights the minimums)
==================================================
AIC BIC FPE HQIC
--------------------------------------------------
0 25.83 25.84 1.649e+11 25.83
1 11.75 11.79 1.270e+05 11.77
2 11.42 11.49* 9.106e+04 11.45
3 11.42 11.52 9.084e+04 11.45
4 11.39 11.52 8.837e+04 11.44*
5 11.39 11.55 8.814e+04 11.45
6 11.38 11.57 8.744e+04 11.45
7 11.37 11.59 8.704e+04 11.46
8 11.38 11.63 8.758e+04 11.47
9 11.38 11.65 8.726e+04 11.48
10 11.36 11.67 8.614e+04 11.48
11 11.37 11.70 8.633e+04 11.49
12 11.35* 11.72 8.525e+04* 11.49
13 11.36 11.76 8.593e+04 11.51
14 11.36 11.79 8.613e+04 11.52
15 11.37 11.82 8.637e+04 11.54
--------------------------------------------------
AIC-selected lag: 12
BIC-selected lag: 2
HQIC-selected lag: 4
这里选择滞后项为4。
4.4 VAR分析
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.api import VAR
from statsmodels.stats.stattools import durbin_watson
# 读取数据
data = pd.read_csv('data1.csv')
# 选择 PG, LWG, TKS 作为 VAR 模型的变量
data_var = data[['PG', 'LWG', 'TKS']]
# 创建 VAR 模型
model = VAR(data_var)
# 拟合 VAR 模型，滞后项选择为 4
var_result = model.fit(4)
# 输出 VAR 模型的结果
print(var_result.summary())
# Durbin-Watson 检验
dw = durbin_watson(var_result.resid)
for col, val in zip(data_var.columns, dw):
print(f'Durbin-Watson statistic for {col}: {val}')
# ------------------ 特征根稳定性检验 ------------------
print("\n特征根检验：")
is_stable = var_result.is_stable()
print(f"模型是否稳定: {is_stable}")
# 如果模型稳定，特征根的模应当全部小于 1
print("特征根：")
print(var_result.roots)
# 绘制特征根图
plt.figure(figsize=(6, 6))
plt.scatter(var_result.roots.real, var_result.roots.imag, marker='o', label='Roots')
plt.title('Roots of the VAR Model')
plt.xlabel('Real Part')
plt.ylabel('Imaginary Part')
plt.axvline(x=0, color='k', linestyle='--')
plt.axhline(y=0, color='k', linestyle='--')
plt.grid(True)
plt.legend()
plt.show()
# ------------------ 脉冲响应分析 ------------------
irf = var_result.irf(10) # 设置脉冲响应分析的预测期为10
irf.plot(orth=True)
plt.title("Impulse Response Functions (Orthogonalized)")
plt.show()
# ------------------ 方差分解 ------------------
fevd = var_result.fevd(10) # 进行方差分解，预测 10 期
fevd.plot(figsize=(10, 6))
plt.title("Variance Decomposition")
plt.show()
# ------------------ 预测未来10期 ------------------
forecast_data = var_result.forecast(data_var.values[-4:], steps=10) # 预测未来 10 期
forecast_index = pd.date_range(start=len(data_var), periods=10, freq='Q') # 时间序列索引
# 将预测结果转换为DataFrame
forecast_df = pd.DataFrame(forecast_data, index=forecast_index, columns=['PG', 'LWG', 'TKS'])
print("\n未来10期的预测结果：")
print(forecast_df)
# 绘制预测结果
plt.figure(figsize=(10, 6))
for i, col in enumerate(forecast_df.columns):
plt.subplot(3, 1, i + 1)
plt.plot(data_var.index, data_var[col], label=f'Actual {col}')
plt.plot(forecast_index, forecast_df[col], label=f'Forecasted {col}')
plt.title(f'{col} Forecast for the next 10 periods')
plt.legend(loc='best')
plt.grid(True)
plt.tight_layout()
plt.show()
最佳滞后阶数：4
未来10期的预测结果：
PG LWG TKS
1970-03-31 00:00:00.000001649 214.651002 5701.425583 440.042314
1970-06-30 00:00:00.000001649 215.180854 5699.876328 439.969412
1970-09-30 00:00:00.000001649 215.557975 5697.933981 440.061312
1970-12-31 00:00:00.000001649 215.914655 5696.978733 439.291193
1971-03-31 00:00:00.000001649 216.184004 5694.033977 438.596610
1971-06-30 00:00:00.000001649 216.393839 5690.734004 437.941639
1971-09-30 00:00:00.000001649 216.583981 5687.305953 437.242164
1971-12-31 00:00:00.000001649 216.749471 5683.843872 436.547578
1972-03-31 00:00:00.000001649 216.895330 5680.375887 435.877311
1972-06-30 00:00:00.000001649 217.029351 5676.965992 435.215078
4.5 检验通过后的VAR模型
VAR Model Expressions
Equation for PG
Equation for LWG
Equation for TKS
4.6 结果分析
（1）从长期来看，普钢价格指数、螺纹钢期货价格、铁矿石价格指数具有长期稳定的协整关系。钢铁价格的变动受到铁矿石价格和螺纹钢期货价格的影响，普钢价格的变动比铁矿石价格波动要小。
（2）从脉冲响应函数来看，铁矿石指数对普钢指数的影响时滞为4期。中国铁矿石价格指数由国内和国外铁矿石指数两部分构成，削弱了国外铁矿石指数对钢铁由于价格变化的影响。钢铁价格复杂的传导机制，对于钢铁产业转型升级带来了不小挑战，要通过多种方式沿着产业链来传导价格，避免过度涨价和跌价，保持钢铁产业价格温和变动。
（3）从方差分解结果来看，普钢指数的波动大部分受自己滞后项影响，其他变量的贡献比例不尽相同。计算结果表明铁矿石价格对钢铁价格的波动产生的影响最大，其后是期货价格。说明在钢铁这种大宗产品价格传导机制中，核心原料必须掌控，强化铁矿石定价权的争夺。
（4）着力发展期货市场。铁矿石市场定价权主导的是发达国家，我国处于价格被动接受者地位。期货市场能扩大钢铁产业上下游期货市场的交易规模，提高钢铁市场国际化程度，吸引更多国内外投资者的广泛参与。期货市场的开发，是提升我国对钢铁产业这类大宗商品的影响力与争取定价权的主要途径。
（5）完善钢铁产业发展联盟，建立钢铁产业采购基金。出台钢铁产业的宏观调控政策，以国家信誉背书，统一管理采购进口铁矿石。国家利用储备货币进行铁矿石期货交易，增强铁矿石采购的主动权、话语权、决策权；可将货币政策金融工具纳入政策体系内，提高外汇的使用率，共同应对铁矿石市场的风险。提升产业的智能化、数字化发展，精准监控整个产业的发展动势，催生产业数字化，推动数字经济和钢铁产业深度融合，打造具有国际竞争力的数字钢铁产业集群。
总结
VAR模型（向量自回归模型）作为一种多变量时间序列分析方法，因其灵活性和应用广泛性在经济学和金融学中得到了广泛使用。VAR模型通过捕捉多个时间序列变量之间的相互依赖关系，能够在不预设因果关系的情况下，分析这些变量之间的动态变化。它被用于预测经济变量、评估政策影响以及进行冲击响应分析。其优势在于能够综合多个变量的历史信息，做出更精确的预测和分析。
尽管VAR模型具有强大的分析能力，但也存在一些局限性。例如，模型的滞后阶数选择可能影响结果的准确性，同时VAR模型无法有效处理非平稳数据。不过，研究者可以通过改进的模型形式，如结构VAR（SVAR）和误差修正模型（VECM），解决这些问题并提高模型的解释能力。在实际应用中，研究者应根据数据特征和分析需求，选择合适的VAR模型，并结合经济理论、实证研究和经验分析，对结果进行合理解读，从而做出更精确的预测和有效的政策建议。
