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

█ARIMA定义

ARIMA的完整写法为ARIMA(p,d,q)

?其中p为自回归系数，代表数据呈现周期性波动 ?d为差分次数，代表数据差分几次才能达到平稳序列

?q为移动平均阶数，代表数据为平稳序列，可以用移动平均来处理。

█平稳性检测方法

?

?方法一：时序图

序列始终在一个常数值附近随机波动，且波动范围有界，且没有明显的趋势性或周期性，所以可认为是平稳序列。下图明显不是一个平稳序列

proc gplot data=gdp; plot gdp*year=1 ;

symbol c=red i=join v=star; run;

?

?方法二：自相关图

自相关系数会很快衰减向0，所以可认为是平稳序列。

proc arima data= gdp;

identify var=gdp stationarity =(adf=3) nlag=12; run;

?

?ADF单位根检验（精确判断）

三个检验中只要有一个Pr

proc arima data= gdp;

identify var=gdp stationarity =(adf=3) nlag=12; run;

█白噪声检验

Pr>卡方<0.05即可认定为通过白噪声检验。 proc arima data= gdp;

identify var=gdp stationarity =(adf=3) nlag=12; run;

█非平稳序列转换为平稳序列

方法一：将数据取对数。

方法二：对数据取差分dif函数

data gdp_log; set gdp;

loggdp=log(gdp);

cfloggdp=dif(loggdp); run;

/**对数数据散点图**/ proc gplot;

plot loggdp*year=1 ;

symbol c=black i=join v=star; run;

/* 一阶差分对数数据散点图*/ proc gplot;

plot cfloggdp*year=1;

symbol c=green v=dot i=join;

run;

从上图中可以看出，一阶差分后序列已经变成平稳的了，因此，数列需要做一阶差分

█转换完毕后再验证

下面代码中的（1）就代表1阶差分，adf=3则代表平稳性检验0-3，

/* 一阶差分对数数据的自相关图、偏自相关图、纯随机性检验、单位根检验 */ proc arima data=gdp_log;

identify var=loggdp(1) stationarity =(adf=3) nlag=12;

run;

用QLB统计量作的?2检验结果表明：对数差分后的GDP序列的QLB统计量的P值为0.0045（<0.05），故序列为非白噪声序列。

单位根检验结果表明：对数差分后的GDP序列有常数均值、无趋势的二阶自回归模型的Tau统计量的P值小于0.0573，故序列基本可以确定为平稳序列，并可初步考虑用ARMA（2，q）模型对它们进行拟合。

█模型定阶

/** 定阶 **/

proc arima data=gdp_log;

identify var=loggdp(1) nlag=6 minic p=(0:2) q=(0:4); run;

/* minic为一定范围模型定阶——相对最优模型识别 */

采用相对最优模型识别，根据上述分析及序列的自相关和偏自相关图，适当选择m = 4，n = 2，使用indentify命令中的minic p = (0: n) q = (0: m)短语进行相对最优模型定阶。

结果显示（图6.10），在p = 1，q = 4时，BIC函数值最小。执行ARIMA过程的Estimate p = 1 q = 4命令做参数检验，结果未能通过参数检验。让q在0～3之间取值，通过反复测试，只有ARMA(1, 3)模型与ARMA(1, 0)模型通过参数检验及模型检验，其检验结果及参数估计如图6.11所示。

█参数估计

/** 参数估计 **/ proc arima data=gdp_log; identify var=loggdp(1); estimate p=1 q=4; run;

/* SAS支持三种估计，默认为条件最小二乘估计，要制定可增加选项： METHOD=ML 极大似然估计 METHOD=ULS 最小二乘估计

METHOD=CLS 条件最小二乘估计

输出项的含义见王燕 P104*/;

从上面可以看出，在p=1 q=4时，通不过检验。

p=1 q=3 和

p=1 q=0时能通过检验

从上面2个模型的检验结果可以看到，它们均为有效模型，但ARMA(1, 0)模型的AIC为-67，SBC为-65均比ARMA(1, 3)的AIC与SBC小，根据AIC准则和SBC准则，ARMA(1, 0)应该更有效，所以应选择前者作为预测模型。

GDP对数序列模型的口径为：

?logxt?0.155955?其中，xt表示GDP序列，模型可写为：

?t1?0.49853B

logxt?1.49853logxt?1?0.49853logxt?2?0.078207??t

█预测

/**参数估计及预测**/ proc arima data=gdp_log;

identify var=loggdp(1) nlag=16; estimate p=1 q=0;

forecast lead=4 id=year out=results;

run;

█还原预测数值并画图

/**绘制预测图**/ data results; set results; y=exp(loggdp);

estimate1=exp(forecast); el95=exp(l95); eu95=exp(u95); run;

proc gplot data=results;

plot y*year=1 estimate1*year=2 el95*year=3 eu95*year=3/overlay; symbol1 c=black i=none v=star; symbol2 c=red i=join v=none;

symbol3 c=green i=join v=none l=2;

run;

█另一种确定p、q的方式

proc arima data= gdp;

identify var=gdp stationarity =(adf=3) ;

run;

直接对gdp求arima模型，可已看出acf是拖尾，而pacf是1阶截尾，所以最好是p=1，q=0

█确定p、q的方式理论

由于ARMA(p,q)模型可以转化为无穷阶移动平均模型，所以ARMA(p,q)模型的自相关系数不截尾。同理，由于ARMA(p,q)模型也可以转化为无穷阶自回归模型，所以

ARMA(p,q)模型的偏自相关系数也不截尾。总结AR(p)模型、MA(q)模型和ARMA(p,q)模型的自相关系数和偏自相关系数的规律，见表6.1所示。 模型 自相关系数?k 拖尾 偏自相关系数?kk AR(p) MA(q) ARMA(p,q) p阶截尾 拖尾 拖尾 q阶截尾 拖尾

█模型优化指标

当一个拟合模型在指定的置信水平?下通过了检验，说明了在这个置信水平?下该拟合模型能有效地拟合时间序列观察值的波动。但是这种有效的拟合模型并不是惟一的。如果同一个时间序列可以构造两个拟合模型，且两个模型都显著有效，那么应该选择哪个拟合模型用于统计推断呢？通常采用AIC和SBC信息准则来进行模型优化。

1. AIC准则

AIC准则是由日本统计学家赤池弘次（Akaike）于1973年提出，AIC全称是最小信息量准则（an information criterion）。AIC准则是一种考评综合最优配置的指标，它是拟合精度和参数未知个数的加权函数：

(6.68) AIC=－2ln(模型中极大似然函数值)+2(模型中未知参数个数) 使AIC函数达到最小值的模型被认为是最优模型。

2. BIC准则

AIC准则也有不足之处：如果时间序列很长，相关信息就越分散，需要多自变量复杂拟合

??2)，模型才能使拟合精度比较高。在AIC准则中拟合误差等于nln(?即拟合误差随样本容量n放大。但是模型参数个数的惩罚因子却与n无关，权重始终为常数2。因此在样本容量n趋于

无穷大时，由AIC准则选择的拟合模型不收敛于真实模型，它通常比真实模型所含的未知参数个数要多。

为了弥补AIC准则的不足，Akaike于1976年提出BIC准则。而Schwartz在1978年根据Bays理论也得出同样的判别准则，称为SBC准则。SBC准则定义为：

(6.69) SBC=－2ln(模型中极大似然函数值)+ln(n)(模型中未知参数个数) 它对AIC的改进就是将未知参数个数的惩罚权重由常数2变成了样本容量n的对数ln(n)。在所有通过检验的模型中使得AIC或SBC函数达到最小的模型为相对最优模型。之所以称为相

对最优模型是因为不可能比较所有模型。

表6.2 河南省历年国民生产总值数据

年份(Year) 1978 1979 1980 1981 1982 1983 1984 1985 1986 1987 1988 1989 1990 1991 生产总值(亿元) (GDP) 162.92 190.09 229.16 249.69 263.30 327.95 370.04 451.74 502.91 609.60 749.09 850.71 934.65 1045.73 人均生产总值(元)(PGDP) 232.3 266.7 316.7 340.1 353 432.9 481.6 579.7 635.3 755.8 909.9 1012.3 1090.6 1201.2 年份(Year) 1992 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 生产总值(亿元) (GDP) 1279.75 1662.76 2224.43 3002.74 3661.18 4079.26 4356.60 4576.10 5137.66 5640.11 6168.73 7048.59 8815.09 人均生产总值(元)(PGDP) 1452.3 1867.4 2475.2 3312.8 4007.4 4430.1 4695.1 4893.7 5444 5923.6 6436.5 7570.2 9469.9

█附：完整代码

data gdp;

infile datalines; input year gdp pgdp;

format gdp BEST12.2 pgdp BEST12.2; datalines;

1978 162.92 232.3 1979 190.09 266.7 1980 229.16 316.7 1981 249.69 340.1 1982 263.30 353 1983 327.95 432.9 1984 370.04 481.6 1985 451.74 579.7 1986 502.91 635.3 1987 609.60 755.8 1988 749.09 909.9 1989 850.71 1012.3 1990 934.65 1090.6 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 1998 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.9 ; run;

/** 原始数据散点图 **/ proc gplot data=gdp; plot gdp*year=1 ;

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_作观察值到横轴悬垂线，none_不连 */

proc arima data= gdp;

identify var=gdp stationarity =(adf=3) ; run;

/** 原始数据对数、差分变换 **/ data gdp_log; set gdp;

loggdp=log(gdp);

cfloggdp=dif(loggdp); run;

/**对数数据散点图**/ proc gplot;

plot loggdp*year=1 ;

symbol c=black i=join v=star; run;

/* 一阶差分对数数据散点图*/ proc gplot;

plot cfloggdp*year=1;

symbol c=green v=dot i=join; run;

/* 一阶差分对数数据的自相关图、偏自相关图、纯随机性检验、单位根检验 */ proc arima data=gdp_log;

identify var=loggdp(1) stationarity =(adf=3) nlag=12; run;

/* loggdp(1)这里的数1为差分阶数 */

/** 定阶 **/

proc arima data=gdp_log;

identify var=loggdp(1) nlag=6 minic p=(0:2) q=(0:4); run;

/* minic为一定范围模型定阶——相对最优模型识别 */

/** 参数估计 **/ proc arima data=gdp_log; identify var=loggdp(1);

estimate p=1 q=0 OUTEST=b outstat=c; run;

/* SAS支持三种估计，默认为条件最小二乘估计，要制定可增加选项： METHOD=ML 极大似然估计 METHOD=ULS 最小二乘估计

METHOD=CLS 条件最小二乘估计 输出项的含义见王燕 P104*/;

/**参数估计及预测**/ proc arima data=gdp_log;

identify var=loggdp(1) nlag=16; estimate p=1 q=0;

forecast lead=4 id=year out=results; run;

/**绘制预测图**/ data results; set results; y=exp(loggdp);

estimate1=exp(forecast); el95=exp(l95); eu95=exp(u95); run;

proc gplot data=results;

plot y*year=1 estimate1*year=2 el95*year=3 eu95*year=3/overlay; symbol1 c=black i=none v=star; symbol2 c=red i=join v=none;

symbol3 c=green i=join v=none l=2; run;

本文来源：https://www.bwwdw.com/article/aqgr.html