含缺失值时间序列的ARMA模型拟合，本文主要内容关键词为：缺失论文,序列论文,模型论文,时间论文,ARMA论文，此文献不代表本站观点，内容供学术参考，文章仅供参考阅读下载。
在处理时间序列数据的工作中，序列太短会导致建模困难，取较长序列则可以保证拟合模型的可靠性。而实际情形是，序列涵盖的历史值越长，就越有可能含有缺失值。正如一般线性模型分析中对缺失值问题的研究，对这类时间序列数据的处理亦有两方面的考虑：一是缺失值的内插，二是模型参数的估计。随着随机过程理论的逐步完备，时间序列分析中缺失值问题的解决措施也日趋成熟[1～3]。利用状态空间方法描述时序状态，进而使用Kalman滤波技术，是同时解决这两个问题的有效手段[4]。
一、基本原理
（一）ARMA时序的状态空间表达
下式为零均值ARMA（p，q）模型的表达式：
其中ε[,t]表示不相关的随机噪声，方差为σ[2]， 可将这一过程的状态定义为阶数m＝Max（p，q+1）的列向量。
┌x（t│t）
┐
│x（t+1│t） │
Z（t）＝│· │
│· │
│· │
└ x（t+m－1│t） ┘
这里，x（t+j│t）表示依靠截止至时刻t时数据，对t+j 时刻取值的估计，即提前期ι＝j的j步预测，也可表示为E｛x（t+j）│x（s），s≤t｝。
基于t+1时刻前数据与t时刻前数据所做的j步预测有如下关系：
该递推公式表明，二者差值仅与t+1时刻的随机输入ε[,t+1]有关，不妨设：
x（t+j│t+1）－x（t+j│t）＝g[,j]ε[,t+1]
（1）
g[,j]由如下递推公式求得
┌g[,1]＝1
│j-1（2）
└g[,j]＝β[,j-1]+∑a[,k]g[,j-k]
k=1
其中，当j＞q，则β[,j]＝0。
状态向量最后一个元素为：
p
x（t+m│t+1）＝∑ a[,k]x（t+m－k│t）+g[,m]ε[,t+1]（3）
k=1
由式（1）与（3）得状态方程的矩阵表达式
Z（t+1）＝FZ（t）+Gε[,t+1] （4）从而构成一个Markov过程。F为m×m阶状态转移函数，F中最后一行为自回归参数。向量G的元素由自回归和移动平均参数算得（见式2）。
ARMA时序的状态空间表达中，另有观测值方程：
y（t）＝HZ（t）+v（t）（5）
v（t）是观测误差，是一随机变量，不同时刻取值互不相关，与｛ε[ ,t]｝也不相关，且具零均值与恒方差R。将ARMA 时序的结构特征用（4）与（5）两方程描述，进而可得如下对数似然函数：
该函数是求解a[,1]，a[,2]，…，a[,p]，β[,1]，β[,2]，…，β[,q]及R的依据。
（二）初始状态协方差阵
初始状态向量是零向量，而初始状态协方差阵的向量元素为：
P[,0j]（0│0）＝C（j）对于j≥i＞0，
i-1P[,ij]（O│O）＝C（j－i）－ ∑ γ[,k]γ[,k+j-i]
k=0
协方差C[,i]的求解则是通过解如下线性方程组获得：
┌C[,0]＝a[,1]C[,1]+a[,2]C[,2]+…+a[,p]C[,p]+γ[,0]+
│ β[,1]γ[,1]+β[,2]γ[,2]+…+β[,q]γ[,q]
│C[,1]＝a[,1]C[,2]+a[,2]C[,3]+…+a[,p]C[,p-1]+β[,1]
│γ[,0]+β[,2]γ[,1]+β[,3]γ[,2]+…+β[,q]γ[,q-1]
│……
└C[,p]＝a[,1]C[,p-1]+a[,2]C[,p-2]+…+a[,p]C[,0]+β[,
p]γ[,0]+β[,p+1]γ[,1]+β[,p+2]γ[,2]+…+β[
,q]γ[,q-p]
方程组中的互协方差γ[,i]乃由递推公式求解。
（三）缺失值存在时的参数估计与缺失值估计
当有缺失值存在时，不妨设y（t+1）缺失， 则该时刻的状态向量与协方差阵形式分别被修正为：
Z（t+1│t+1）＝Z（t+1│t） （7）
P（t+1│t+1）＝P（t+1│t） （8）
式（6）对数似然函数中对应t+1时刻的求和项被略去。 即依照对数似然函数（6）进行参数估计时所利用的信息，仍由原始序列提供。
早期对ARMA时序缺失值的插值问题是从最小二乘的角度加以考察。而在状态空间为研究框架的分析中，缺失值的内插值及相应的误差方差值可以方便地用与Kalman滤波相联系的递归平滑方法导出[5]。例如，在一个AR（p）过程中，若T时刻缺失，此前与此后均有至少p个以上连续的观测值，缺失值Z[,T]的最优插值由下式给出：
二、实例分析
实际工作中，基于Jones提供的算法，利用SAS 6.12 软件系统可以方便地完成时序的参数估计与缺失值的内插[6]。
表1 1987～1997 年山西医科大学第二医院设备科输液器出库量（单位：只）
年份/季度
出库量
年份/季度
出库量
年份/季度
出库量
出库量
1987.1 1500 1989.427000 1992.3 67000
22000
1987.2 1000 1990.127444 1992.4
... 47400
1987.3 1000 1990.236000 1993.1 76800
34625
1987.4 1500 1990.353700 1993.2
600
38200
1988.118000 1990.421400 1993.3 85000
30800
1988.2 3000 1991.141400 1993.4 52000
46300
1988.330000 1991.2 6000 1994.1 29000
44000
1988.424500 1991.336000 1994.2 30000
36800
1989.418000 1991.449000 1994.3 67500
47100
1989.233200 1992.122500 1994.4 40400
49620
1989.325960 1992.240500 1995.1 49300
40000
山西医科大学第二医院设备供应科1987年1季度～1997年4季度一次性输液器出库量如表1，在1992年第4季度有一次缺失值。利用SAS 6.12软件系统对其ARMA模型拟合，并对拟合残差白噪声检验，迟滞lag≥1时，残差序列自相关函数值与零的差别均无统计学意义（P＞0.05）， 拟合效果可靠。求得拟合模型为ARIMA（2，1，0）（1，1，0）s，自回归参数α[,1，1]＝－1.064，α[,1，2]＝－0.522，季节自回归参数α[,2，1]＝－0.523，各参数值均有统计学意义（P＜0.005）。 所得拟合曲线如图1。该图还给出了提前期ι＝1～12的预测值及其95%置信区间，表2中列出ι＝1～2的预测值及误差分析。1992年第4季度的估计值为37297只。
表2 输液器出库量预测结果
年度/季度 实际值 预测值95%可信服
绝对误差 相对误差(%)
1998.158400
49878 (18156,81600) -8522 -14.59
1998.240000
39701
(7914,71489) -299 -0.75
三、讨论
当平稳时序中存在缺失值时，使用状态空间的Markov表达及Kalman递归滤波算法进行ARMA模型拟合并无技术上的困难。模型识别的判别依据是对数似然函数或AIC准则等。
图1 输液器出库量序列的ARMA模型拟合及预测
缺失值存在时形如（7）与（8）的修正，可适用于有连续若干缺失值的情形。不过，当j→∞，Z（t+j│t）趋于零向量；P（t+j│t ）趋于P（0│0），即初始状态协方差阵。因此， 大块缺失数据之后进行的递归运算，等效于在新的起点（此时刻之前的历史值几乎不提供后续递归运算可资利用的信息）重新开始递归运算。
本文提供的算法，其应用前提是时间序列满足平稳可逆条件，即：方程1－
＝0与1+
＝0的根均在单位圆外，当此条件未能满足时，可以文献[4]中提供的参数修正方法，进而对之进行非线性最优估计。