南邮时间序列实验报告

南京邮电大学

实验报告

实验名称：AR(p)模型的建立

ARMA模型

课程名称：应用时间序列

班级学号：B11080404

姓名：陈海霞

开课时间：2013 /2014 学年第二学期

实验 一 AR(p)模型的建立

一、 实验题目

设{}t ε是均值为0，方差为4的白噪声序列，(4)AR 模型的自回归系数为：

12340.9, 1.4,0.7,0.6a a a a =-=-=-=-，

(1) 在计算机上模拟产生一个符合此模型的长为505的序列片断 (2) 用以上的前500个数据对时间序列进行建模

(3) 用递推预测法预测后5个数据，与真实数据作比较，检验预测效果。

二、 实验原理及结果分析

（1）AR （4）模型为：12340.9 1.40.70.6t t t t t t X X X X X ε----=----+，{}t ε是WN （0,4）。 首先产生706个白噪声序列{}t ε，在设定序列t {X }前4个值为0的情况下，产生后面706个t X 。取201X 到706X 之间505个数据，作为符合此模型的片段序列。 长为505的序列片断为

Columns 1 through 7

-1.6486 -1.2827 0.4863 3.6305 -2.1616 2.3075 -3.2600 Columns 8 through 14

-4.0617 6.6829 1.2724 -3.2541 -2.2762 6.0537 -0.8093 Columns 15 through 21

-2.7988 -0.4980 4.0827 -2.8195 0.1556 2.6400 -2.0984 Columns 22 through 28

3.0469 -3.0930 5.0175 -1.8583 -2.6975 3.9419 -2.4583 Columns 29 through 35

-0.9778 4.6628 -4.1501 3.4180 0.7833 -6.6836 4.8753 Columns 36 through 42

3.1626 -5.6393 2.7466 -0.0131 -0.5746 5.0387 -5.4597 Columns 43 through 49

-1.1448 9.3258 -7.2781 -1.8489 9.0824 -8.0117 -2.1088 Columns 50 through 56

4.9452 -3.2509 2.1783 -3.7674 -3.9564 6.2281 0.7752 Columns 57 through 63

-2.0820 -1.2922 -1.2861 3.7972 4.8990 -7.6666 -1.5087 Columns 64 through 70

6.9074 -4.6859 -4.4424 2.6415 5.3233 -3.3085 -4.9210 Columns 71 through 77

3.2859 3.3753 -0.4918 -6.0565 2.0189 6.5229 -6.2671 Columns 78 through 84

-5.2132 8.1316 0.7140 0.4299 -4.0847 -2.3634 6.0798 Columns 85 through 91

2.7587 -6.7040 2.6006 2.3261 -0.6494 -0.7489 -1.2697 Columns 92 through 98

0.3373 1.1502 -3.9484 -0.3369 5.3149 -2.9744 -4.1054 Columns 99 through 105

6.9682 0.0271 -5.9869 4.7977 -2.3928 3.5680 0.3574 Columns 106 through 112

-5.0968 4.0452 -0.2734 0.0166 1.7894 -5.3719 2.9166 Columns 113 through 119

1.2599 -3.1610 -0.7623 1.2415 0.9579 1.9063 -3.8318 Columns 120 through 126

-2.2437 5.4332 -2.5049 -2.7726 6.5762 -3.3202 -4.3710 Columns 127 through 133

5.4127 -1.7038 -2.3969 2.9173 -2.2327 0.8692 4.7882 Columns 134 through 140

-6.3225 -0.9963 7.5435 -2.2728 -2.8918 1.1488 -2.1442 Columns 141 through 147

5.0079 -1.1340 -4.8897 2.8174 2.9483 -3.0528 -1.5609 Columns 148 through 154

-0.8641 1.2427 6.4252 -1.9073 -9.6514 5.8646 5.8059 Columns 155 through 161

-6.4211 2.8790 -0.1245 -4.1447 3.2837 1.3619 -2.7480 Columns 162 through 168

-1.9474 3.9604 0.3485 -0.3001 -2.3046 0.8466 1.1656 Columns 169 through 175

-0.4441 1.8892 0.3142 -3.1781 6.2901 -2.7625 -1.8085 Columns 176 through 182

3.4508 -3.0910

4.4346 -0.4656 -6.2508 4.6428 2.8445 Columns 183 through 189

-4.4477 0.6162 0.5373 -1.8424 3.0585 -3.9518 -0.0967 Columns 190 through 196

4.7312 -2.5664 -0.1160 -1.8677 -1.5911 9.5630 0.2328 Columns 197 through 203

-9.6110 1.4483 7.3692 -1.7266 -8.0665 3.7535 4.4337 Columns 204 through 210

-1.7092 -2.3077 -7.7363 6.4960 4.7370 -8.1310 -2.8658 Columns 211 through 217

6.2017 -0.5482 -0.1900 -3.0736 0.8053 1.3056 -3.7892 Columns 218 through 224

0.8969 4.2618 -1.0159 -4.6282 5.0698 -0.1972 -2.8631 Columns 225 through 231

1.4848 -1.3444 1.7330

2.0872 -4.6884 -4.3782 12.1706 Columns 232 through 238

-4.6890 -3.2277 5.5255 -9.7148 6.0399 6.6920 -9.4818 Columns 239 through 245

1.8886 4.3605 -3.4107

2.0981 1.9325 -

3.9338 1.0377 Columns 246 through 252

1.9502 -0.4143 -0.8270 -

2.2374 0.2333 4.3005 -0.2618 Columns 253 through 259

-2.0929 0.8246 -1.6636 3.5307 -2.1545 -2.6078 1.2247 Columns 260 through 266

4.4987 -1.7981 -2.1969 0.3333 1.8965 -1.6881 -2.6223 Columns 267 through 273

1.1559 6.1562 -3.6269 -5.2409 3.7809 0.6854

2.4372 Columns 274 through 280

-0.7770 -4.6102 3.9574 4.4686 -3.6919 -0.1215 -2.2027 Columns 281 through 287

0.7822 6.9943 -3.6771 -3.6098 3.2648 0.9603 -1.3694 Columns 288 through 294

-0.3355 2.2325 -0.5881 -1.0960 1.2105 -3.7500 6.7471 Columns 295 through 301

-2.0846 -1.5453 1.7557 -0.7898 -0.3974 4.2752 -3.6040 Columns 302 through 308

-1.9831 5.1484 0.2379 -5.9040 -0.0958 6.4623 -1.1515 Columns 309 through 315

-7.8338 6.0314 1.2688 -2.4958 3.1931 -2.0916 -2.9175 Columns 316 through 322

4.7176 1.1561 -1.9735 -0.2712 2.3305 0.7569 -1.5717 Columns 323 through 329

-2.5402 3.1030 2.8523 -2.1630 -2.7434 1.4132 4.5201 Columns 330 through 336

-0.3166 -5.3042 -0.1754 5.1126 -0.0785 -2.5813 -7.4676 Columns 337 through 343

10.1845 -0.3370 -8.9094 2.7457 2.5486 -2.3195 2.6970 Columns 344 through 350

0.1606 -5.5921 2.9602 1.0063 0.0901 -2.5625 -0.5199 Columns 351 through 357

1.6619 0.1863 -0.0225 0.2041 -3.9329 1.0144

2.9189 Columns 358 through 364

-2.9282 1.2670 -0.8981 -2.3965 -0.4082 2.7991 1.0941 Columns 365 through 371

-1.4877 -0.5528 -1.1695 0.6171 4.9160 -9.2752 4.0470 Columns 372 through 378

5.7613 -9.2379 3.9880 1.2773 -1.4542 5.6733 -3.0832 Columns 379 through 385

-5.2407 4.2054 0.6001 -2.9886 0.7615 -3.7765 4.8305 Columns 386 through 392

3.3343 -8.0033 -0.2696 6.2272 -2.8318 -1.8286 2.5073 Columns 393 through 399

-0.2841 -2.8274 4.1081 -2.2792 -3.9379 0.4100 1.5676 Columns 400 through 406

2.1587 -

3.4707 -

4.8650 8.7234 3.2379 -9.9828 1.9502 Columns 407 through 413

7.3327 -5.2973 -2.2734 4.3372 -1.6377 -0.6445 0.4255 Columns 414 through 420

-0.7174 3.5823 -1.2427 -1.2148 3.0170 -1.9085 -1.8497 Columns 421 through 427

1.5624 1.5779 0.6372 -1.9837 -0.0816 -0.0610 -0.5280 Columns 428 through 434

-2.6395 6.4001 -3.6248 -2.8943 7.1429 -3.5746 -5.6710 Columns 435 through 441

3.3058 5.6545 -5.4648 -1.3714

4.7310 0.6276 -2.3320 Columns 442 through 448

-0.3020 0.3274 1.6821 -1.2582 -2.3358 2.1568 0.3293 Columns 449 through 455

-0.7688 1.1805 -1.5014 1.7361 -0.3484 0.6664 -2.4582 Columns 456 through 462

0.6786 -0.0852 -1.1447 0.2410 1.9136 -1.8712 -2.0340 Columns 463 through 469

3.1176 0.0792 -2.6877 -1.0560

4.1052 -0.1145 -2.0853 Columns 470 through 476

2.3643 -4.0368 2.3417 0.7316 -

3.1192 2.7407 -2.7758 Columns 477 through 483

-2.1756 6.2223 -2.1149 -1.8809 -2.0565 3.2973 3.4009 Columns 484 through 490

-2.3080 -3.8672 2.6712 2.5792 -3.8149 1.1462 -1.9989 Columns 491 through 497

-0.3990 4.7531 -1.8109 -2.1938 1.8063 3.5909 -0.3040 Columns 498 through 504

-6.2670 5.0454 3.3294 -1.7735 -1.6727 -1.3086 3.5864 Column 505

4.8878

(2)取505个数据建立模型时，观察自相关系数和偏相关系数

自相关系数

偏相关系数

根据图可以看到偏相关系数是四阶截尾的。所以，可以建立AR(4)模型。 (3) AR(4)模型参数的估计 通过Yule-Walker 方程估计参数

n n n a γ1

-Γ=

结果为

所以模型为：

4-t 321X 5432.06633.03775.18945.0----=---t t t t X X X X

（4）白噪声2χ检验

对独立同分步的白噪声t ε，用j ρ

?表示基于观测值数据N εεε...,21的样本自相关系数，对于很大的N ，有

)?...??()(?222212m N m ρρρχ

+++= 原假设：{}t H ε:0是独立白噪声；对立假设{}t H ε:1是相关序列

给定检验水平05.0=α，查2

χ分布表得到临界值αλ，满足

αλχα=>))((2m P

当计算结果αλρρρ

>+++)?...??(2

2221m N 时，应该拒绝原假设，t ε不是白噪声， 当αλρρρ

<+++)?...??(2

2

22

1m N 时，应该接受原假设。t ε是白噪声。 实验结果显示，是接受原假设。

（5）递推预测

通过建立的模型，用递推预测法预测后5个数据yy ，真实数据为YY

定义误差比为

YY

YY

yy -=

e 画出误差图

看出误差是在可控制范围内

三、 实验代码

a1=-0.9;a2=-1.4;a3=-0.7;a4=-0.6;x(1)=0;x(2)=0;x(3)=0;x(4)=0;r=zeros(1,30); e=randn(1,706);

e=2*(e-mean(e))/std(e); for i=5:710

x(i)=a1*x(i-1)+a2*x(i-2)+a3*x(i-3)+a4*x(i-4)+e(i-4); end

for i=1:505

Y(i)=x(i+200); end Y;

y=Y(1:505);

figure(1);autocorr(y); figure(2);parcorr(y); ybar=mean(y); y=y-ybar;

gamma0=cov(y); for j=1:4

gamma(j)=sum(y(1:500-j).*y(1+j:500))/500; end

for i=1:4 for j=1:4 if(i==j)

Ga(i,j)=gamma0; else if(i>j)

Ga(i,j)=gamma(i-j);

else Ga(i,j)=gamma(j-i); end end end end

a=(inv(Ga)*gamma(1:4)')'; for i=501:505

y(i)=y(i-1)*a(1)+y(i-2)*a(2)+y(i-3)*a(3)+y(i-4)*a(4)+ybar; end

yy=y(501:505) YY=Y(501:505)

实验 二 ARMA 模型及预测

一、实验题目

为课本附录8的化工浓度数据建立ARMA 模型，并作一步和两步预测

二、实验原理及结果分析

（1）判断模型是否平稳

对原始数据处理，画出时间序列图

根据发现序列不平稳，所以做一阶差分 1--=t t t X X Y 得出结果是平稳的

（2）确定平稳模型

通过看自相关系数图和偏系数图来确定模型

偏相关系数

自相关系数

由图可知，y 的自相关系数一阶截尾，很有可能为MA （1）模型。 （3）采用AIC 函数对MA 模型定阶，并对参数进行估计

根据图中观察得到，MA 模型不会超过MA(5),假定m=5，计算5个AIC 值，比较看哪个最小，进而确定模型。根据matlab 编程得到的AIC 值为

AIC 值最小的为第一个，即确定模型为MA （1）模型，参数估计得

即15276

.0--=t t t y εε （4）白噪声2χ检验

对独立同分步的白噪声t ε，用j ρ?表示基于观测值数据N εεε...,21的样本自相关系数，对于很大的N ，有

)?...??()(?222212m N m ρρρχ

+++= 原假设：{}t H ε:0是独立白噪声；对立假设{}t H ε:1是相关序列 给定检验水平05.0=α，查2χ分布表得到临界值αλ，满足

αλχα=>))((2m P

当计算结果αλρρρ

>+++)?...??(22221m N 时，应该拒绝原假设，t ε不是白噪声， 当αλρρρ

<+++)?...??(2

2

22

1m N 时，应该接受原假设。t ε是白噪声。 实验结果显示，是接受原假设。

（5）模型预测

利用信息递推估计方法，可以进行一步预测，预测出

197y 的值，进而计算出

19817.4001 X =

因为这是MA （1）模型，新息的预测只有一步，所以198y 无法预测，即199X 也无法

计算，所以

1990X =

三、实验代码

load 'data.txt'

plot(data)

title('原始数据序列图'); for i=1:196

y(i)=data(i+1)-data(i); end plot(y)

title('一阶差分后的时间序列图'); figure(3);

autocorr(y); %作自相关系数图 figure(4);

parcorr(y); %作偏相关系数图 mu=mean(y);

gamma0=sum((y-mu).^2)/196; for k=1:30

gamma(k)=sum((y(1:196-k)-mu).*(y(1+k:196)-mu))/196; end

sigma_2=zeros(5,1);

AIC=zeros(5,1);

for q=1:5

C=zeros(q,1);

C(1,1)=1;

gamma_q=zeros(q,1);

for i=1:q

gamma_q(i,1)=gamma(i);

end

for k=2:20

Gamma_k=zeros(k,k);

for i=1:k

for j=1:k

if(i==j)

Gamma_k(i,j)=gamma0;

else if (i>j)

Gamma_k(i,j)=gamma(i-j);

else Gamma_k(i,j)=gamma(j-i);

end

end

end

end

omega=zeros(q,k);

for j=1:q

for t=1:k

omega(j,t)=gamma(j+t-1);

end

end

pi{k-1}=omega*inv(Gamma_k)*omega';

end

i=1;

while abs(det(pi{i})-det(pi{i+1}))>0.00001

i=i+1;

end

sigma_2(q,1)=gamma0-C'*pi{i}*C; % sigam平方的估计AIC(q,1)=log(sigma_2(q,1))+2*q/196; %定阶公式

end

AIC

rou=gamma(1)/gamma0;

b=(1-sqrt(1-4*rou^2))/(2*rou)

%检验白噪声

e0=0;e(1)=-b*e0+y(1);

for i=2:196

e(i)=-b*e(i-1)+y(i);

end

for k=1:186

ep(k)=e(k+10);

end

ebar=mean(ep);

for j=1:20

rou(j)=sum((ep(1:186-j)-ebar).*(ep(1+j:186)-ebar))/sum((ep-ebar).^2);

chi2(j)=186*sum(rou.^2);

fws(j)=chi2inv(0.95,j);

end

figure(5);

plot(chi2);

hold on

plot(fws);

for k=1:196

gamma(k)=sum((y(1:196-k)-mu).*(y(1+k:196)-mu))/196;

end

theta=[];v0=gamma0;theta(1,1)=gamma(1)/v0;

v(1)=gamma0-theta(1,1).^2*v0;

for n=2:196

theta(n,n)=gamma(n)/v0;

for m=1:n-1

k=n-m;

theta(n,m)=(gamma(m)-theta(k,k)*theta(n,n)*v0-(theta(k,1:k-1).*theta(n,m+1:n-1))*flipud(v(1:k-1)'))/v(k);

end

v(n)=gamma0-theta(n,n)^2*v0-(theta(n,1:m).^2)*v(1:n-1)';

end

for i=186:196

y1(i-185)=theta(i-1,1)*e(i-1);

a(i-185)=data(i);

X1(i-185)=y1(i-185)+a(i-185);

end

%预测

y(197)=theta(196,1)*e(196);

data(198)=y(197)+data(197);

gamma(197)=sum((y(1:0)-mu).*(y(1+197:196)-mu))/197;

theta(197,197)=gamma(197)/v0;

for m=1:197-1

k=197-m;

theta(197,m)=(gamma(m)-theta(k,k)*theta(197,197)*v0-(theta(k,1:k-1).*theta(197,m+1:197-1))*f lipud(v(1:k-1)'))/v(k);

end

e(197)=-b*e(196)+y(197);

y(198)=theta(197,1)*e(197);

data(199)=y(198)+data(198);

data(198:199)

y(197)=theta(196,1)*e(196);

data(198)=y(197)+data(197);

gamma(197)=sum((y(1:0)-mu).*(y(1+197:196)-mu))/197;

theta(197,197)=gamma(197)/v0;

for m=1:197-1

k=197-m;

theta(197,m)=(gamma(m)-theta(k,k)*theta(197,197)*v0-(theta(k,1:k-1).*theta(197,m+1:197-1))*f lipud(v(1:k-1)'))/v(k);

end

e(197)=-b*e(196)+y(197);

y(198)=theta(197,1)*e(197);

data(199)=y(198)+data(198);

data(198:199)