程序Word文档格式.docx
《程序Word文档格式.docx》由会员分享,可在线阅读,更多相关《程序Word文档格式.docx(12页珍藏版)》请在冰豆网上搜索。
![程序Word文档格式.docx](https://file1.bdocx.com/fileroot1/2023-1/23/ef0486a7-29ed-49e8-98fd-75f44e9e4ebb/ef0486a7-29ed-49e8-98fd-75f44e9e4ebb1.gif)
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)