基于Bayes推断的ARIMAX模型的建模分析

Modeling Analysis of ARIMAX Model Based on Bayes Inference

摘要: ARIMAX模型在多个领域有着重要应用。但是针对ARIMAX模型的参数估计均是经典统计方法，利用贝叶斯估计是一个值得探究的问题。考虑到共轭先验分布的性质特点，对于ARIMAX模型中的系数的先验分布为正态分布，噪声项的方差先验分布为逆伽马分布假定，本文给出了参数的后验分布，并使用Gibbs采样的方式，给出各个参数的一个估计。模拟试验的结果表明，本文的估计方法具有很好的功效，借助参数的迭代图表明文中使用的方法具有稳健性。使用创业板数据与上证数据做出实证分析，并发现文中给出的方法不仅具有很好的解释性，同时能够提取出完整的数据信息。

Abstract: The ARIMAX model has important applications in various fields. However, parameter estimation for the ARIMAX model has traditionally relied on classical statistical methods, and exploring Bayesian estimation presents a worthwhile direction. Considering the properties of conjugate prior distributions, this study assumes a normal prior distribution for the coefficients in the ARIMAX model and an inverse gamma prior for the variance of the noise term. The posterior distributions of the parameters are derived, and estimates for each parameter are obtained using Gibbs sampling. Simulation results demonstrate the effectiveness of the proposed estimation method. Iterative plots of the parameters indicate the robustness of the method. An empirical analysis using data from the ChiNext Index and the Shanghai Composite Index further reveals that the proposed method not only offers strong interpretability but also effectively captures comprehensive information from the data.

1. 引言

Box和Tiao (1975) [1]针对非平稳时间序列数据，首先将ARMA模型(Autoregressive Moving Average, ARMA)结合差分运算，提出求和移动平均自回归模型(Autoregressive Integrated Moving Average, ARIMA)，但是很多数据不仅仅受自身的影响，还要受到其他多个变量的影响。为此Engle和Granger (1987) [2]在ARIMA的基础上，引入外生变量$\left\{x_{1t}\right\},\left\{x_{2t}\right\},\cdot \cdot \cdot ,\left\{x_{kt}\right\}$ ，提出ARIMAX模型(Autoregressive Integrated Moving Average with Exogenous Variables, ARIMAX)。研究人员发现，使用ARIMAX模型在经济学与金融领域往往有良好的预测以及拟合效果，例如Peter和Silvia (2012) [3]使用人均GDP作为响应变量，失业率作为外生变量，发现ARIMAX的解释性要优于ARIMA模型；Sharma和Parkhar (2021) [4]借助ARIMAX模型对股票市场预测进行敏感度分析。但是在使用ARIMAX模型应用之前，各项系数的确定是一个值得探究的问题。传统的估计手段是Engle和Granger (1987) [2]使用协整回归的方式进行检验，不同于传统的方式，本文参考夏强和刘金山(2011) [5]对TAR模型门限变量使用贝叶斯的思想进行检验的手段，对ARIMAX模型的系数进行估计。假定先验分布已知，对于计算后验概率的方式，本文采取马尔可夫蒙特卡洛模拟(MCMC)方法中的Gibbs采样的方式计算。

本文剩余章节安排如下：第二章介绍ARIMAX模型与各个参数的后验分布，第三章为数值模拟，包含数值分析情形与实证数据分析，第四章是对文章的总结与概括。

2. 模型及估计方法

2.1. ARIMAX模型介绍

假设响应序列$\left\{y_t\right\}$ 和输入序列(自变量序列)$\left\{x_{1t}\right\},\left\{x_{2t}\right\},\cdot \cdot \cdot ,\left\{x_{kt}\right\}$ 平稳，且响应序列和输入序列之间具有相关关系。响应序列$\left\{y_t\right\}$ 和输入序列$\left\{x_{1t}\right\},\left\{x_{2t}\right\},\cdot \cdot \cdot ,\left\{x_{kt}\right\}$ 最终建立的线性模型为：

$y_t={\beta }_0+{\displaystyle \sum _{i=1}^k\frac{{\Psi }_i\left(B\right)}{{\Phi }_i\left(B\right)}}{B}^{r_i}{x}_{it}+\frac{\Psi \left(B\right)}{\Phi \left(B\right)}{a}_t$ (1)

其中$t\in \{1,2,\cdot \cdot \cdot ,n\}$ 为时间参数，$B$ 代表一阶滞后算子，$r_i$ 表示第$i$ 个变量的延迟阶数，$a_t$ 表示零均值白噪声序列，以及${\Phi }_i(B)$ 表示第$i$ 个自变量序列$\left\{x_{it}\right\}$ 的$p_i$ 阶自回归系数多项式：

${\Phi }_i\left(B\right)=1-{\varphi }_{i1}B-{\varphi }_{i2}{B}^2-\cdot \cdot \cdot -{\varphi }_{i{p}_i}{B}^{p_i},i=1,\cdot \cdot \cdot ,k$

且${\Psi }_i(B)$ 是第$i$ 个自变量序列$\left\{x_{it}\right\}$ 的$m_i$ 阶移动平均系数多项式：

${\Psi }_i\left(B\right)={\psi }_{i0}-{\psi }_{i1}B-\cdot \cdot \cdot -{\psi }_{i{m}_i}{B}^{m_i},i=1,\cdot \cdot \cdot ,k$

经过一系列变换，模型(1)可以转变为

$y_t={\varphi }_1{y}_{t-1}+{\varphi }_2{y}_{t-2}+\cdot \cdot \cdot +{\varphi }_p{y}_{t-p}+{\displaystyle \sum _{i=1}^k{\Theta }_i\left(B\right)B^{l_i}{x}_{it}}+{\epsilon }_t$ (2)

此时，式中${\Theta }_i(B)={\theta }_{i0}+{\theta }_{i1}B+\dots +{\theta }_{i{q}_i}{B}^{q_i}$ 为第$i$ 个自变量$\left\{x_{it}\right\}$ 的$q_i$ 阶移动平均系数多项式，$q_i$ 为第$i$ 个自变量$x_{it}$ 对变量$y_t$ 的有效作用时长，${\epsilon }_t$ 服从标准正态分布。现在模型中，未知的参数为$\left\{y_t\right\}$ 的自回归系数，

记$\Phi ={\left({\varphi }_1,\cdot \cdot \cdot ,{\varphi }_p\right)}^{\prime }$ ；以及自变量序列的移动平均系数，记为${\Theta }_i={\left({\theta }_{i_0},\cdot \cdot \cdot ,{\theta }_{i{q}_i}\right)}^{\prime },i=1,\cdot \cdot \cdot ,k$ 。

2.2. ARIMAX模型参数的贝叶斯推断

为了对参数进行贝叶斯推断，首先要有对参数的先验分布的认识，这样才能得到参数的条件后验分布(Gelfand和Smith (1990) [6])，考虑到共轭先验分布的性质特点[7]，假定参数向量$\Phi$ 服从多元正态分布$N_p\left(\mu ,\Sigma \right)$ ；${\Theta }_i$ 之间互相独立，并且服从多元正态分布$N_{q_i}\left({\alpha }_i,{\Xi }_i\right)$ ，噪声项的方差${\sigma }^2$ 服从逆伽马分布$LG\left(\beta ,\gamma \right)$ 。

为书写方便起见，将模型(2)中的每个观测按行摆放，且不妨设$t$ 足够大，记

$Y_0=\left(\begin{array}{ccc}{y}_{t-1}& \cdots & y_{t-p}\\ ⋮& \ddots & ⋮\\ y_{n-1}& \cdots & y_{n-p}\end{array}\right),X_i=\left(\begin{array}{ccc}{x}_{i,t-l_i}& \cdots & x_{i,t-l_i-q_i}\\ ⋮& \ddots & ⋮\\ x_{i,n-l_i}& \cdots & x_{i,n-l_i-q_i}\end{array}\right),i=1,\cdot \cdot \cdot ,k$

并记因变量$Y=(y_t,\cdot \cdot \cdot ,y_n)\prime$ ，于是模型(2.2)可化为：

$Y=Y_0\Phi +{\displaystyle \sum }_{i=1}^k{X}_i{\Theta }_i+\epsilon$

其中$\epsilon =({\epsilon }_t,\cdot \cdot \cdot ,{\epsilon }_n{)}^{\prime }$ ，且服从多元正态分布$N_{n-t+1}\left(0,{\sigma }^2{I}_{n-t+1}\right)$ 。因此可以得到参数的联合先验分布：

$\pi (\Phi ,{\Theta }_i,{\sigma }^2)\propto \frac{1}{{\sigma }^n}\exp \left\{-\frac{1}{2{\sigma }^2}\left[{\left(Y-Y_0\Phi -{\displaystyle \sum _{i=1}^k{X}_i{\Theta }_i}\right)}^{\prime }\left(Y-Y_0\Phi -{\displaystyle \sum _{i=1}^k{X}_i{\Theta }_i}\right)\right]\right\}$

于是可以得到参数${\sigma }^2$ 条件后验分布为逆伽马分布：

$P\left({\sigma }^2|\Phi ,{\Theta }_i\right)~LG\left(\beta +\frac{n}{2},\gamma +E^{\prime }E\right)$ (3)

其中

$E=Y-Y_0\Phi -{\displaystyle \sum _{i=1}^k{X}_i{\Theta }_i}$

参数$\Phi$ 的后验分布为

$P(\Phi |{\sigma }^2,{\Theta }_i)~{\mathbb{N}}_p({\mu }^{\ast },{\Sigma }^{\ast })$ (4)

其中

${\mu }^{\ast }={\left(Y_0{\Sigma }_0^{-1}{Y}_0+{\Sigma }^{-1}\right)}^{-1}\left(Y_0{}^{\prime }{\Sigma }_0^{-1}E+\mu \right),{\Sigma }^{\ast }={\left(Y_0{\Sigma }_0^{-1}{Y}_0+{\Sigma }^{-1}\right)}^{-1}$

参数${\Theta }_i,i=1,\cdot \cdot \cdot ,k$ 的后验分布为正态分布

$P({\Theta }_i|{\sigma }^2,\Phi )~N_{q_i}\left({\alpha }_i{}^{\ast },{\Xi }_i{}^{\ast }\right)$ (5)

其中

${\alpha }_i^{\ast }={\left(X_i{}^{\prime }{\Sigma }_0^{-1}{X}_i+{\Xi }_i^{-1}\right)}^{-1}\left(X_i{}^{\prime }{\Sigma }_0^{-1}E+{\Xi }_i^{-1}{\alpha }_i\right),{\Xi }_i{}^{\ast }={\left(X_i{}^{\prime }{\Sigma }_0^{-1}{X}_i+{\Xi }_i^{-1}\right)}^{-1}$

可以发现各个参数的后验分布都是标准的分布，那么就可以使用Gibbs采样[8]的方法，得到各个参数的估计量。

3. 数值分析

3.1. 试验模拟

将模型设置为：

$y_t=0.5y_{t-1}+0.4y_{t-2}+0.5x_{t-1}+0.3x_{t-2}+{\epsilon }_t$

${\epsilon }_t$ 服从标准正态分布。自变量序列通过下述方式迭代产生：

$x_t=0.5x_{t-1}+{\rho }_t$

Figure 1. Diagram of parameter $\Phi ,\Theta ,{\sigma }^2$ iteration

图1. 参数$\Phi ,\Theta ,{\sigma }^2$ 的迭代图

${\rho }_t$ 服从标准正态分布，初始值$x_1$ 为一个(0,1)的均匀分布的随机数。根据模型(2)，可以观察到各个参数的真实值为${\sigma }^2=1,\Phi ={\left(0.5,0.4\right)}^{\prime },\Theta ={\left(0.5,0.3\right)}^{\prime }$ 。由于超参数对实验结果的选择影响不大，见文献Casella (2001) [9]与Doss (2012) [10]，于是可以设置噪声项的方差${\sigma }^2$ 的先验分布为$LG\left(1.5,1\right)$ ，系数$\Phi$ 与$\Theta$ 均为标准多元正态分布。采用Gibbs采样的方法，估计各个参数。但是使用Gibbs抽样时的一个问题是，在什么时候抽样可以停止，即收敛，为此我们设置迭代一万次，并舍弃前五千次的数据，使用后五千次的数据进行后验分析，并绘制迹图。算法流程如下：

算法1. 试验流程

(1) 生成100个数据，记为$y_1,y_2,\cdot \cdot \cdot ,y_{100}$ ，将第3到第100个数据作为因变量向量$Y={\left(y_3,y_4,\cdot \cdot \cdot ,y_{100}\right)}^{\prime }$ ，设置初始值

$Y_0={\left(\begin{array}{l}{y}_1,y_2,\cdot \cdot \cdot ,y_{98}\\ y_2,y_3,\cdot \cdot \cdot ,y_{99}\end{array}\right)}^{\prime }$ ，$X_0={\left(\begin{array}{l}{x}_1,x_2,\cdot \cdot \cdot ,x_{98}\\ x_2,x_3,\cdot \cdot \cdot ,x_{99}\end{array}\right)}^{\prime }$ ；

(2) 利用(3)抽取${\sigma }^2$ ，利用(4)抽取$\Phi$ ，利用(5)抽取参数$\Theta$ ；

(3) 重复(1)与(2) 10000次；

(4) 剔除前5000次数据，求后5000次的均值，得到参数的估计${\widehat{\sigma }}^2,\widehat{\Phi },\widehat{\Theta }$ ；

(5) 重复(1)至(4) 100次，对得到的数据求均值以及标准差。

表1中给出各个参数的估计值，可以发现，通过Gibbs取样，在进行一百次抽样试验中，可以把误差控制在0.05以内。图1为使用Gibbs采样时，迭代的迹图，在迭代5000次之后发现参数已经是收敛的，因此认为估计是稳健的。

Table 1. Bayes estimates of parameter ${\sigma }^2,\Phi ,\Theta$ under 100 simulations

表1. 参数${\sigma }^2,\Phi ,\Theta$ 贝叶斯估计值在100次模拟下的估计值

参数	真值	实验结果	误差	标准差
${\sigma }^2$	1.0	1.0811	0.0811	0.1855
$\Phi$	0.50.4	0.45490.4242	−0.0451 −0.0309	0.07330.0718
$\Theta$	0.50.3	0.48670.3492	−0.01330.0492	0.10480.1279

在表2中，给出协整的方式估计的参数值，通过对比可以发现使用Bayes的方式进行估计，大部分情况下优于经典方法(Engle和Granger (1987) [2])使用协整的方式进行估计。

Table 2. Co-integration estimates of parameter ${\sigma }^2,\Phi ,\Theta$ co-integration estimates under 100 simulations

表2. 参数${\sigma }^2,\Phi ,\Theta$ 协整估计值在100次模拟下的估计值

参数	真值	实验结果	误差	标准差
${\sigma }^2$	1.0	1.129	0.129	0.7975
$\Phi$	0.50.4	0.57690.3757	0.0769 −0.0343	0.25560.1613
$\Theta$	0.50.3	0.43570.3419	−0.06430.0419	0.18510.1186

3.2. 实证分析

借助ARIMAX模型，对金融数据进行分析建模的有很多，例如Badshah等人(2023) [11]，Zolfaghari和Gholami (2021) [12]。本文选取自2019年2月18日至2025年3月5日的创业板指数与上证指数数据，选取上证指数作为外生变量序列，创业板指数作为响应变量序列。为使数据更加平稳，对数据采取对数差分的手段。针对自回归阶数与移动平均阶数$p,q$ 的确定，采用信息准则AIC、BIC、AICc最小的方式进行确定，发现这是一个带有外生变量的ARIMA(1,0,0)的模型。外生变量的滞后阶数采用逐步回归的方式，基于AIC最小的准则，发现外生变量的滞后阶数为0。在进行参数的估计时，对先验分布设置为${\sigma }^2$ 服从逆伽马分布$LG\left(1.5,1\right)$ ，其余两个参数均服从标准正态分布。得到的拟合结果为

$Y_t=0.0915Y_{t-1}+1.3496X_t+{\epsilon }_t$

Figure 2. Diagram of autocorrelation in measure error

图2. 测量误差的自相关图

图2给出模型拟合的残差的自相关函数图。表3给出针对实证数据模型拟合的误差输出，从图2以及残差的自相关系数为−0.0017十分接近于0得知，可以通过白噪声检验，同时也表明ARIMAX模型对于创业板指数与上证指数的解释性很强，同时将模型中的具有同响应变量具有相关信息的数据提取得十分充分。

Table 3. Errors of model

表3. 模型误差

残差方差	AIC	AICc	BIC	平均误差	均方误差	ACF
0.0001	−9131.55	−9131.53	−9115.68	−0.0001	0.0107	−0.0017

4. 结论

从设定的模型出发，基于贝叶斯推断，文中给出一种新的针对ARIMAX模型参数估计的手段，这个方法主要利用Gibbs采样计算后验概率分布，从后验分布中抽取的均值就是参数的一个估计。在模拟试验中发现，文中的方法的功效很好，并且估计的方法具有稳健性，使用这种方法在实证分析中发现对上证指数与创业板指数具有很高的解释性。