王燕时间序列分析第五章SAS程序

第一题
data yx_51;
input x@@;
difx=dif(x);
t=1+_n_-1;
cards;
304 303 307 299 296 293 301 293 301 295 284 286 286 287 284 282 278 281 278 277 279 278 270 268 272 273 279 279 280 275 271 277 278 279 283 284 282 283 279 280 280 279 278 283 278 270 275 273 273 272 275 273 273 272 273 272 273 271 272 271
273 277 274 274 272 280 282 292 295 295 294 290 291 288 288 290 293 288 289 291 293 293 290 288 287 289 292 288 288 285 282 286 286 287 284 283 286 282 287 286 287 292 292 294 291 288 289
;
proc gplot;
plot x*t=1 difx*t=2;
symbol1c=red v=circle i=join;
symbol2c=yellow v=star i=join;
run;
proc arima;
identify var=x(1);
estimate p=1;
run;
结果如下
时序图：
一阶差分后时序图:
通过原始数据的时序图可以明显看出，此序列非平稳，因而对序列进行一阶差分。

从一阶差分后的自相关图可以看出，一阶差分后的序列的自相关系数一直都比较小，始终控制在二倍标准差以内，可以认为一阶差分后的序列始终都在零轴附近波动，因而可以认为一阶差分后的序列为随机性很强的平稳序列，另外通过一阶差分后的时序图也可以看出，一阶差分后的序列平稳，且LB统计量对应
的P值大于α=0.05，因而认为一阶差分后的序列为白噪声序列。

由于一阶差分后的序列为平稳的白噪声序列，因而此时间序列拟合ARIMA （0,1,0）模型，即随机游走模型，模型为：
X t=x t-1+εt
所以下一期的预测值为289
第二题
data yx_52;
input x@@;
t=1949+_n_-1;
difx=dif(x);
cards;
5589.00 9983.00 11083.00 13217.00 16131.00 19288.00 19376.00 24605.00
27421.00 38109.00 54410.00 67219.00 44988.00 35261.00 36418.00 41786.00
49100.00 54951.00 43089.00 42095.00 53120.00 68132.00 76471.00 80873.00
83111.00 78772.00 88955.00 84066.00 95309.00 110119.00 111893.00 111279.00
107673.00 113495.00 118784.00 124074.00 130709.00 135635.00 140653.00 144948.00
151489.00 150681.00 152893.00 157627.00 162794.00 163216.00 165982.00 171024.00
172149.00 164309.00 167554.00 178581.00 193189.00 204956.00 224248.00 249017.00
269296.00 288224.00 314237.00 330354.00
;
proc gplot;
plot x*t=1 difx*t=2;
symbol1c=orange v=circle i=none;
symbol2c=blue v=star i=join;
proc arima;
identify var=x(1);
estimate q=1;
forecast lead=5id=t;
run;
从时序图可以看出，时间序列非平稳，且随着时间而呈现明显的上升趋势，因而对序列采用一阶差分：
通过原始数据的时序图可以明显看出，此序列非平稳，随着时间呈现上升趋势，因而对序列进行一阶差分。

从一阶差分后的自相关图可以看出，一阶差分后的序列的自相关系数一阶截尾，拟合ARIMA(0，1，1)模型，得到模型：
X t-X t-1=(1+0.48349B) εt
残差的检验显示，残差序列通过白噪声检验，参数显著性检验显示参数显著，说明模型拟合良好，对序列相关信息提取充分。

得到2009~2013年铁路货运量的预测结果如下：
第三题;
data yx_53;
input x@@;
difx=dif(dif12(x));
t=intnx('month','01jan1973'd,_n_-1);
format t date.;
cards;
9007.00 8106.00 8928.00 9137.00 10017.00 10826.00 11317.00 10744.00
9713.00 9938.00 9161.00 8927.00 7750.00 6981.00 8038.00 8422.00 8714.00 9512.00
10120.00 9823.00 8743.00 9129.00 8710.00 8680.00 8162.00 7306.00 8124.00
7870.00 9387.00 9556.00 10093.00 9620.00 8285.00 8433.00 8160.00 8034.00
7717.00 7461.00 7776.00 7925.00 8634.00 8945.00 10078.00 9179.00 8037.00
8488.00 7874.00 8647.00 7792.00 6957.00 7726.00 8106.00 8890.00 9299.00 10625.00
9302.00 8314.00 8850.00 8265.00 8796.00 7836.00 6892.00 7791.00 8129.00 9115.00
9434.00 10484.00 9827.00 9110.00 9070.00 8633.00 9240.00
;
proc gplot;
plot x*t=1 difx*t=2;
symbol1c=coral v=circle i=join;
symbol2c=blue v=star i=join;
run;
proc arima;
identify var=x(1,12);
estimate p=1q=(1)(12) ;
run;
从时序可以看出，序列呈现出周期性的变化趋势，且序列非平稳，因而对序列进行一阶12步差分。

一阶12步差分后的时序图显示，自相关系数在延迟12阶时显著大于2倍标准差范围，说明差分后的序列仍然含有显著的季节效应，因而考虑拟合乘积季节模型;
根据自相关系数在延迟1阶和延迟12阶时显著大于2倍标准差范围，偏自相关系数延迟1阶显著大于2倍标准差范围，考虑拟合ARIMA(1,1,1)×(0,1,1)12。

得到模型：
∇∇12x t=(1 − 0.87376 B)
(1 − 0.49078 B)
( 1 − 0.53808 B12)εt
残差的检验显示，残差序列通过白噪声检验，参数显著性检验显示参数显著，说明模型拟合良好，对序列相关信息提取充分。