时间序列分析SAS 非平稳序列1 旅游收入

绘出原序列图如下，观察可见其有非常明显的非平稳特征：
status 30000
20000
10000
0 60Q1 62Q1 64Q1 66Q1 68Q1 70Q1 72Q1 74Q1 76Q1 78Q1 80Q1 82Q1 84Q1 t 86Q1 88Q1 90Q1 92Q1 94Q1 96Q1 98Q1 00Q1 02Q1 04Q1 06Q1 08Q1 10Q1
参数检验结果显示各参数均显著。整个模型的可决系数 R 2 高达 0.9999，正态性检验的 P 值为 0.2051， 说明模型通过了正态性检验。认为模型拟合成功。
ln xt = 1.003ln xt −1 + ε t i ,i , d 最终模型口径为 ε t = ht et , et ~ N (0, 0.0036) h = 0.00172 + 0.4728ε 2 t t −1
5 60Q1 62Q1 64Q1 66Q1 68Q1 70Q1 72Q1 74Q1 76Q1 78Q1 80Q1 82Q1 84Q1 t 86Q1 88Q1 90Q1 92Q1 94Q1 96Q1 98Q1 00Q1 02Q1 04Q1 06Q1 08Q1 10Q1
作为一个比较，最后尝试一下 Auto-Regressive 延迟因变量回归模型（即不对原序列取自然对数，直接 试拟合延迟应变量自回归模型） 。 在 data 步中加入 lagx=lag(status); 再添加命令行 proc autoreg data=example;
model logx=lagx/lagdep=lagx nlag=1 noint garch=(q=1);
输出的最终拟合结果如下：
拟合结果显示残差序列延迟 1 阶的自相关项系数（AR1）不显著。而 GARCH1 项系数也不显著，故将命 令行改写为 model logx=lagx/lagdep=lagx noint garch=(q=1);即尝试拟合 ARCH(1)模型。 模型最终拟合结果如下：
model logx=lagx/lagdep=lagx nlag=25 backstep dwprob archtest; output out=out p=xp;
尝试建立带有延迟因变量的回归模型 ln xt = a + b ln xt −1 + ut 。 输出对数序列图如下：
logx 11
10
9
8
参数估计显著。模型的表达式为xt = 1.0117 xt −1 + ut , ut
proc arima; identify var=status(1); estimate p=0 q=0; forecast lead=5 id=t out=results;
常数项拟合结果可以认为是显著的，对残差的白噪声检验结果也通过（可以认为对差分序列拟合随机游
xt = 133.27179 + xt −1 + ε t 走模型，其残差序列为白噪声）。模型的具体形式为 E (ε t ) = 0, Var (ε t ) = σ ε2 , E (ε t ε s ) = 0, s ≠ t Ex ε = 0, ∀s < t s t
长期趋势图 trend：
TREND 30000
季节因子 season：
SEASON 104
103
20000
102
101
10000
100
99
0 60Q1 62Q1 64Q1 66Q1 68Q1 70Q1 72Q1 74Q1 76Q1 78Q1 80Q1 82Q1 84Q1 t 86Q1 88Q1 90Q1 92Q1 94Q1 96Q1 98Q1 00Q1 02Q1 04Q1 06Q1 08Q1 10Q1
参数估计结果显示各参数均显著。如下：
观察残差序列自相关图。由于样本容量较大，所以输出延迟 25 阶的自相关图。观察可见残差序列具有短 期的自相关性。
同时 Q 和 LM 统计量检验结果显示残差序列具有显著的异方差性。如下：
综合考虑残差序列的自相关性和异方差性检验结果，尝试拟合不含截距项的 AR(1)—GARCH(1,1)模型。 加入命令行：
model status=lagx/lagdep=lagx; output out=out p=xp;
指令系统建立带有延迟因变量的回归模型 xt = a + bxt −1 + ut .输出结果如下：
由 Durbin h 统计量的 P 值为 0.3479，这显示残差序列不存在显著的相关性，不需要对残差序列继续拟 合自回归模型。 注意到参数检验结果显示截距项不显著，故在命令行中加入 noint。输出拟合结果为：
作出的拟合图线如下：
status 30000
20000
10000
0 60Q1 62Q3 65Q2 68Q1 70Q4 73Q3 76Q2 79Q1 81Q4 84Q3 t 87Q2 90Q1 92Q4 95Q3 98Q2 01Q1 03Q4 06Q3 09Q2
对未来十期的预测结果如下：
接下来考虑对残差信息的充分提取，同时要考虑残差是否有异方差性，所以考虑对原序列拟合 GARCH 模型。尝试拟合完成后，GARCH 模型的正态性检验不能通过。故考虑对原序列取自然对数再进行分析。 在 data 步加入命令行： logx=LOG(status); lagx=lag(logx)； 写命令行： proc autoreg data=example;
7
6
5 60Q1 62Q1 64Q1 66Q1 68Q1 70Q1 72Q1 74Q1 76Q1 78Q1 80Q1 82Q1 84Q1 t 86Q1 88Q1 90Q1 92Q1 94Q1 96Q1 98Q1 00Q1 02Q1 04Q1 06Q1 08Q1 10Q1
DW 检验结果显示残差序列具有比较显著的正自相关性。如下：
110
100
90
80 60Q1 62Q1 64Q1 66Q1 68Q1 70Q1 72Q1 74Q1 76Q1 78Q1 80Q1 82Q1 84Q1 t 86Q1 88Q1 90Q1 92Q1 94Q1 96Q1 98Q1 00Q1 02Q1 04Q1 06Q1 08Q1 10Q1
最后输出的原序列点和拟合图线图：
STATUS 30000
20000
10000
0 60Q1 62Q1 64Q1 66Q1 68Q1 70Q1 72Q1 74Q1 76Q1 78Q1 80Q1 82Q1 84Q1 t 86Q1 88Q1 90Q1 92Q1 94Q1 96Q1 98Q1 00Q1 02Q1 04Q1 06Q1 08Q1 10Q1
接下来进入主要部分：随机分析。先考虑 ARIMA 模型，对原序列进行一阶差分。在 input 语句后加入
difx=dif(status);
加入命令行：proc arima;
identify var=status(1);
一阶差分图：
difx 3000
2000
1000
0
-1000
-2000
-3000
-4000 60Q1 62Q1 64Q1 66Q1 68Q1 70Q1 72Q1 74Q1 76Q1 78Q1 80Q1 82Q1 84Q1 t 86Q1 88Q1 90Q1 92Q1 94Q1 96Q1 98Q1 00Q1 02Q1 04Q1 06Q1 08Q1 10Q1
98 60Q1 62Q1 64Q1 66Q1 68Q1 70Q1 72Q1 74Q1 76Q1 78Q1 80Q1 82Q1 84Q1 t 86Q1 88Q1 90Q1 92Q1 94Q1 96Q1 98Q1 00Q1 02Q1 04Q1 06Q1 08Q1 10Q1
随机波动图 irr：
IRR 120
编写输出命令行如下： proc gplot data=out;
plot logx*t=2 xp*t=3/overlay; symbol2 v=star i=nnone c=black; symbol3 v=none i=join c=red; run;
输出最终拟合图线如下：
logx 11
10
9
8wenku.baidu.com
7
6
一阶差分序列的白噪声检验结果为：
由差分序列延迟 6、12、18、24 阶的QLB统计量的 P 值都很大（以延迟 6 阶为例，P 值为 0.6883，说 明可以以 68.83%的置信水平，认为差分序列为白噪声序列） ，可以认为原序列一阶差分序列是白噪声序列。 因此这是一个典型的随机游走模型。尝试拟合ARIMA（0,1,0）模型，改写命令行为：
美国 1960 年第一季度到 2008 年第四季度每季度旅游收入序列 朱洵 PB08204039
在进行非平稳序列的随机性分析之前，先做个准备工作——用比较完备的 X-11 过程进行确定性分析。 编写命令行如下：
data example; input status@@; t=intnx('quarter','1jan1960'd,_n_-1); format t yyq4.; cards; 224 227 237 231 228 237 240 242 247 250 231 229 245 253 253 264 292 297 303 315 319 337 347 377 374 394 410 412 412 392 414 428 433 431 451 460 502 514 514 513 574 579 575 603 623 609 624 678 673 695 713 736 836 817 860 899 972 969 1021 1070 1187 1127 1144 1239 1324 1443 1518 1457 1437 1521 1553 1639 1673 1788 1839 1883 2044 2075 2058 2264 2547 2534 2737 2770 3108 3102 3228 3475 3280 3178 3039 2896 2647 2662 2999 2639 3860 4552 4319 4446 4363 4604 4391 4406 5004 4820 5250 5313 5340 5774 6096 6353 6711 7138 7722 7863 8535 8679 9220 9770 1043510543 10489 11539 10529 12064 12533 13260 13619 13716 13562 13845 14204 14469 14486 14716 14399 14714 14493 14810 14863 15041 16357 17133 16814 17722 16194 19079 18351 18276 18595 18204 17926 18181 17296 17922 18121 18346 18818 19516 20486 21013 20343 20558 20735 19619 17374 14165 16158 16155 16287 18005 15809 14336 16180 18034 17863 18667 18703 19313 19930 20926 20411 20532 20906 21284 21475 22055 21818 23154 25241 26499 27200 28193 28864 26212 ; proc x11 data=example; quarterly date=t; var status; output out=out b1=status d10=season d11=adjusted d12=trend d13=irr; data out; set out; estimate=trend*season/100; proc gplot data=out; plot status*t=1 estimate*t=2/overlay; plot adjusted*t=3 season*t=3 trend*t=3 irr*t=3; symbol1 c=blue v=star i=none; symbol2 c=red v=none i=spline; symbol3 c=blue v=star i=spline; run;