学院名称理学院
专业Iff级 应用统廿学14-2
学号1
齐鲁工业大学实验报告应麵
课程名称《应用时间序列分折实酚》指导教帥董实验日期2017630
•可修編.
院(系)理学院专业册级境计14 - 2 实醴地点机电楼C428
学生XX艳雪学号 _____ 同组人无
实騎顶目名称ARIMA模里、确定性分析法，名元时间序列建模
一、 实验目的和要求
1•熟悉非平稳序列的确定性分析法：趋势分析、季节效应分析、综合分析 2•熟悉差分平稳序列的建模步骤。
3•掌握单位根检验、枷整检验、动态回旧模里的理立。
二、 实验原理
1. 序列的各种变化部脳于四大因素的绘合时长期趣势(Trend ) ,1 环波动(Circle )，季节性变化(Season )，机波动(Immediate ) •常個设它们有如下 的相互模型:
加法模型XT+q+y + 厶 乘法模型X,=7；CfSfZ,
混合模型模型结构不唯一
2•非平稳序列如果能通il适当阶数的差分后实观平稳，就可以对差分后序列 S II ARMA模里合了，所以ARIMA模里是差分运算与ARMA模里的组合
=e(B)s,
3•单位根检验：
(1 ) DF 检验；(2 ) ADF 检验；(3)PP 检验; 4•动态回1月模塑ARIMAX
如果两个非平稳序列之同具有协整关系,则先建立它们的回归模型，再对平梯的 残
•可修編.
差序列建立ARMA模型。
9
±0(B)B， 台® (B)
y.=
£i 0(B) = 6(B)®
三、实验容
1、 P202贡：第7题(X11因索分解注) 2、 P155页：第3题(乘枳季节模里) 3、 P240页：第4题
出口为心，进口为儿，回答以下问趣
(1 ) B IB xf,儿的时序图，用单位根检验序列它们的平WK; (2) 对In旺,In”分别«1合模里(提示：建立ARIMA模里)；
(3) 考察in”,In乞的协整关系,建立lux关于In旺的协整模塑,同时建立淚差 修正模型。
•可修編.
EL实验过程
（-）P202页：第7 g （X11因素分解法） 1•绘胃序列时序图。（程序见附录）
由上图可得季节序列的掾帽酿序列水平的变化而变化，所以季节效应与茜势效应不狼立，采 用乘选模里：xt = rf X5t X it
2•进Ax-11季节调整模型
经过三个阶段共十步的亜夏迭代后，得到如下的抓合效果图:
•可修編.
显然，孩地区妍牛的月度产妍量序列具有显普的季节变动特征。
（二）P155页：第3题（乘积季节模型）
1.绘制序列时序图。绘制时序图,如图1所示（程序见附录1）。
图1美国月度事枚死亡人数序列时序图
时序图显示该序列具有以年为周期的李节效应。
2 •差分平稳化：对原序列作1阶12步差分，帝望提取原序列季节效应，差 分后序列时序图如图2所示。
•可修編.
difl 12
2 WO
1973
1973
W73
1974
!974
1974
W75
1975 •WS 1976
ties
1976
1976
197?
1977
1977
1978
-1WO “ _____ ______ ______ ___ _____ I ................................... I .................. .... ....... ........... .... ........ ... ........................ ............ ................ ..... .... ...............
|«I Vl|l|aw ril V||VVVIT||I »II TT1|«1¥I »II IVITVII Ijl VTVTTVJV II ITVIJI IVTITV |fl VIVTVjT VVIVTIJI IV Till |ri VTVTVjll »1I I||
1973
1973
!9?9
图2美国月度事枚死亡人数1阶12形差分后序列时序图
时序图显示差分后序列类似平隐。
3 •模塑定阶：考察差分后序列自相关图，如图3,进一步确定平IBUHIBi, 并估计《1合模型的阶数。
图3美国月度事枚死亡人ft 12^差分后序列自Hl关图 自相关图显示延迟12阶自相
关系数显普大于2倍标准差围，这说明差分后序列 中的组含着非常显普的季节效应。延迟1阶的自相关系数也大于2倍的标准差, 这说明差分后序列还具有姬期相关性。
观察偏自相关图，如图4,得到的结论和上面的结论一致。
•可修編.
美国月度事枚死亡人数1阶12步差分后序列IS0HI关图
图5序列白課声检验
图5显示，原序列延迟各阶LB貌廿量的P值小于显普性水平0.05,所以拒 绝原假设，序列不通过白噪声检验。
根弼差分后序列的自相关图和偏自相关图的111贯，拥合乘枳季节模型 ARIMA （p,d、q）x（P、D, Q）l2 o
自相关图显示，12阶以的自相关系数1阶截尾，偏自相关图显示，12阶以的偏 自相关系数1阶截尾，所以尝试使用ARMA（1, 0）模塑提取差分后序列的短期自 相关信息。
再考虑季节自相关特征，这时考察延迟12阶、24阶等以周期长度为单位的 自相关系数和偏自相关系数的特征。自相关图显示延迟12阶自相关系数显普非 零，而偏自相关图显示延迟12阶偏自相关系数显著非零，这时用1112步为周期 的ARMAdJh模塑提取差分后序列的季节自相关信息。
•可修編.
Auto reg ressive Factors
Factor 1： 1 - 0.48078
Moving Average Factors
Factor 1:
郴⑴
Factor 2： 1 - 0.53808 B删(12)
1 - 0.87376 B
图6拥合模里
综合前面的差分信息，我们要#1合的乘枳季节模型为AN/MA（（1,1J）|2。
使用条件最小二乘估it方法,确定该模型的口径为:
叫严上坐上如巴
（1一酬）
(l-0.87376B)(l-0.53808B12)「
(1-0.490783) 匂
•可修編.
1丄O）X
5 •模里检验：对序列合AMMA（l」,O）x（l丄1九模型，模型艮模型参数的显普 性检豔如图7、8所示。
图7模型参数的显普性
由图7知，》1合效果显示模型参数显普。
Autocorrelation Check of Res i dua I s
To Lag
6 12 18 24
Chi-
Square
Pr >
DF 3
3 15 21
ChiSq 0.4435 0.4146 0.3206 0.3789
Autocorre lat ions ---------------------------------------------
2.G8
3.25 16.97 22.34 -0.045 -0.211 -0.006 -0.030
0.055 -0.037 0.160 -0.013
0.137
0.079 -0.045 -0.071 -0.123 -0.043 -0.086 0.027 0.046 0.186 0.152 0.190 0.020 0.051 -0.182 -0.067
图8歿差白噪声检验
对抑合模型进行白嗥声检脸,结果显示p值31大于显著性水平0・05•接受原假设, 歿差序列通11白噪声检騎，模里显普，说明模型m合良好，对序列相关信息提取 充分。
将序列#1合值和序列观察值联合作图,如图9所示。
•可修編.
图9美国月度事放死亡人数扭合效果图 说明:图中,点为序
列观察值;曲裁为序列81合值。
从图9可以直观地看出该乘枳季节模型对原序列的#1合效果良好。
（三）P240页：第4题
1.画岀儿的时序图，用单位根检验序列的平稳性; 输出时序图如图1所示（程序见附录2）。
图1我国出口总额Xt、进口总额yt时序图
图1巾，黑色为岀口总额Xt序列时序图，红色为进口总额yt序列时序图。 从图1巾可以看出出口总额xt序列.进口总额yt序列均显著非平稳，迪个直观 判Bi还可以通il单位根检验验证。同时时序图显示这两个序列具有某种同变关 系O
对我国岀口总制序列xt进行ADF检验，单位根检验结果如图2所示。
•可修編.
图2出口总额xt白噪声、单位根检验
检騎结果显示，无论考虑何种类塑的模型，检騎统计量的P值均显普大干 0.05的显普性水平，所以可以认为中国我国岀口总额序列xt显普非平稳，目这 兀种处理均不能实现残差序列平稳。
对我国进口总额序列VtjJ行ADF检验,单位根检验结果如图3所示。
•可修編.
The AREMA Procedure
Autocorrelation Check for White Noise
Chi- Pr > Square DF
ChiSq ................................................. Autocorrelat i ons1
124.98 6
<0001 0.856 0.708 0.575 0.455 0.347 fl.273 12
135.76
12
<.0001
0.227
0.191
0.15G
0.123
0.114
fl.090
Augrne nted Dickey-Ful I er Unit Root Tests Ti»pe Lae® Rho Pr < Rho
Tau Pr < Tau
F
Pr > F
Zero Meo.n
0
9.2683 0.9999 13.62 0.9999
1 8.4853 0.3393 1.68 0.9762
Single Mean
0 8.36T5 0.3333 11.32 0.3333 33.07 0.0UI0
1 7.7085 0.8989 1.53 0.9892 1.86 0.6021 Trend
0 7.8871 0.9999 7.93 0.9999 74.SA 0.0010
1
4.9889
0.9999
1.00
0.9999
2.63 0.6560
图3进口总额yt白噪声、单位根检验
同出口序列xt的检验结果一样，在显普11水平取为0.05 N ,可以认为我国口序列yt非平稳，USA种处理均不能实现歿差序列平稳。
显然，SMf序列的ADP检騎结果与根据时序图得到的直观判斷完全一致
2. 对分别抑合模型（提示:建立ARIMA模型）；
对我国出口对数序列Inxt和进口对数序列Inyt绘制时序图，如图4所示。
•可修編.
进
图4我国岀口忠额Xt、迪口忠額yt皿对数时序图
图4巾，黒色线代表我国出口对数序列Inxt,红色线代表我国进口对数序列Inyto 时序图显示这两个对数序列有显普的上升范衿，为典塑的非平稳序列。同时时序 图显示泄两个序列具有杲种同变关系。
因为序列呈现岀近做线性范矜，所以选择1阶差分。1阶差分后岀口对数序 列Inxt时序图如图5所示。
图5对数序列Lnx差分时序图
时序图显示，Inxt差分后序列在均值附近比较稳定地彼动。为了进一步确定平稳 性，考察差分后序列的自相关图，如图6所示。
•可修編.
图6对数序列Lnxt差分后自柚关图
自相关图显示序列有很強的姬期相关性，所以可以初步认为Inxtl阶差分后 序列平隐。
对平稳的1阶差分序列进行白蝶声检騎，白蝶声检騎结果如图7所示。
图7lnxt —阶差分后序列白噪声枪验
•可修編.
在检騎的显普性水平取为0.05 «条件下，由于延迟6阶、12阶的Pffifi小 T0.05,所以Inxt差分后的序列不能视为白蝶声序列，即差分后序列“含着不 容忽視的相关信息可以提取。
对平稳非白蝶声差分序列抓合ARMA模里,1阶差分后序列的自相关图（见图 6）已经显示该序列有不截尾的性质。再考察其偏自相关系数的性质，如图8所 爪O
图8对数序列Lnxt差分后倫自柚关图
偏自相关图显示出1阶截尾也所以考虑用AR（1）模型抓合Inxtl阶差分后序 列。考虑到前面已经进行的1阶差分运算，实际上是用AQM/UhO）模里抵合原 序列。对序列#1合A/?/MAU,O）模里，模里参数及模里的显普性检騎如图9、10 所示。
•可修編.
a
图9模塑参数显普性检验
由图9知，系数显普性检验显示两参数均显普。 对残差序列进行白蝶声检騎，检豔结果如图10所示。
图10残差白礫声检验
显然，#1合检验统廿量的P值都显普大于显普性检验水平0.05,可以认为残差序即为白噪声序列，模型显普，这说明A/?/MAU,O）模里对Inxt序列建模成功。
•可修編.
列
图门模型
在条件最小二乘估廿原理下,81合结果为:
刊\"。皿9+匸需才2
将对数序列#1合值Inxt和对数序列观寮值InxtR合作图,如图12所示。
图12对数序列Lnxtffl合效果图 说明：图中，星号为序
列观察值；曲线为拥合值。
从图可以直观地看岀该AKIMBO）模型对原序列的Jfl合效果良好。
因为对数序列Inyt呈现出近111线性赶势，所以选择1阶差分。1阶差分后进
•可修編.
口对数序列Inyt时序图如图13所示。
difIny
-fl
1950 I960 197Q 1980 year
1W) 2000 20^0
@13对数序列Lny差分时序图
时序图显示,Inyt差分后序列在均值附近比较稳定地iOJo为了进一步确罡平稳 性，考察差分后序列的自相关图，如图14所示。
图14对数序列Lnyt差分后自相关图
自相关图显示序列有很强的短期相关性,所以可以侨步认为Inytl阶差分后 序列平隐。
对平隐的1阶差分序列1JH亍白噪声检验，白蝶声检醴结果血图15所示。
•可修編.
图15 Inyt -阶差分后序列白噪声检验
在检豔的显著性水平取为0.05的条件下,由于延迟6阶的P值小于0.05, 所以Inyt差分后的序列不能视为白蝶声序列，即差分后序列还组含着不容忽视的 相关信息可以提取。
对平隐非白噪声差分序列合ARMA模里,1阶差分后序列的自相关图（见图 14）已经显示该序列有1阶截尾的性质。再考察其偏自相关系数的性质，如图 16所示。
图16对数序列Lnyt差分后倫自柚关图
倫自相关图显示该序列1阶截尾的性飢所以考虑用AR⑴模型81合Inyt1 阶差分后序列。考虑到前面已经进行的1阶差分运算，实际上是用ARIMW丄0) 模型抓合原序列。对序列拥合AMMA1,1,O)模塑，模里参数及模里的显普性检验 如图17、18所示。
•可修編.
a
图17模里参数显着性检验
由图170,系数显普性检验显示两参数均显普。
•可修編.
对残差序列进行白繰声检验,检脸结果如图18所示。
图18残差白課声检验
显然，#1合检验统廿量的P值都显普大于显著II检验水平0.05,可以认为残差序 列即为白噪声序列，模型显普。这说明AQAM1,1,O）模里对该序列建模成功。
图19模里
•可修編.
在条件最小二乘估廿原J!下，jfl合结果为:
5 272+匚窗护
将对数序列抑合值Inyt和对数序列观寮值lnyt»合作图,如图20所示。
@20对数序« Lnyt ffl合效果图 说明：图中，星号为序
列现察值；曲找为ffl&ISo
从图20可以直观地看岀该丄0)模里对原序列的#1合效果良好。
3. 考察lnJz,lnxf的协整关系,建立hi几关于in旺的协整模塑,同时建立淚差 修正模里。
对我国出口对数序列Inxt fll进口对数序lnyt绘制时序图，如图4所示。可以发 规时序图显示迪两个序列具有某种同变关系，可以考虑建立ARIMAX模塑。对 Inxt、lnyt、Inxtl阶差分、(Vlny,)序列分别进行单位根检验(ADF )o输岀结果 则图21 24所示。
图21对数序列InXtl阶单位根检验
检验结果显示，无论考虑何种类型的模型，检验貌it量的P值均显普大于 0.05的显普性水平，所以可以认为中国我国岀口总额对数序列Inxt显普非平隐,
•可修編.
a
A2A种处理均不能实现歿差序列平隐。
图22对敛序列Inytl阶单位根枪验
同岀口对数序列Inxt的检騎结果一样，在显普性水平取为0.05 IH ,可以认 为我国进口对数序列Inyt非平稳,且这穴种处理均不能实现残差序列平稳。
图23 {Vlnxf}1附单位根枪验
检验第果显示，无论考虑何种类型的模型,检脸绒计量的P值均显著小于0.05 的显普性水平，拒绝原假设，所以可以从为中国我国出口总額{叫}对序列显 著平稳。
•可修編.
a
图24 (VlnyJI阶单位根枪验
同出口 (Vlnxf)序列的检验结果一样，在显普牲水平取为0.05BL可以从为 我国进口 {Vlny,}序列平稳，USA种处理均能实现残差序列平稳。
•可修編.
利用最小二乘估廿,回1月模型输出结果如图25所示。
图25回IH模里结果
构造出的回归模型如下:
In yz = 0.99265 In xt + ®
图26对数序列Lnyt与对数序列Inxt之间的相关图
相关图显示对数序列Inyt在延迟阶数为零时与对数序列Inxt相关关系最大。 因此可以将对数序列Inyt与对数序列Inxt同期建模。
•可修編.
a
图27残差单位根检验
歿差序列平隐，说明对数序列Inyt与对数序列Inxt之间具有协整的关系,我们可 以大IB的在这两个对数序列之同建立动盗回归模型而不必担心虚假回IJ3问题。
图28模型参数显著性检勁，无常数顶
考察歿差序列白噪声检验第果，如图29所示。
•可修編.
a
图29歿差序列白噪声枪验
输岀结果显示，延迟各阶LB筑计量的P值都大于显普性水平0.05,可以认 为残差序列为白噪声序列,对序列相关信息提取充分。
图30模型
•可修編.
a
根据输出的模型#1合结果可知，最后的#1合模型口径为:
1
In yt =0・991791n 兀 +
1-0.69934B®
•可修編.
将对数序列#1合值Inyt和对数序列观寮值InytJK合作图，如图31所示。
year
图31对数序« Lnyt ffl合效果图
说明：图中，星号为序列现察值；曲线为抓合H。
从图可以直观地看岀该ARIMAX模里对原序列的#1合效果良好。 将序列#1合值yt利序列观察值ytIR合作图,如图32所示。
•可修編.
说明：图中，星号为序列现察值；曲线为拥合值。
从图可以直观地看岀该ARIMAX模型对原序列的#1合效果良好。 沟造误差修正模型：
为了研究我国的进出口总额的短期波动特征，我们利用差分序列
M锻差序列
｛vhl
”｝和 以及前
｛ECM｝
^，构造ECM模里：
ECM
一
=In0.99265
输出给果如图33所。
Paranat?r Est imates
Variable diCInx I a^ecn
DF 1 1
Parameter Esl i rralc 1.0219? -0.31055
Standard
Error 0.08092 0.10004
l Value 12.63 ■3.10
P
「 > Itl
<.0001 0.0030
图33误差修正模型显著11检射
由图33输出结果结果知ECM模型为:
V In ” = 1.02197 Vln 兀一 0.31055 ECM一 + st
附录
椁序1:
data example4_7; input x;
tdntnxCquarter1 ,'1 jan1978'd,_n_-1); format tyyq4.; cards;
589 561 640 656 727 697 640 599 568 577 553 582
•可修編.
600 566 653 673 742 716 660 617 583 587 565 598 628 618 688 705 770 736 678 639 604 611 594 634 658 622 709 722 782 756 702 653 615 621 602 635 677 635 736 755 811 798 735 697 661 667 645 688 713 667 762 784 837 817 767 722 681 687 660 698 717 696 775 796 858 826 783 740 701 706 677 711 734 690 785 805 871 845 801 764 725 723 690 734 750 707 807 824 886 859 819 783 740 747 711 751
procxl 1 data=example4_7; quarterly date=t; var x;
output out=out b1=xd10=season d11=adjusted d12=trend d13=irr; data out; set out;
estimate=trend *season/100; procoplotdata=out;
plot season*t=2 adjusted*t=2 trend*t=2 irr*t=2; plot x*t=1 estimate *t=2/overlay; symbol1c=black i=join v=star; symbol2c=red i=join v=none w=2; run;
程序2：
data ti5_3; input x;
dif1J2=dif12(dif(x));
time=intnx(,month,,,1jan1973,d1_n_-1); format time year4.; 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
；run;
procgplot;
plot x*time=1 dif1_12*time=2; symbollc=coral v=circle i=join; symbol2c=blue v=star i=join;run; procarima; identifyvar=x(1,12);
estimatep=1q=(1)(12); forecastlead=Oid=time out=out;run; procoplotdata=out; plot x*time=1 forecast*time=2 /overlay;
•可修編.
symbol1c=black i=nonev=dot h=0.2; symbol2c=red i=join v=none; run;
程序3: data ti6_4; input year xt yt; lnxt=log(xt); lnyt=log(yt); diflnx=dif(lnxt); diflny=dif(lnyt); cards;
195020 21.3 195124.2 35.3 195227.1 37.5
195334.8 46.1 195440 44.7 195548.7 195655.7 195754.5
195978.1 196063.3 196147.7 196247.1
196455.4 196563.1
196758.8 196857.6 196959.8 197056.8 197168.5 61.1 53 50 195867 61.7
71.2 65.1 43 33.8 196350 35.7
42.1 55.3 196666 61.1
53.4 50.9 47.2 56.1 52.4
.
•可修編
197282.9 64 1973116.9 103.6 1974.4 152.8
1975143 147.4 1976134.8 129.3 1977.7 132.8
1978167.6 187.4 1979211.7 242.9 1980271.2 298.8 1981367.6 367.7 1982413.8 357.5 1983438.3 421.8 1984580.5 620.5 1985808.9
1257.8 19861082.1 1498.3 198714701614.2 19881766.7
2055.1 198919562199.9
19902985.8 2574.3 19913827.1 3398.7
19924676.3 4443.3 19935284.8 5986.2 199410421.8 9960.1
199512451.8 11048.1
199612576.4
11557.4 199715160.7 11806.5 199815223.6 11626.1
199916159.8 13736.5 200020634.4 18638.8 200122024.4 20159.2
200226947.9 24430.3 200336287.9 34195.6 200449103.3 46435.8
200562648.1 54273.7
200677594.6
63376.9
200793455.6 73284.6
2008100394.9 79526.5 ;run; procgplotdata=ti6_4;
plot xt*year=1 yt*year=2/overlay; symbol1c=black i=join v=none;
symbol2c=red i=join v=nonew=2l=2;run;procarimadata=ti6_4;
identifyvar=xt stationarity=(adf=1);
.
•可修編
identifyvar=yt stationarity=(adf=1 );run; procoplotdata=ti6_4;
plot lnxt*year=1 lnyt*year=2/overlay; plot diflnx*year=1 diflny*year=2; symbol1c=black i=join v=circle;
symbol2c=red i=join v=star;run; procarima; identifyvar=lnxt(1); estimatep=1;
forecastlead=Oid=year out=out1; identifyvar=lnyt(1); estimatep=1;
forecastlead=Oid=year out=out2;run; procgplotdata=out1 ;
plot lnxt*year=1 forecast*year=2 /overlay; symbol1c=black i=nonev=star; symbol2c=red i=join v=none;run; procoplotdata=out2;
plot lnyt*year=1 forecast*year=2 /overlay; symbol1c=black i=nonev=star; symbol2c=red i=join v=none;run; procarimadata=ti6_4;
identifyvar=lnxt stationarity=(adf=1); identifyvar=lnyt stationarity=(adf=1); identifyvar=lnxt(1) stationarity=(adf=1); identifyvar=lnyt(1) stationarity=(adf=1 );run; procreodata=ti6_4; model lnyt=lnxt/noint;
outputout=out residual=residual;run; procarimadata=out;
identifyvar=residual stationarity=(adf);run; procarimadata=ti6_4;
identifyvar=lnyt crosscorr=(lnxt); estimatep=1input=lnxt noint;
forecastlead=Oid=yearout=result;run; data result; set result; yt=exp(lnyt);
estimate=exp(forecast);run; procoplotdata=result;
plot lnyt*year=1 forecast*year=2 /overlay; plot yt*year=1 estimate *year=2 /overlay; symboll c=black i=non ev=star; symbol2c=red i=join v=none;run; data ti6_4j
.•可修編
set ti6_4;
ecm=lnyt- 0.99265 *lnxt; lag_ecm=lag(ecm);
diflnx=dif(lnxt); diflny=dif(lnyt);run; procreodata=ti6_4j model diflny=diflnx lag_ecm /noint; run;
•可修編.
因篇幅问题不能全部显示，请点此查看更多更全内容