SAS学习系列39时间序列分析ⅢARIMA模型可编辑修改word版.docx

1、SAS学习系列39时间序列分析ARIMA模型可编辑修改word版39. 时间序列分析ARIMA 模型随着对时间序列分析方法的深入研究，人们发现非平稳序列的确定性因素分解方法（如季节模型、趋势模型、移动平均、指数平滑等） 只能提取显著的确定性信息，对随机性信息浪费严重，同时也无法对确定性因素之间的关系进行分析。而非平稳序列随机分析的发展就是为了弥补确定性因素分解方法的不足。时间序列数据分析的第一步都是要通过有效手段提取序列中所蕴藏的确定性信息。Box 和 Jenkins 使用大量的案例分析证明差分方法是一种非常简便有效的确定性信息的提取方法。而 Gramer 分解定理则在理论上保证了适当阶数的差

2、分一定可以充分提取确定性信息。（一）ARMA 模型即自回归移动平均移动模型，是最常用的拟合平稳时间序列的模型，分为三类：AR 模型、MA 模型和ARMA 模型。一、AR(p)模型p 阶自回归模型1.模型：xt = 0 + 1 xt-1 + p xt- p + t其中， 0 ，随机干扰序列 t 为 0 均值、 2 方差的白噪声序列（p E( t s ) = 0 , ts），且当期的干扰与过去的序列值无关，即 E(xtt)=0.由于是平稳序列，可推得均值 = 01 - - - . 若 0 = 0 ，称为中心化的 AR(p)模型， 对于非中心化的平稳时间序列， 可以令 = (1 - - - )， x

3、* = x- 转化为中心化。记 B 为延迟算子， (B) = I - B - - Bp 称为 p 阶自回归多项式，则 AR(p)模型可表示为： p (B)xt = t .2.格林函数用来描述系统记忆扰动程度的函数，反映了影响效应衰减的快慢程度（回到平衡位置的速度），Gj 表示扰动 t-j 对系统现在行为影响的权数。例如，AR(1)模型（一阶非齐次差分方程）， Gj= j ,j = 0,1, 2, 1模型解为 xt = Gj t- j .j=03.模型的方差 221对于 AR(1)模型，Var( xt ) = GjVar( t- j ) = .4.模型的自协方差j=01 - 2对中心化的平稳模型

4、，可推得自协方差函数的递推公式：用格林函数显示表示： (k ) = G G E( - - - ) = 2 G + Gi j t j t k j j k j i=0 j=0 j=0对于 AR(1)模型，1 1 (k ) = k (0) = k5.模型的自相关函数递推公式：2 11 - 21 1对于 AR(1)模型， (k ) = k (0) = k .平稳 AR(p)模型的自相关函数有两个显著的性质：（1）拖尾性指自相关函数 (k)始终有非零取值，不会在 k 大于某个常数之后就恒等于零；（2）负指数衰减i随着时间的推移，自相关函数 (k)会迅速衰减，且以负指数 k（其中 i 为自相关函数差分方程

5、的特征根）的速度在减小。6.模型的偏自相关函数自相关函数 (k)实际上并不只是 xt 与 xt-k 之间的相关关系，它还会受到中间 k-1 个随机变量 xt-1, , xt-k+1 的影响。为了能剔除了中间 k-1 个随机变量的干扰，单纯测度 xt 与 xt-k 之间的相关关系，引入了滞后 k 偏自相关函数（PACF），计算公式为：其中，滞后k 偏自相关函数实际上等于k 阶自回归模型第k 个回归系数 kk ：两边同乘以 xt-k，求期望再除以 (0) 得到取前 k 个方程构成的方程组：称为 Yule-Walker 方程，可以解出 kk .可以证明平稳 AR(p)模型，当 kp 时， kk型的偏

6、自相关函数具有 p 步截尾性。= 0 . 即平稳 AR(p)模注：实际上样本的随机性使得偏自相关函数不是严格截尾，例如上面两图都 1 阶显著不为 0，1 阶之后都近似为 0.二、MA(q)模型q 阶移动平均模型1.模型：其中， 0 ，随机干扰序列 为 0 均值、 2 方差的白噪声序列（q t E( t s ) = 0 , ts）。若 =0，称为中心化的 MA(q)模型，非中心化的 MA(q)模型可以t t通过 x* = x - 转化为中心化。q记 B 为延迟算子， (B) = I - B - - Bq 称为 q 阶自移动平均系数多项式，则中心化 MA(q)模型可以表示为 xt = q (B)

7、t .2.模型的方差3.模型的自协方差只与滞后阶数 k 相关，且 q 阶截尾。当 k=0 时，当 1kq 时，当 kq 时， (k ) = 0 .4.模型的自相关函数： (k ) = (k ) （q 阶截尾性） (0)5.模型的滞后 k 阶偏自相关函数（中心化）可以证明滞后 k 阶偏自相关函数具有拖尾性。6.模型的可逆性以 MR(1)为例，模型： x = - 或xt = tt t 1t-11 - 1B模型： x = - 1 或xt = t t t-111 - 1 B t 11 1 1它们的自相关函数 = - / (1 + 2 ) 相同（即相同的自相关函数对应不同的回归模型），为了保证对应的唯一

8、性，需要增加约束条件，即MR(q)模型的可逆性条件。观察两个模型的第二种表示：当| 1 | 1时，模型不收敛、模型收敛。表示成收敛形式的 MR(q)模型称为可逆 MR(q)模型。一个自相关函数只对应唯一一个可逆 MR(q)模型。三、ARMA(p, q)模型自回归移动平均模型1.模型其中， 0 ， 0 ，随机干扰序列 t 为 0 均值、 2 方差的白噪声p q 序列（ E( t s ) = 0 , ts），且当期的干扰与过去的序列值无关，即 E(xtt)=0.若 0 =0 ，则称为中心化的 ARMA(p,q)模型。引入延迟算子，中心化的 ARMA(p,q)模型可表示为： p (B)xt = q

9、(B) t .显然，AR(p)和 MA(q)模型是 ARMA(p,q)模型的特例。2.数字特征（1）均值： E( xt ) = 0 ;1 - - - （2）自协方差函数： (k ) = 2 G G，其中 Gi 为格林函数； i i+k i=0（3）自相关函数： (k ) = (k ) = (0)GiGi+ki=02ii=03.模型的初步定阶对于平稳非白噪声序列，计算出样本自相关系数（ACF）和偏自相关系数（PACF），根据其性质估计自相关阶数 p 和移动平均阶数q ，称为 ARMA(p,q)模型的定阶。kk可以推导出：样本自相关函数 (k ) 和偏自相关函数 都近似服1从正态分布N (0, )

10、 .n取显著水平 =0.05，若样本自相关系数和样本偏自相关系数在最初的 k 阶明显大于 2 倍标准差，而后几乎 95%的系数都落在 2 倍标准差的范围内，且非零系数衰减为小值波动的过程非常突然，通常视为k 阶截尾；若有超过 5%的样本相关系数大于 2 倍标准差，或者非零系数衰减为小值波动的过程比较缓慢或连续，通常视为拖尾。4.参数估计对非中心化的 ARMA(p,q)模型x = + q (B) .t (B) t参数 可用样本均值来估计总体均值（矩估计法），初步定阶估计出自相关阶数 p 和移动平均阶数q 后，模型共有 p+q+1 个未知参数： , , , , , , 2 .（1）参数的矩估计用时

11、间序列样本数据计算出延迟1 阶到p+q 阶的样本自相关函数 (k ) ，延迟 k 阶的总体自相关函数为 k ( 1 , , p , 1 , , q ) . 用计算 出的样本自相关函数来估计总体自相函数，得到 p+q 个联立方程组：从中解出 1 , , p , 1 , , q 的值作为未知参数估计值 , , , , , .ARMA(p,q)模型的两边同时求方差，并把前面的参数的估计值代入， 可得白噪声序列的方差估计为：（2）参数的极大似然估计当总体分布类型已知时，极大似然估计是常用的估计方法。其基本思想是，认为样本来自使该样本出现概率最大的总体。因此，未知参数的极大似然估计，就是使得似然函数（即

12、联合密度函数）达到最大值的参数值：在时间序列分析中，序列的总体分布通常是未知的。为了便于分析和计算，通常假设序列服从多元正态分布，它的联合密度函数是可导的。在求极大似然估计时，为了求导方便，常对似然函数取对数， 然后对对数似然函数中的未知参数求偏导数，得到似然方程组。理论上，只要求解似然方程组即可得到未知参数的极大似然估计。但在实际上是使用计算机经过复杂的迭代算法求出未知参数的极大似然估计。两种估计的比较：矩估计的优点是不要求知道总体的分布，计算量小，估计思想简单直观。但缺点是只用到了样本自相关系数的信息，序列中的其他信息被忽略了，这导致估计精度一般较差。因此，它常被作为极大似然估计和最小二乘

13、估计的迭代计算的初始值。极大似然估计的优点是充分应用了每一个观察值所提供的信息， 因而它的估计精度高，同时，还具有估计的一致性、渐近正态性和渐近有效性等优良统计性质，是一种非常优良的参数估计方法。（3）参数的最小二乘估计使 ARMA(p,q)模型的残差平方和达到最小的那组参数值：通过计算机借助迭代方法求出。由于充分利用了序列的信息，该方法估计精度最高。在实际运用中，最常用的是条件最小二乘估计，假定时间序列过去未观察到序列值等于序列均值，可得到残差的有限项表达式：于是残差平方和达到最小的那组参数值为：5.模型和参数的显著性检验ARMA(p,q)模型中，使用 QLB 统计量检验残差序列的自相关性，

14、 为了克服 DW 检验的有偏性，Durbin 在 1970 年提出了修正的 Durbin h 统计量： 其中，n 为观察值序列的长度， 2 为延迟因变量系数的最小二乘估计的方差。参数的显著性检验是要检验每一个模型参数是否显著非零。若某个参数为零，模型中包含这个参数的乘积项就为零，可以简化模型。因此，该检验的是为了精简模型。原假设 H0：某未知参数 j=0；H1：j0. 可以构造出检验未知参数显著性的 t(n-m)检验统计量，其中 m 为参数的个数。6.模型优化当一个拟合模型在置信水平 下通过了检验，说明了在该置信水平下该拟合模型能有效地拟合时间序列观察值的波动。但是这种有效的拟合模型并不是惟一

15、的。如果同一个时间序列可以构造两个拟合模型，且两个模型都显著有效，那么应该选择哪个拟合模型用于统计推断呢？通常采用 AIC 和 SBC 信息准则来进行模型优化。（1）AIC 准则最小信息量准则由日本统计学家赤池弘次（Akaike）于 1973 年提出，是一种考评综合最优配置的指标，它是拟合精度和参数未知个数的加权函数：AIC=2ln(模型中极大似然函数值)+2(模型中未知参数个数)使其达到最小值的模型被认为是最优模型。（2）BIC/SBC 准则 AIC 准则的不足：若时间序列很长，相关信息就越分散，需要多自变量复杂拟合模型才能使拟合精度比较高。在 AIC 准则中拟合误差等于n ln( 2 )

16、，即随样本容量 n 增大，但模型参数个数的惩罚因子（始终=2）却与 n 无关。因此在样本容量 n 趋于无穷大时，由 AIC 准则选择的拟合模型不收敛于真实模型，它通常比真实模型所含的未知参数个数要多。为了弥补 AIC 准则的不足，Akaike 于 1976 年提出 BIC 准则。而 Schwartz 在 1978 年根据贝叶斯理论也得出同样的判别准则，称为SBC 准则。SBC 准则定义为：SBC=2ln(模型中极大似然函数值)+ln(n)(模型中未知参数个数)即将未知参数个数的惩罚权重由常数 2 变成了 ln(n)。在所有通过检验的模型中使得 AIC 或 SBC 函数达到最小的模型为相对最优模

17、型（因为不可能比较所有模型）。7.模型预测即利用时间序列已观察到的样本值对时间序列在未来某个时刻的取值进行估计。常用的预测方法是线性最小方差预测。根据 ARMA(p,q)模型的平稳性和可逆性，可以用格林函数的传递形式和逆转函数的逆转形式等价描述该序列：右式代入左式得： xt = Gi I j xt-i- j = Gi I j xt-i- j = Ci xt-1-ii=0 j=0 i=0j=0i=0可见，xt 是历史数据 xt-1, xt-2, 的线性函数。对于任意一个将来时刻 t+l，也可以用上式预测，但 xt+l-1, , xt+1 未知。根据线性函数的可加性，所有未知信息都可以用已知信息的

18、线性函数表示出来，并用该线性函数进行估计：用et (l) = xt+l - xt+l 来衡量预测误差，最常用的预测原则是预测误差的方差最小法：在线性预测方差最小法下得到的估计值 xt+l 是在序列 xt, xt-1, 已知的情况下得到的条件无偏最小方差估计值。且预测方差只与预测步长 l 有关，而与预测起始点 t 无关。预测步长 l 越大预测值的方差越大，因此只适合于短期预测。在正态假定下，估计值 xt+l 的 1- 的置信区间为：（二）ARIMA 模型混和自回归移动平均模型一、原理也称 Box-Jenkins 模型，用来处理单变量同方差的非平稳时间序列，通过差分法或适当的变换转化为平稳序列，再

19、使用 ARMA 模型。注：残差的条件方差是异方差的时间序列，适合用 GARCH 模型。ARIMA(p,d,q)模型的形式如下：(B)d x= (B) 或 d x = (B) t其中， dt t t= (I - B）d 为 d 阶差分，(B)为平稳可逆 ARMA(p,q)模型的自回归和移动平均系数多项式。可见，ARIMA 模型的实质就是差分运算与 ARMA 模型的组合。任何非平稳序列只要通过适当阶数的差分实现平稳，就可以对差分后序列进行 ARMA 模型的拟合了。dd 阶差分后的序列可表示为：其中， Ci 为组合数，即 d 阶差分后序列等于原来序列的若干序列值的某种加权和。二、建模步骤分为三个阶段

20、：识别阶段、估计阶段和预测阶段。1.识别阶段使用 identify 语句来指定响应变量序列并且识别候选 ARIMA 模型。一般先对序列进行非线性、差分和平稳性检验，可能对序列进行差分，然后计算自相关系数 ACF、逆自相关系数 IACF、偏自相关系数 PACF 和互相关系数。此阶段的输出通常会建议一个或多个可拟合的 ARIMA 模型。如果模型确定，还可以检验样本自相关系数 SACF 和样本偏自相关系数 SPACF，以分出模型的类型。2.估计阶段使用 estimate 语句来指定 ARIMA 模型去拟合在前面 identify 语句中指定的响应变量，并且估计该模型的参数。estimate 语句也生

21、成诊断统计量从而帮助判断该模型的适用性。关于参数估计值的显著性检验可以指出模型里的一些项是否不需要：拟合优度统计量 R2 可帮助比较该模型和其他模型的优劣；白噪声残差检验可指明残差序列是否包含可被其他更复杂模型采用的额外信息，如果诊断检验表明模型不适用，则可尝试另一个模型然后重复估计和诊断。3.预测阶段使用 forecast 语句来预测时间序列的未来值，并对这些来自前面estimate 语句生成的 ARIMA 模型的预测值产生置信区间。（三）PROC ARIMA 过程ARIMA 过程采用 Box-Jenkins 方法建立模型，是集一元时间序列模型判定、参数估计和预测为一体的多功能综合工具。当

22、ARIMA 模型包括其他时间序列作为输入变量时，有时也被称为 ARIMAX 模型。ARIMA 模型还支持干预或中断时间序列模型、误差的多元回归分析、任意复杂程度的有理转移函数模型。基本语法：proc arima data=数据集 out=输出数据集; where 条件表达式;identify var=变量() ; estimate ;forecast ;说明：（1）where 语句指定用于分析的时间间隔，通常条件表达式是有关日期变量的条件表达式，例如：31dec98d 日期变量 31dec99d（2）identify 语句主要完成时间序列的差分计算，样本 ACF、IACF 和 PACF 函数的

23、计算、卡方检验统计量和白噪声自相关检验的 p 值的计算。主要选项： var=变量(d1,dk)是必选项，指定要分析的时间序列变量， 按括号内列出的差分周期列表来计算时间序列的滞后差分。例如：var=X(1) 为对滞后 1 项的序列差分，即 Xt-Xt-1; var=X(2) 为对滞后 2 项的序列差分，即 Xt-Xt-2; var=X(1,1) 为 X 进行二阶差分，即(Xt-Xt-1)-(Xt-1-Xt-2); nlag=数字指定计算自相关的滞后数，其值应大于 p+d+q，小于观测数，默认值为 24； crosscorr=(干预变量(d1)列出有 var=指定的响应序列的交叉相关变量。干预变

24、量在交叉相关变量中。交叉相关变量的差分由圆括号内的差分滞后数确定。（3）estimate 语句对已执行的 identify 语句中的响应变量规定一个模型，主要选项： p=(p1,p2,)(p1,p2,)定义一个在 p 中指定的滞后处具有自回归参数的模型，p 的默认值为 0; q=(q1,q2,)(q1,q2,)定义一个在 q 中指定的滞后处具有滑动平均参数的模型，q 的默认值为 0。如果 p=和 q=都没有指定，则拟合随机模型； noconstant在模型中舍弃常数项 ; noint在该模型中不拟合截距参数； method=ml | uls | cls指定估计时使用的方法，分别为极大似然方法、

25、无条件最小二乘法、有条件最小二乘法，默认为 cls; outest=数据集将参数估计值输出到指定的数据集； outmode=数据集将模型和参数估计值输出到指定的数据集； outstat=数据集将模型诊断统计量输出到指定的数据集； plot可以绘制残差自回归函数等；（4）forecast 语句利用 estimate 语句所产生的参数估计生成时间序列的预测值，主要选项： alpha=设置预测置信限的大小，上下置信限的置信水平为 1-，默认值为 0.05; lead=n指定要计算的多步向前预测值的次数，默认为 24； back=n指定在数据末尾前 n 个观测值开始进行分步预测，默认为 0; inte

26、rval=时间间隔指定观测之间的时间间隔，常用的时间间隔为 year、qtr、month、week、weekday、day、hour、minute、second; id=变量名指明输入数据集中一个变量，用于识别与观测有关的时间周期； out=数据集将预测值和其他值输出到一个指定的数据集中。例 1 有一组 1949 年至 1961 年国际航线旅客月度人数的记录：YEAR1234567891011121949112118132129121135148148136119104118195011512614113512514917017015813311414019511451501781631721

27、781991991841621461661952171180193181183218230242209191172194195319619623623522924326427223721118020119542041882352272342643022932592292032291955242233267269270315364347312274237278195628427731731331837441340535530627130619573153013563483554224654674043473053361958340318362348363435491505404359310337

28、19593603424063964204725485594634073624051960417391419461472535622606408461390432使用 ARIMA 过程进行建模和预测。（一）读入数据、绘制时间序列图、判断平稳性，代码： data arimad01; date=intnx(month,31dec1948d,_n_);input x ;format date monyy5.; datalines;112118132129121135148148136119104118115126141135125149170170158133114140145150178163172178199199184162146166171180193181183218230242209191172194