ARIMA预测原理以及SAS实现代码

1、-ARIMA定义ARIMA的完整写法为ARIMA(p,d,q)A其中p为自回归系数，代表数据呈现周期性波动Ad为差分次数，代表数据差分几次才能到达平稳序列Aq为移动平均阶数，代表数据为平稳序列，可以用移动平均来处理。获得观察值序列分析结束差分运算拟合ARMA模型平稳性检测方法A方法一：时序图序列始终在一个常数值附近随机波动，且波动范围有界，且没有明显的趋势性或 周期性，所以可认为是平稳序列。下列图明显不是一个平稳序列proc gplot data=gdp;plot gdp*year=l ;symbol c=red i=join v=star;run;symbols c=green i=join

2、 v=none 1=2;run;4一二介仔会介介处i-r-iI1IIr1975198019851990/ /: /: /:/:/ /; :/ X/ ./ |11r- 1I-1Ir t | v1995200020052010另一种确定p、q的方式proc arima data= gdp;identify var=gdp stationarity = (adf=3) ;run;直接对gdp求arima模型，可已看出acf是拖尾，而pacf是1阶截尾，所以最好是p=l ,q=0300002500020000150001000050000gdp”的趋势和相关分析1.0-0.5b o.o -0.5-1.

3、0-1.0-0.5b o.o -0.5-1.0-0123456滞后11111110123456滞后0123456滞后确定p、q的方式理论由于/及比4(0,/模型可以转化为无穷阶移动平均模型，所以/QV勿(p,q)模型的自相关系数不截尾。同理，由于ZRK4()国)模型也可以转化为无穷阶自回归模型，所以4EH4(p,q)模型的偏自相关系数也不截尾。总结/H(p)模型、M4(q)模型和模型的自相关系数和偏自相关系数的规律，见表6.1所示。模型自相关系数偏自相关系数。尿AR(p)拖尾P阶截尾MA(q)阶截尾拖尾ARMA(p,q)拖尾拖尾模型优化指标当一个拟合模型在指定的置信水平。下通过了检验，说明了在

4、这个置信水平口下该拟合 模型能有效地拟合时间序列观察值的波动。但是这种有效的拟合模型并不是惟一的。如果同 一个时间序列可以构造两个拟合模型，且两个模型都显著有效，那么应该选择哪个拟合模型 用于统计推断呢？通常采用AIC和SBC信息准那么来进行模型优化。AIC准那么AIC准那么是由日本统计学家赤池弘次（Akaike）于1973年提出，AIC全称是最小信息量准 那么（an information criterion）。AIC准那么是一种考评综合最优配置的指标，它是拟合精度和参 数未知个数的加权函数：AIC=-21n （模型中极大似然函数值）+2 （模型中未知参数个数）（6.68）使AIC函数到达最

5、小值的模型被认为是最优模型。BIC准那么AIC准那么也有缺乏之处：如果时间序列很长，相关信息就越分散，需要多自变量复杂拟合 模型才能使拟合精度比拟高。在AIC准那么中拟合误差等于ln（3；）,即拟合误差随样本容量 放大。但是模型参数个数的惩罚因子却与无关，权重始终为常数2。因此在样本容量趋于 无穷大时，由AIC准那么选择的拟合模型不收敛于真实模型，它通常比真实模型所含的未知参 数个数要多。为了弥补AIC准那么的缺乏，Akaike于1976年提出BIC准那么。而Schwartz在1978年根据 Bays理论也得出同样的判别准那么，称为SBC准那么。SBC准那么定义为：SBC-21n（模型中极大似

6、然函数值）+ln（n）（模型中未知参数个数）（6.69）它对AIC的改进就是将未知参数个数的惩罚权重由常数2变成了样本容量n的对数ln（n）。在 所有通过检验的模型中使得AIC或SBC函数到达最小的模型为相对最优模型。之所以称为相 对最优模型是因为不可能比拟所有模型。表6. 2河南省历年国民生产总值数据附：完整代码年份(Year)生产总值（亿元）（GDP）人均生产总值（元）（PGDP）年份(Year)生产总值（亿元）（GDP）人均生产总值（元）（PGDP）1978162.92232.319921279.751452.31979190.09266.719931662.761867.4198022

7、9.16316.719942224.432475.21981249.69340.119953002.743312.81982263.3035319963661.184007.41983327.95432.919974079.264430.11984370.04481.619984356.604695.11985451.74579.719994576.104893.71986502.91635.320005137.6654441987609.60755.820015640.115923.61988749.09909.920026168.736436.51989850.711012.3200370

8、48.597570.21990934.651090.620048815.099469.919911045.731201.2data gdp;infile datalines;input year gdp pgdp;format gdp BEST12. 2 pgdp BEST12. 2; datalines;1991 1045. 73 1201.2 1992 1279. 75 1452. 3 1993 1662. 76 1867.4 1994 2224.43 2475.2 1995 3002. 74 3312.8 1996 3661.18 4007.4 1997 4079.26 4430.1 1

9、998 4356.60 4695.1 1999 4576. 10 4893. 7 2000 5137.66 5444 2001 5640. 11 5923. 6 2002 6168.73 6436.5 2003 7048.59 7570.2 2004 8815.09 9469.91978162. 92232.31979190.09266. 71980229. 16316.71981249. 69340. 11982263. 303531983327. 95432.91984370.04481.61985451. 74579.71986502.91635.31987609. 60755.8198

10、8749. 09909.91989850. 711012. 31990934. 651090.6run;/*原始数据散点图*/proc gplot data=gdp;plot gdpyear=l ;symbol c=red i=join v=star;run;/*注symbol常用参数C图形颜色，red红色，black黑色，green绿色，blue蓝色，pink洋红等*/ /* V一观测值的图形，star *, dot_. , cicle_圆圈，diamond 菱形，none 不标 */ /* I一观察值的链接方式，join_线连，spline_光滑连接，needle_作观察值到横轴悬 垂线，

11、none_不连*/proc arima data= gdp;identify var=gdp stationarity =(adf=3);run;/*原始数据对数、差分变换*/ data gdplog;set gdp;loggdp=log(gdp);cfloggdp=dif(loggdp);run;/*对数数据散点图*/proc gplot;plot loggdp*year=l ;symbol c=black i=join v=star;run;/* 一阶差分对数数据散点图*/proc gplot;plot cfloggdp*year=l;symbol c=green v=dot i=join

12、;run;/* 一阶差分对数数据的自相关图、偏自相关图、纯随机性检验、单位根检验*/ proc arima data=gdp_log;identify var=loggdp(1) stationarity = (adf=3) nlag=12;run;/* loggdp(l)这里的数1为差分阶数*/*定阶*/proc arima data=gdp_log;identify var=loggdp(1) nlag=6 minic p=(0:2) q=(0:4);run;/* minic为一定范围模型定阶相对最优模型识别*/*参数估计*/ proc arima data=gdplog; identif

13、y var=loggdp(1); estimate p=l q=0 0UTEST=b outstat=c;run;/* SAS支持三种估计，默认为条件最小二乘估计，要制定可增加选项：METH0D=ML极大似然估计METHOD=ULS最小二乘估计METHOD=CLS条件最小二乘估计 输出项的含义见王燕P104V；/ *参数估计及预测*/ proc arima data=gdp_log; identify var=loggdp(1) nlag=16; estimate p=l q=0;forecast lead=4 id=year out=results; run;/*绘制预测图*/data re

14、sults;set results;y=exp(loggdp);estimatel=exp(forecast);e195=exp (195);eu95=exp (u95);run;proc gplot data=results;plot y*year=l estimatelyear=2 el95*year=3 eu95*year=3/overlay;symbol1 c=black i=none v=star;symbo12 c=red i=join v=none;symbols c=green i=join v=none 1=2;run;fldp100008000600040002000 0p

15、roc arima data= gdp;identify var=gdp stationarity =(adf=3) nlag=12;run;A ADF单位根检验（精确判断）三个检验中只要有一个PrRho小于0.05即可认定为平稳序列，主要是stationarity=(adf=3)起作用proc arima data= gdp;identify var=gdp stationarity = (adf=3) nlag=12;run;增广Dickev-Ful2r单位根检羲类型滞后RholPr RhoTauPr 卡方丰方自相关667.426.00010.8370.7160.6100.5060.404

16、0.3121273.2312.00010.2140.1120.010-0.080-0.148-0.201非平稳序列转换为平稳序列方法一：将数据取对数。方法二：对数据取差分dif函数data gdp_log;set gdp;loggdp=log(gdp);cfloggdp=dif(loggdp);run;/*对数数据散点图*/proc gplot;plot loggdp*year=l ;symbol c=black i=join v=star;run;loggdp1 0 T-/* 一阶差分对数数据散点图*/proc gplot;plot cfloggdp*year=l;symbol c=gree

17、n v二dot i二join;run;cfloggdp I 1 I 厂 * L I * I 1下面代码中的(1)就代表1阶差分，adf=3那么代表平稳性检验0-3, /* 一阶差分对数数据的自相关图、偏自相关图、纯随机性检验、单位根检验*/ proc arima data=gdp log;identify var=loggdp(1) stationarity =(adf=3) nlag=12;run;白噪声的目相关侬,同卡周苜重度I Pfjp密I自相关一6 18.816 COW40O,479 0.16TT0,161 -0.293 -0.309 -0.34712 20.45121011550 -

18、0.139 0.044 0.107 0.068 -0.017 -0.046增广Dickey.Fuller单位根检验类型滞后RhoPr RhoTauPr FZero MeanM .8396|0.3437-0.850.33791-1.54520.3815)po.790.36492pt 06100.4567-0.520.47933p0.44950.5685-0.286.5744Single Meanfl-12.53320.0434-2.700.08773.670.20451-15.19870.0166g.533.20031122-39.7723一.0001is3/M575)M.370.08813-1

19、58.1750.000129。4.340.0897Trend0-12.48910.1999-2.610.27803.490.52281p4.97180.1004p.420.35983.080.59442-40.3809r.oooil0.19204.090.4171|3-166.0470.0001-2.910.17864.300.3817用Qlb统计量作的2检验结果说明：对数差分后的GDP序列的Qlb统计量的。值为0.0045 (卡方目相关6 18.810.0045 0.479 0.161 -0.161 -0.293 -0.309 -0.347Minimum Information Criter

20、ionLagsMAO, MA1| MA 2MA 3MA 4AR0-5.50235-5.92658-5.91546-6.03307-6.03373ARi|-5.9469-5.87077-6.01308-6.36963厂6.376AR 2-6.20883-6.09519-5.97303-6.3134-6.28134误差的麒型：AR(6)最小表值：BIC(1,4)= 6376 采用相对最优模型识别，根据上述分析及序列的自相关和偏自相关图，适中选择?二4, 二 2,使用indentify命令中的minic p = (0: n) q = (0: m)短语进行相对最优模型定阶。结果显示(图6.10),在P

21、=1, 0=4时，BIC函数值最小。执行ARIMA过程的Estimate p = l q =4命令做参数检验，结果未能通过参数检验。让q在。3之间取值，通过反复测 试，只有ARMA(1,3)模型与ARMA(1, 0)模型通过参数检验及模型检验，其检验结果及参数估 计如图6.11所示。参数估计/*参数估计*/proc arima data=gdp_log; identify var=loggdp (1); estimate p=l q=4;/* SAS支持三种估计，默认为条件最小二乘估计，要制定可增加选项:METHOD=ML极大似然估计METHOD=ULS最小二乘估计METHOD二CLS条件最小

22、二乘估计输出项的含义见王燕P104V；条件最小二乘估讦参数估计标准误差t值近似Pr |t|滞后MU3.610.0017MA1,10.42352| 5.231530.080.93631MA1,20.03084 2.941690.010.9917MA1,30.43527p 2.787040.160.87753MA1940.11026| 0,551330.200.84354AR1,10.98559 5.117540.190.8492从上面可以看出，在p=l q=4时,通不过检验。 条件最小二乘估计参数I估计标准误差t值近似Pr |t|滞后0.163580.02679F6.11 |t|滞后MU0.15

23、5960.023596.61.00010amT0.498530.181732.740.01131；常18祜计二f 0.078207：:方差估计 * 0.004036： :标准误差估计0,06353：AIC-p67.6208 ：SBC| -65.1046：滋差数26；AIC和SBC不包括对数行列式。P=1 q=0时能通过检验从上面2个模型的检验结果可以看到，它们均为有效模型，但ARMA(1,O)模型的AIC为-67, SBC为-65均比ARMA(1, 3)的AIC与SBC小,根据AIC准那么和SBC准那么，ARMA(1, 0)应该更有 效，所以应选择前者作为预测模型。AIU 机 个臼他 E ST

24、U 工、。参数估计相关性 爹数工MU AR1J MU JT.000 0.069 AR1,1 0.069 1.000残差的目相关检查至滞后卡方自由度Pr 卡方目相关1 6|6.0750.29970.0540.062-0.189-0.191-0.081-0.3017.221107806-0.0540.0650.0880.080-0.0140.0681811.99170.8006-0.087-0.1890.0390.0450.144-0.0522415.94230.85770.0000.0260.103-0.075-0.0660.026:变蚩可丽前片的模型：:估计均值 0.155955：:差分期间1

25、；自回归因因子 1:jV 0.49853GDP对数序列模型的口径为：rVlogx, =0.155955 +1-0.498538其中，为表示GDP序列，模型可写为：log xt = 1.49853 log -0.49853 logxz_2 + 0.078207 + st预测/*参数估计及预测*/proc arima data=gdp_log;identify var=loggdp(1) nlag=16;estimate p=l q=0;forecast lead=4 id=year out=results;run;year loggdp 尬 FORECAST 顾STD/)L95顾吧5顾 RESI

26、DUAL119785.0932592828:219795.24749764425.24921470750.06353023575.12469773365.3737316814-0.001717063319805.43442044955.40259706460.06353023575.27808009075.52711403850.0318233854*19815.52022014845.60581398460.06353023575.48129701075.7303309585-0.085593836119825.57329406655.64120097030.06353023575.5166

27、8399635.7657179442-0.067906904619835.79286115775.6779601660.06353023575.55344319215.80247713990. 1149009918T19845.91361110795.98052878770.06353023575.85601181376.1050457616-0.06691768819856.11310679316.05201561320.06353023575.927498639360610911799919866.22041122766.29076826250.063530235

28、76.16625128866.4152852364-0.0703570351019876.41280300446.35211276520.06353023576.22759579126.47662973910.06069023921119886.61885913656.58692297550.06353023576.46240600166.71143994940.031936161:1219896.74607129496.79979117320.06353023576.67527419936.9243081471-0.053719878；1319906.84017212776.88769739

29、280.06353023576.76318041897.0122143667-0.0475252651419916.95247048516.96529129820.06353023576.84077432437.0898082721-0.0128208131151992708666163370.06353023576.96214465987.21117860760.06775839171619937.41623415137.33330481050.06353023577.20878783667.45782178440.08292934081719947.7072559

30、827.62496311640.06353023577.50044614257.74948009030.1819958.00728048417.93054580960.06353023577.80602883578.05506278350.07673467451919968.20554077888.23505839670.06353023578.11054142288.3595753706-0.029517618 ?2019978.31367087848.38258637120.06353023578.25806939738.5071033452-0.068915493212219988.37944721