arma

2019-12-02 本文已影响0人围城17

<center>1.前言</center>

本文设计了一个ARMA模型估计算法，该算法首先将ARMA模型转换为卡尔曼滤波器状态空间模型，在给定模型阶数 $0 \leq p,q \leq L$ 下，先使用Hannan–Rissanen算法，从观测数据 $y_t=\{y_0,y_1,\cdots,y_{N-1}\}$ 得到模型的初始估计，再以该初始估计为起点，使用非线性最优化得到模型的极大似然估计。得到模型的最优估计后，使用BIC定阶准则选取合适的模型阶数，最后使用卡尔曼滤波器对时间序列进行预报。本文使用的 $ARMA(p,q)$ 模型定义为:
$y_t = \phi_1 y_{t-1}+\cdots+\phi_r y_{t-r}+a_t + \theta_1 a_{t-1} + \cdots + \theta_{t-r}a_{t-r} \tag{1}$
其中 $r=max\{p,q+1\}$ ,对于 $i>p$ 有 $\phi_i=0$ ,对 $j>q$ 有 $\theta_j=0$ , $a_t$ 为零均值的白噪声，即 $E\left[a_t\right]=0,E\left[a_t a_s\right] \begin{cases} \sigma_a^2, \quad t=s \\ 0, \quad t \neq s \end{cases},E\left[a_s y_t\right]=0$

<center>2.卡尔曼滤波器设计</center>

本文使用的卡尔曼滤波器模型定义如下 $^{\left[1\right]}$ :
$\begin{align} x_{t+1} &= Ax_t+Ra_{t+1} \tag{2}\\ y_t &= Zx_t \tag{3} \end{align}$

其中 $x_t \in R^{r\times 1}$ 为状态向量， $y_t \in R^{1 \times 1}$ 为ARMA模型序列数据， $A \in R^{r \times r}$ 为状态转移矩阵， $R \in R^{r\times 1}$ ， $Z \in R^{1 \times r}$ 为观测矩阵，这些向量和矩阵定义如下：
$A=\left[ \begin{matrix} \phi_1 &1 &0 &0 &\cdots &0 \\ \phi_2 &0 &1 &0 &\cdots &0 \\ \phi_3 &0 &0 &1 &\cdots &0 \\ \vdots &\vdots &\vdots &\vdots &\ddots &\vdots \\ \phi_{r-1} &0 &0 &0 &\cdots &1 \\ \phi_{r} &0 &0 &0 &\cdots &0 \\ \end{matrix} \right]; R=\left[ \begin{matrix} 1 \\ \theta_1 \\ \theta_2 \\ \vdots \\ \theta_{r-1} \end{matrix} \right]; Z = \left[ \begin{matrix} 1 &0 &0 &\cdots &0 \end{matrix} \right];$
则最后一行可表示为:
$x_{r,t+1} = \phi_r x_{1,t} + \theta_{r-1}a_{t+1} \tag{4}$
延迟 ${r-1}$ 得到:
$x_{r,t-r+2} = \phi_r L^{r-1} x_{1,t} + \theta_{r-1}L^{r-1}a_{t+1} \tag{5}$
其中 $L^r x_t=x_{t-r}$ ,取 $r-1$ 行得到：
$x_{r-1,t+1} = \phi_{r-1} x_{1,t} + x_{r,t} + \theta_{r-2}a_{t+1} \tag{6}$
延迟 $r-2$ 得到:
$x_{r-1,t-r+3} = \phi_{r-1} L^{r-2} x_{1,t} + x_{r,t-r+2} + \theta_{r-2}L^{r-2}a_{t+1} \tag{7}$
代入(5)得到：
$x_{r-1,t-r+3} = \phi_{r-1} L^{r-2} x_{1,t} + \phi_r L^{r-1} x_{1,t} + \theta_{r-1}L^{r-1}a_{t+1} + \theta_{r-2}L^{r-2}a_{t+1} \tag{8}$
或表示为:
$x_{r-1,t-r+3} = \left(\phi_{r-1} L^{r-2}+ \phi_r L^{r-1} \right)x_{1,t} x_{1,t}+ \left(\theta_{r-1}L^{r-1} + \theta_{r-2}L^{r-2}\right)a_{t+1} \tag{9}$
取 $r-2$ 行得到:
$x_{r-2,t+1} = \phi_{r-2} x_{1,t} + x_{r-1,t} + \theta_{r-3}a_{t+1} \tag{10}$
延迟 $r-3$ 得到:
$x_{r-2,t-r+4} = \phi_{r-2} L^{r-3} x_{1,t} + x_{r-1,t-r+3} + \theta_{r-3}L^{r-3}a_{t+1} \tag{11}$
带入(9)得到：
$x_{r-2,t-r+4} = \left[\phi_{r-2} L^{r-3}+\phi_{r-1} L^{r-2}+\phi_{r} L^{r-1} \right]x_{1,t} + \left[\theta_{r-1}L^{r-1}+\theta_{r-2}L^{r-2}+\theta_{r-3}L^{r-3}\right]a_{t+1} \tag{12}$
依次类推在第一行有：
$x_{1,t+1} = \left[\phi_{1}+ \cdots + \phi_{r-2}L^{r-3}+\phi_{r-1} L^{r-2}+\phi_{r} L^{r-1} \right]x_{1,t} + \left[\theta_{r-1}L^{r-1}+\theta_{r-2}L^{r-2}+\theta_{r-3}L^{r-3}+\cdots+1\right]a_{t+1} \tag{13}$
或:
$(1-\phi_{1}L - \phi_{2}L^2 - \cdots - \phi_r L^{r})x_{1,t+1}= \left[\theta_{r-1}L^{r-1}+\theta_{r-2}L^{r-2}+\theta_{r-3}L^{r-3}+\cdots+1\right]a_{t+1} \tag{14}$
由式(4)和 $Z$ 的定义有:
$y_t = x_{1,t} \tag{15}$
带入(14)得到ARMA模型：
$(1-\phi_{1}L - \phi_{2}L^2 - \cdots - \phi_r L^{r})y_t= \left[\theta_{r-1}L^{r-1}+\theta_{r-2}L^{r-2}+\theta_{r-3}L^{r-3}+\cdots+1\right]a_t \tag{16}$
证明(2)(3)式等价于(1)式。
在建立卡尔曼滤波器模型后，便可使用递推公式进行计算：
$\begin{align} \hat{x} &=A\hat{x}_{t-1}+H_t(y_t-ZA\hat{x}_{t-1}) \tag{17} \\ P^{'}_t &=A P_{t-1} A^T + Q \tag{18}\\ H_t &= P^{'}_tZ^T(ZP^{'}_tZ^T)^{-1} \tag{19} \\ P_t &= (I-H_tZ)P^{'}_t \tag{20} \end{align}$
其中 $Q=RR^T\sigma_a^2$ 为噪声的协方差矩阵, $H_k$ 为增益矩阵

<center>3.模型的极大似然估计</center>

当卡尔曼滤波器得到估计值 $\hat{x}_t$ ,可以得到新息 $\Delta_t=y_t - Z\hat{x}_t$ ，新息的方差为 $\omega_t=E\left[(y_t-Z\hat{x}_t)(y_t-\hat{x}_t)^T\right]=E\left[(Zx_t-Z\hat{x}_t)(Zx_t-\hat{x}_t)^T\right]=ZP^{'}Z^T$ ,则观测向量 $\{y_0,y_1,...y_{N-1}\}$ 的似然函数表示为 $^{\left[1\right]}$ ：
$L=\prod_{t=0}^{N-1}(2 \pi \omega_t)^{-1/2}exp\left(-\frac{\Delta_t^2}{2\omega_t}\right) \tag{21}$
取对数得到：
$l=-\sum_{t=0}^{N-1}\left[ln(\omega_t)+\Delta_t^2/\omega_t\right] \tag{22}$
此时似然函数的参数为 $\phi_i,\theta_j,\sigma_a$ ,将 $Q$ 中的 $\sigma_a$ 移入 $P^{'}_t$ 得到
$l=-\sum_{t=0}^{N-1}\left[ln(\sigma_a^2 \omega_t)+ \frac{\Delta_t^2}{\omega_t\sigma_a^2}\right] \tag{23}$
以 $\sigma_a^2$ 为参数对 $l$ 最大化得到
$\sigma_a^2=\frac{1}{N}\sum_{t=0}^{N-1}\Delta_t^2/\omega_t \tag{24}$
带入式(22)得到
$\begin{align} l &=-\sum_{t=0}^{N-1}\left[ln(\sigma_a^2)+ln(\omega_t)+ \frac{\Delta_t^2}{\omega_t\sigma_a^2}\right] \\ \tag{25} &=-\left[Nln(1/N)+N +Nln\sum_{t=0}^{N-1}\Delta_t^2/\omega_t+\sum_{t=0}^{N-1}ln\omega_t\right] \\ \end{align}$
去除常数项得到：
$l=-\left[Nln\sum_{t=0}^{N-1}\Delta_t^2/\omega_t+\sum_{t=0}^{N-1}ln\omega_t\right] \tag{26}$
该函数是参数 $\phi_i,\theta_j$ 的非线性函数，需要使用非线性优化算法求函数的最大值。

<center>4.初始参数确定</center>

由于初始参数的选择对优化算法的效率和结果有较大影响，因此需要对模型的参数进行初步估计。本文使用了Hannan–Rissanen线性回归算法 $^{[2][3]}$ 进行初步估计。该算法分为两步:

先将观测序列当做高阶 $AR(p)$ 模型：
$y_t = \hat{\phi}_{h1}y_{t-1} + \cdots + \hat{\phi}_{hh}y_{t-h} + a_t \tag{27}$
使用尤尔-瓦尔克方程求解模型参数，并使用BIC准则确定模型的阶数 $h$ ：
$BIC(h)=log(\hat{\sigma}_h^2 +hlogN/N \tag{28}$
其中 $\hat{\sigma}_h^2$ 为估计的噪声方差， $N$ 观测序列长度。
则估计噪声 $\hat{a}_t$ 可表示为
$\hat{a}_t = y_t - \hat{\phi}_{h1}y_{t-1} + \cdots + \hat{\phi}_{hh}y_{t-h} \tag{29}$
在得到估计噪声序列后，计算参数向量 $\beta=(\phi,\theta)$ 在线性回归模型序列 $(X_{t-1},...,X_{t-p},\hat{a}_{t-1},...,\hat{a}_{t-q}),t=h+1+q,...,N-1$ 的最小二乘解，也就是最小化：
$S(\beta)=\sum_{t=h+1+q}^{N-1}(y_t - \phi_1 y_{t-1}-\cdots - \phi_p X_{t-p} - \theta_1 \hat{a}_{t-1}-\cdots-\theta_q \hat{a}_{t-q})^2 \tag{30}$
该方程实际也是通过尤尔-瓦尔克方程求解。

实验证明该估计结果较接近最优化结果。

<center>4.BIC定阶</center>

由于ARMA模型不可能唯一确定，模型的阶数也并不唯一 $^{\left[3\right]}$ ，需要根据实际情况选择相应的准则确定最合适的阶数，如FPE准则,AIC准则,BIC准则等,本文使用了BIC准则选取模型阶数，对于ARMA模型，BIC准则表示为:
$\begin{align*} BIC &= (N-p-q)ln\left[N\hat{\sigma}_a^2/(n-p-q)\right] + N\left(1+ln\sqrt{2\pi}\right) \tag{31} \\ &+(p+q)ln\left[\left(\sum_{t=0}^{N-1}y_t^2-N\hat{\sigma}_a^2\right)/(p+q)\right] \tag{32} \end{align*}$
其中 $\hat{\sigma}_a^2$ 为白噪声 $a_t$ 方差的极大似然估计。模型的阶数：
$p,q = \underset{0\leq p,q \leq L}{minarg }BIC(p,q)$
其中 $L$ 为模型的最高阶数。

<center>5.关于模型的跟踪</center>

由于算法中使用了卡尔曼滤波器进行递推运算，因此在得到序列的新观测数据后可以加入该滤波器进行递推计算并进行新一轮优化，由于系统具有一定的惯性，因此新的模型参数数在上一次参数的附近，因此可以用上一次的结果作为优化的起点。

<center>5.算法流程</center>

由于算法涉及大量矩阵运算，因此程序使用了开源矩阵运算库eigen，非线性优化部分使用了google开源的非线性优化库ceres。

程序先使用generate_arma_data函数生成长度为200的arma序列，接着调用arma_estimator_hr函数(Hannan–Rissanen算法)求解模型的初步估计，最后由arma_estimator_mle函数求解极大似然估计，得到最优估计后根据BIC准则计算BIC结果，得到最小的 $BIC(p,q),1\leq p,q \leq 4$ ：程序运行结果如下：

:p=1,q=1,BIC=799.489580544931
:p=1,q=2,BIC=806.092794685896
:p=1,q=3,BIC=812.357735217179
:p=1,q=4,BIC=817.981515183498
:p=2,q=1,BIC=805.966371063034
:p=2,q=2,BIC=812.101541520479
:p=2,q=3,BIC=817.981515183498
:p=2,q=4,BIC=823.657602062738
:p=3,q=1,BIC=812.101541520479
:p=3,q=2,BIC=817.981515183498
:p=3,q=3,BIC=823.657602062738
:p=3,q=4,BIC=829.163703893365
:p=4,q=1,BIC=817.904762296210
:p=4,q=2,BIC=823.581092066115
:p=4,q=3,BIC=829.087436787408
:p=4,q=4,BIC=834.447889914660
Finish. p=1,q=1,coefs=0.774653 0.499346

coefs为最终结果， $\phi=0.774653$ (实际为0.8）， $\theta=0.499346$ (实际为0.5）

<center>8.参考文献

[1].Hevia, C., World, T., & Decrg, B. (2008). Maximum Likelihood Estimation of an ARMA ( p , q ) Model. October, (October), 1–6.

[2].McDougall A. Algorithms in time series[J]. 1988.

[3].Shanmugam, Ramalingam. Introduction to Time Series and Forecasting[J]. Technometrics, 39(4):426-426.