程序Word文档格式.docx
《程序Word文档格式.docx》由会员分享，可在线阅读，更多相关《程序Word文档格式.docx（12页珍藏版）》请在冰豆网上搜索。
Qs=1-6/（n*（n^2-1））*sum（（t-Rt）.^2）
t=Qs*sqrt（n-2）/sqrt（1-Qs^2）
t_0=tinv（0.975,n-2）
2
2自相关系数与偏相关系数计算程序
例2随机模拟下列序列n=10000
利用模型数据研究上述序列自相关特性
randn（'
state'
sum（clock））;
elps=randn（1,10000）;
x
（1）=0;
forj=2:
10000
x（j）=0.8*x（j-1）+elps（j）-0.4*elps（j-1）;
end
y=（x-mean（x））;
gama0=var（x）;
L=10;
forj=1:
L
gama（j）=y（j+1:
end）*y（1:
end-j）'
/10000;
rho=gama/gama0
rho=
Columns1through9
0.52270.41610.32760.25400.20310.16160.11770.09750.0822
Column10
0.0641
>
rh02=autocorr（x）
rh02=
1.00000.52280.41610.32760.25400.20320.16160.11770.0975
Columns10through18
0.08220.06410.05440.03510.02650.02350.02350.01520.0192
Columns19through21
0.01750.01920.0134
自相关系数图：
n=length（rho）
bar（t,rho）
n=length（rho2）
bar（t,rho2）
计算理论自相关函数的Matlab程序为
rho=@（k,p1,p2）（1-p1*p2）*（p1-p2）/（1+p2^2-2*p1*p2）*p1.^（k-1）;
a=rho（[1:
10],0.8,0.4）
例3讨论例14所述时间序列的偏相关特性
f（1,1）=rho
（1）;
fork=2:
s1=rho（k）;
s2=1
forj=1:
k-1
s1=s1-rho（k-j）*f（k-1,j）;
s2=s2-rho（j）*f（k-1,j）;
f（k,k）=s1/s2;
f（k,j）=f（k-1,j）-f（k,k）*f（k-1,k-j）;
|
Error:
Illegaluseofreservedkeyword"
end"
.
f
f=
0.51950
0.30240.4180
pcorr=diag（f）'
pcorr=
0.51950.4180
pcorr2=parcorr（x）
?
Undefinedfunctionormethod'
parcoor'
forinputargumentsoftype'
double'
pcorr2=
1.00000.51960.20280.07070.03820.04160.02560.0130-0.0196
0.00920.0149-0.0029-0.01680.0047-0.0132-0.0035-0.00700.0018
-0.02370.01460.0015
偏相关系数图：
例2：
某化工生产过程每2小时的浓度读数如表如下所示
化工生产过程浓度数据
17.016.616.316.117.116.916.817.417.117.0
16.717.417.217.417.417.017.317.217.416.8
17.117.417.417.517.417.617.417.317.017.8
17.518.117.517.417.417.117.617.717.417.8
17.617.516.517.817.317.317.117.416.917.3
17.616.916.716.816.817.216.817.617.216.6
17.116.916.618.017.217.317.016.917.316.8
17.317.417.716.816.917.016.917.016.616.7
16.816.716.416.516.416.616.516.716.416.4
16.216.416.316.417.016.917.117.116.716.9
16.517.216.417.017.016.716.216.616.916.5
16.616.617.017.117.116.716.816.316.616.8
16.917.116.817.017.217.317.217.317.217.2
17.516.916.916.917.016.516.716.816.716.7
16.616.517.016.716.716.917.417.117.016.8
17.217.217.417.216.916.817.017.417.217.2
17.117.117.117.417.216.916.917.016.716.9
17.317.817.817.617.517.016.917.117.217.4
17.517.917.017.017.017.217.317.417.417.0
18.018.217.617.817.717.217.4
（1）对这一生产过程建模；
（2）对这一生产过程进行10步预测
A=[17.016.616.316.117.116.916.817.417.117.016.717.417.217.417.417.017.317.217.416.817.117.417.417.517.417.617.417.317.017.817.518.117.517.417.417.117.617.717.417.817.617.516.517.817.317.317.117.416.917.317.616.916.716.816.817.216.817.617.216.617.116.916.618.017.217.317.016.917.316.817.317.417.716.816.917.016.917.016.616.716.816.716.416.516.416.616.516.716.416.416.216.416.316.417.016.917.117.116.716.916.517.216.417.017.016.716.216.616.916.516.616.617.017.117.116.716.816.316.616.816.917.116.817.017.217.317.217.317.217.217.516.916.916.917.016.516.716.816.716.716.616.517.016.716.716.917.417.117.016.817.217.217.417.216.916.817.017.417.217.217.117.117.117.417.216.916.917.016.716.917.317.817.817.617.517.016.917.117.217.417.517.917.017.017.017.217.317.417.417.018.018.217.617.817.717.217.4]
化数据
a=textread（'
hua.txt'
）;
a=nonzeros（a'
r11=autocorr（a）
r12=parcorr（a）
da=diff（a）;
r21=autocorr（da）
r22=parcorr（da）
n=length（da）;
fori=0:
3
forj=0:
spec=garchset（'
R'
i,'
M'
j,'
Display'
'
off'
[coeffX,errorsX,LLFX]=garchfit（spec,da）;
num=garchcount（coeffX）;
[aic,bic]=aicbic（LLFX,num,n）;
fprintf（'
R=%d,M=%d,AIC=%f,BIC=%f\n'
i,j,aic,bic）;
r=input（'
R=?
’）;
m=input（'
M=?
spec2=garchset（'
r,'
m,'
[coeffX,errorsX,LLFX]=garchfit（spec2,da）
[sigmaForecast,w_Forecast]=garchpred（coeffX,da,10）
x_pred=a（end）+cumsum（w_Forecast）
例3国际航线月度旅客总数（1949.1-1960.12单位：
千人）
国际航线月度旅客总数
月份
年份123456789101112
1949112118132129121135148148136119104118
1950115126141135125149170170158133114140
1951145150178163172178199199184162146166
1952171180193181183218230242209191172194
1953196196236235229243264272237211180201
1954204188235227234264302293259229203229
1955242233267269270315364347312274237278
1956284277317313318374413405355306271306
1957315301356348355422465467404347305336
1958340318362348363435491505404359310337
1959360342406396420472548559463407362405
1960417391419461472535622606408461390432
A=[112118132129121135148148136119104118115126141135125149170170158133114140145150178163172178199199184162146166171180193181183218230242209191172194196196236235229243264272237211180201204188235227234264302293259229203229242233267269270315364347312274237278284277317313318374413405355306271306315301356348355422465467404347305336340318362348363435491505404359310337360342406396420472548559463407362405417391419461472535622606408461390432]
对所给时间序列建模，并对时间序列进行2年（24个月）的预报
时间序列图形如下：
从直观上看，这一时间序列有上升趋势，且有明显的12月的季节性，而且波动幅度越来越大，下面采用数据变换的方法使数据平稳，对Xt做对数变换
设Yt为一个乘积型季节性序列，对Yt作差分
用AIC准则对模型定阶，当p=0,q=13时AIC达到最小值-431.5504，由此可见取模型
计算的matlab程序如下
clc,clear
loadhang.txt;
hang=hang'
x0=hang（:
m1=length（x0）;
plot（linspace（49,61,m1）,x0）
x=log（x0）;
s=12;
n=24;
fori=s+1:
m1
y（i-s）=x（i）-x（i-s）;
m2=length（y）;
w=diff（y）;
m3=length（w）;
3
s+1
spec=garchset（'
[coeffX,errorsX,LLFX]=garchfit（spec,w）;
[aic,bic]=aicbic（LLFX,num,m3）;
savebdataxywnm1m2s
（2）建立模型并进行预测。
根据上面确定的模型阶数，建立tW的MA（13）模型，并预报tY的24个月的值。
然而，我们要求的是tX的24个月的预报值。
因为ttY=lnX，并假设Yt是正态序列。
在概率论中有这一事实：
设X是随机变量，又lnX服从正态分布
，则X服从对数正态分布，X的均值是
这个例子中，
是对数正态序列，
代表Yt在时刻k（k=144）的m步预报，
代表预报标准差，
代表Xt在时刻k（k=144）的m步预报值，则有
由此算得
的24个月预报值如下：
466.884410.313458.429495.865508.518580.256685.96636.911495.594
504.109425.486476.18513.984451.762504.6546.089560.091639.183
755.713701.763546.123555.574468.981524.921
计算的matlab程序如下：
loadbdata
spec2=garchset（'
0,'
13,'
[coeffX,errorsX,LLFX]=garchfit（spec2,w）;
[sigmaForecast,w_Forecast]=garchpred（coeffX,w,n）;
yhat=y（m2）+cumsum（w_Forecast）;
x（m1+j）=yhat（j）+x（m1+j-s）;
x_hat=x（m1+1:
end）;
x0_hat=exp（x_hat+sigmaForecast.^2/2）
x0_hat_check=exp（x_hat）
升级会员 