向量自回归 tsa.vector_ar

statsmodels.tsa.vector_ar 包含了用于使用 向量自回归 (VAR)向量误差修正模型 (VECM) 同时建模和分析多个时间序列的方法。

VAR(p) 过程

我们感兴趣的是对一个 \(T \times K\) 多元时间序列 \(Y\) 进行建模,其中 \(T\) 表示观测值数量,\(K\) 表示变量数量。估计时间序列与其滞后值之间关系的一种方法是向量自回归过程

\[ \begin{align}\begin{aligned}Y_t = \nu + A_1 Y_{t-1} + \ldots + A_p Y_{t-p} + u_t\\u_t \sim {\sf Normal}(0, \Sigma_u)\end{aligned}\end{align} \]

其中 \(A_i\) 是一个 \(K \times K\) 系数矩阵。

我们在很大程度上遵循 Lutkepohl (2005) 的方法和符号,我们这里不进行展开。

模型拟合

注意

下面引用的类可以通过 statsmodels.tsa.api 模块访问。

为了估计 VAR 模型,必须首先使用一个具有同质或结构化数据类型的 ndarray 来创建模型。当使用结构化或记录数组时,该类将使用传入的变量名。否则,它们可以显式传入

# some example data
In [1]: import numpy as np

In [2]: import pandas

In [3]: import statsmodels.api as sm

In [4]: from statsmodels.tsa.api import VAR

In [5]: mdata = sm.datasets.macrodata.load_pandas().data

# prepare the dates index
In [6]: dates = mdata[['year', 'quarter']].astype(int).astype(str)

In [7]: quarterly = dates["year"] + "Q" + dates["quarter"]

In [8]: from statsmodels.tsa.base.datetools import dates_from_str

In [9]: quarterly = dates_from_str(quarterly)

In [10]: mdata = mdata[['realgdp','realcons','realinv']]

In [11]: mdata.index = pandas.DatetimeIndex(quarterly)

In [12]: data = np.log(mdata).diff().dropna()

# make a VAR model
In [13]: model = VAR(data)

注意

VAR 类假设传入的时间序列是平稳的。非平稳或趋势数据通常可以通过一阶差分或其他方法转换为平稳数据。对于非平稳时间序列的直接分析,标准稳定的 VAR(p) 模型不适用。

要实际进行估计,请使用所需的滞后阶数调用 fit 方法。或者,您可以让模型根据标准信息准则选择滞后阶数(参见下文)。

In [14]: results = model.fit(2)

In [15]: results.summary()
Out[15]: 
  Summary of Regression Results   
==================================
Model:                         VAR
Method:                        OLS
Date:           Thu, 03, Oct, 2024
Time:                     16:15:42
--------------------------------------------------------------------
No. of Equations:         3.00000    BIC:                   -27.5830
Nobs:                     200.000    HQIC:                  -27.7892
Log likelihood:           1962.57    FPE:                7.42129e-13
AIC:                     -27.9293    Det(Omega_mle):     6.69358e-13
--------------------------------------------------------------------
Results for equation realgdp
==============================================================================
                 coefficient       std. error           t-stat            prob
------------------------------------------------------------------------------
const               0.001527         0.001119            1.365           0.172
L1.realgdp         -0.279435         0.169663           -1.647           0.100
L1.realcons         0.675016         0.131285            5.142           0.000
L1.realinv          0.033219         0.026194            1.268           0.205
L2.realgdp          0.008221         0.173522            0.047           0.962
L2.realcons         0.290458         0.145904            1.991           0.047
L2.realinv         -0.007321         0.025786           -0.284           0.776
==============================================================================

Results for equation realcons
==============================================================================
                 coefficient       std. error           t-stat            prob
------------------------------------------------------------------------------
const               0.005460         0.000969            5.634           0.000
L1.realgdp         -0.100468         0.146924           -0.684           0.494
L1.realcons         0.268640         0.113690            2.363           0.018
L1.realinv          0.025739         0.022683            1.135           0.257
L2.realgdp         -0.123174         0.150267           -0.820           0.412
L2.realcons         0.232499         0.126350            1.840           0.066
L2.realinv          0.023504         0.022330            1.053           0.293
==============================================================================

Results for equation realinv
==============================================================================
                 coefficient       std. error           t-stat            prob
------------------------------------------------------------------------------
const              -0.023903         0.005863           -4.077           0.000
L1.realgdp         -1.970974         0.888892           -2.217           0.027
L1.realcons         4.414162         0.687825            6.418           0.000
L1.realinv          0.225479         0.137234            1.643           0.100
L2.realgdp          0.380786         0.909114            0.419           0.675
L2.realcons         0.800281         0.764416            1.047           0.295
L2.realinv         -0.124079         0.135098           -0.918           0.358
==============================================================================

Correlation matrix of residuals
             realgdp  realcons   realinv
realgdp     1.000000  0.603316  0.750722
realcons    0.603316  1.000000  0.131951
realinv     0.750722  0.131951  1.000000

有多种方法可以使用 matplotlib 可视化数据。

绘制输入时间序列

In [16]: results.plot()
Out[16]: <Figure size 1000x1000 with 3 Axes>
_images/var_plot_input.png

绘制时间序列自相关函数

In [17]: results.plot_acorr()
Out[17]: <Figure size 1000x1000 with 9 Axes>
_images/var_plot_acorr.png

滞后阶数选择

滞后阶数的选择可能是一个难题。标准分析采用似然比检验或基于信息准则的阶数选择。我们已经实现了后者,可以通过 VAR 类访问。

In [18]: model.select_order(15)
Out[18]: <statsmodels.tsa.vector_ar.var_model.LagOrderResults at 0x7f18fd4e8430>

调用 fit 函数时,可以传入最大滞后数量和用于阶数选择的准则。

In [19]: results = model.fit(maxlags=15, ic='aic')

预测

线性预测器是在均方误差方面最优的 h 步超前预测

\[y_t(h) = \nu + A_1 y_t(h − 1) + \cdots + A_p y_t(h − p)\]

我们可以使用 forecast 函数生成该预测。请注意,我们必须指定预测的“初始值”。

In [20]: lag_order = results.k_ar

In [21]: results.forecast(data.values[-lag_order:], 5)
Out[21]: 
array([[ 0.0062,  0.005 ,  0.0092],
       [ 0.0043,  0.0034, -0.0024],
       [ 0.0042,  0.0071, -0.0119],
       [ 0.0056,  0.0064,  0.0015],
       [ 0.0063,  0.0067,  0.0038]])

forecast_interval 函数将生成上述预测以及渐近标准误。可以使用 plot_forecast 函数可视化这些结果。

In [22]: results.plot_forecast(10)
Out[22]: <Figure size 1000x1000 with 3 Axes>
_images/var_forecast.png

类参考

VAR(endog[, exog, dates, freq, missing])

拟合 VAR(p) 过程并进行滞后阶数选择

VARProcess(coefs, coefs_exog, sigma_u[, ...])

类表示一个已知的 VAR(p) 过程

VARResults(endog, endog_lagged, params, ...)

估计具有固定滞后数量的 VAR(p) 过程

估计后分析

对于向量自回归过程,估计后可以获得多个过程属性和额外结果。

LagOrderResults(ics, selected_orders[, vecm])

用于选择模型滞后阶数的结果类。

HypothesisTestResults(test_statistic, ...)

用于假设检验的结果类。

NormalityTestResults(test_statistic, ...)

用于 Jarque-Bera 非正态性检验的结果类。

WhitenessTestResults(test_statistic, ...)

用于残差自相关检验的 Portmanteau 结果类。

脉冲响应分析

脉冲响应在计量经济学研究中很重要:它们是对一个变量的单位脉冲的估计响应。在实践中,它们是使用 VAR(p) 过程的 MA(\(\infty\)) 表示来计算的

\[Y_t = \mu + \sum_{i=0}^\infty \Phi_i u_{t-i}\]

我们可以通过在 VARResults 对象上调用 irf 函数来执行脉冲响应分析。

In [23]: irf = results.irf(10)

这些结果可以使用 plot 函数可视化,以正交化或非正交化形式呈现。默认情况下,会在 95% 的显著性水平下绘制渐近标准误,用户可以修改该水平。

注意

正交化是使用估计误差协方差矩阵 \(\hat \Sigma_u\) 的 Cholesky 分解来完成的,因此解释可能会因变量排序的不同而改变。

In [24]: irf.plot(orth=False)
Out[24]: <Figure size 1000x1000 with 9 Axes>
_images/var_irf.png

请注意,plot 函数很灵活,如果需要,可以只绘制感兴趣的变量。

In [25]: irf.plot(impulse='realgdp')
Out[25]: <Figure size 1000x1000 with 3 Axes>
_images/var_realgdp.png

累积效应 \(\Psi_n = \sum_{i=0}^n \Phi_i\) 可以与长期效应一起绘制,如下所示。

In [26]: irf.plot_cum_effects(orth=False)
Out[26]: <Figure size 1000x1000 with 9 Axes>
_images/var_irf_cum.png

IRAnalysis(model[, P, periods, order, svar, ...])

脉冲响应分析类。

预测误差方差分解 (FEVD)

在 i 步预测中,分量 j 对 k 的预测误差可以使用正交脉冲响应 \(\Theta_i\) 进行分解。

\[ \begin{align}\begin{aligned}\omega_{jk, i} = \sum_{i=0}^{h-1} (e_j^\prime \Theta_i e_k)^2 / \mathrm{MSE}_j(h)\\\mathrm{MSE}_j(h) = \sum_{i=0}^{h-1} e_j^\prime \Phi_i \Sigma_u \Phi_i^\prime e_j\end{aligned}\end{align} \]

这些值可以通过 fevd 函数计算,直到预先设定的最大步数。

In [27]: fevd = results.fevd(5)

In [28]: fevd.summary()
FEVD for realgdp
      realgdp  realcons   realinv
0    1.000000  0.000000  0.000000
1    0.864889  0.129253  0.005858
2    0.816725  0.177898  0.005378
3    0.793647  0.197590  0.008763
4    0.777279  0.208127  0.014594

FEVD for realcons
      realgdp  realcons   realinv
0    0.359877  0.640123  0.000000
1    0.358767  0.635420  0.005813
2    0.348044  0.645138  0.006817
3    0.319913  0.653609  0.026478
4    0.317407  0.652180  0.030414

FEVD for realinv
      realgdp  realcons   realinv
0    0.577021  0.152783  0.270196
1    0.488158  0.293622  0.218220
2    0.478727  0.314398  0.206874
3    0.477182  0.315564  0.207254
4    0.466741  0.324135  0.209124

它们也可以通过返回的 FEVD 对象进行可视化。

In [29]: results.fevd(20).plot()
Out[29]: <Figure size 1000x1000 with 3 Axes>
_images/var_fevd.png

FEVD(model[, P, periods])

计算和绘制预测误差方差分解和渐进标准误差。

统计检验

提供了一些不同的方法来进行关于模型结果的假设检验,以及模型假设(正态性、误差的纯白性/“iid-ness”等)的有效性检验。

格兰杰因果关系

人们经常想知道一个变量或一组变量是否对另一个变量“因果”,对于“因果”的一些定义。在 VAR 模型的背景下,可以认为一组变量在其中一个 VAR 方程中是格兰杰因果关系。我们不会详细介绍格兰杰因果关系的数学原理或定义,而是留给读者自己去理解。 VARResults 对象具有 test_causality 方法,用于执行沃尔德 (\(\chi^2\)) 检验或 F 检验。

In [30]: results.test_causality('realgdp', ['realinv', 'realcons'], kind='f')
Out[30]: <statsmodels.tsa.vector_ar.hypothesis_test_results.CausalityTestResults at 0x7f18f3daf580>

正态性

如本文档开头所述,白噪声成分 \(u_t\) 假设服从正态分布。虽然这种假设不是参数估计一致或渐进正态所必需的,但在残差为高斯白噪声时,有限样本中的结果通常更可靠。为了检验这种假设是否与数据集一致,VARResults 提供了 test_normality 方法。

In [31]: results.test_normality()
Out[31]: <statsmodels.tsa.vector_ar.hypothesis_test_results.NormalityTestResults at 0x7f18f3daf640>

残差的纯白性

为了检验估计残差的纯白性(这意味着没有显著的残差自相关性),可以使用 VARResultstest_whiteness 方法。

HypothesisTestResults(test_statistic, ...)

用于假设检验的结果类。

CausalityTestResults(causing, caused, ...[, ...])

格兰杰因果关系和瞬时因果关系的结果类。

NormalityTestResults(test_statistic, ...)

用于 Jarque-Bera 非正态性检验的结果类。

WhitenessTestResults(test_statistic, ...)

用于残差自相关检验的 Portmanteau 结果类。

结构性向量自回归

有一套匹配的类来处理某些类型的结构性 VAR 模型。

SVAR(endog, svar_type[, dates, freq, A, B, ...])

拟合 VAR,然后估计 A 和 B 的结构性成分,定义为

SVARProcess(coefs, intercept, sigma_u, ...)

类表示已知的 SVAR(p) 过程。

SVARResults(endog, endog_lagged, params, ...)

估计具有固定滞后数量的 VAR(p) 过程

向量误差修正模型 (VECM)

向量误差修正模型用于研究一个或多个永久随机趋势(单位根)的短期偏差。VECM 通过对时间序列向量的差值施加结构来对时间序列向量的差值进行建模,这种结构是由假定的随机趋势数量所隐含的。VECM 用于指定和估计这些模型。

VECM(\(k_{ar}-1\)) 具有以下形式

\[\Delta y_t = \Pi y_{t-1} + \Gamma_1 \Delta y_{t-1} + \ldots + \Gamma_{k_{ar}-1} \Delta y_{t-k_{ar}+1} + u_t\]

其中

\[\Pi = \alpha \beta'\]

[1] 的第 7 章所述。

具有确定性项的 VECM(\(k_{ar} - 1\)) 具有以下形式

\[\begin{split}\Delta y_t = \alpha \begin{pmatrix}\beta' & \eta'\end{pmatrix} \begin{pmatrix}y_{t-1} \\ D^{co}_{t-1}\end{pmatrix} + \Gamma_1 \Delta y_{t-1} + \dots + \Gamma_{k_{ar}-1} \Delta y_{t-k_{ar}+1} + C D_t + u_t.\end{split}\]

\(D^{co}_{t-1}\) 中,我们有位于协整关系内部的(或限制在协整关系内的)确定性项。 \(\eta\) 是相应的估计量。为了将确定性项传入协整关系,可以使用 exog_coint 参数。对于截距和线性趋势这两个特殊情况,存在一种更简单的方法来声明这些项:可以分别将 "ci""li" 传递给 deterministic 参数。因此,对于协整关系内的截距,可以将 "ci" 作为 deterministic 传递,或者将 np.ones(len(data)) 作为 exog_coint 传递(如果 data 作为 endog 参数传递)。这确保了 \(D_{t-1}^{co} = 1\) 对于所有 \(t\)

我们还可以使用协整关系之外的确定性项。这些在上述公式中的 \(D_t\) 中定义,相应的估计量位于矩阵 \(C\) 中。通过将这些项传递给 exog 参数来指定这些项。对于截距和/或线性趋势,我们仍然可以使用 deterministic 来作为另一种方法。对于截距,传递 "co",对于线性趋势,传递 "lo",其中 o 代表 outside

下表显示了 [2] 中考虑的五种情况。最后一列指示为每种情况将哪个字符串传递给 deterministic 参数。

情况

截距

线性趋势的斜率

deterministic

I

0

0

"nc"

II

\(- \alpha \beta^T \mu\)

0

"ci"

III

\(\neq 0\)

0

"co"

IV

\(\neq 0\)

\(- \alpha \beta^T \gamma\)

"coli"

V

\(\neq 0\)

\(\neq 0\)

"colo"

VECM(endog[, exog, exog_coint, dates, freq, ...])

表示向量误差修正模型 (VECM) 的类。

coint_johansen(endog, det_order, k_ar_diff)

VECM 的协整秩的 Johansen 协整检验。

JohansenTestResult(rkt, r0t, eig, evec, lr1, ...)

Johansen 协整检验的结果类。

select_order(data, maxlags[, deterministic, ...])

基于每个可用的信息准则计算滞后阶数选择。

select_coint_rank(endog, det_order, k_ar_diff)

计算 VECM 的协整秩。

VECMResults(endog, exog, exog_coint, k_ar, ...)

用于保存向量误差修正模型 (VECM) 的估计相关结果的类。

CointRankResults(rank, neqs, test_stats, ...)

用于保存协整秩检验结果的类。

参考文献


最后更新:2024 年 10 月 3 日