程序.docx
- 文档编号:7381431
- 上传时间:2023-01-23
- 格式:DOCX
- 页数:12
- 大小:45.26KB
程序.docx
《程序.docx》由会员分享,可在线阅读,更多相关《程序.docx(12页珍藏版)》请在冰豆网上搜索。
程序
1平稳性检验方法
例1美国某城镇1964~1983年间每年发生的偷窃率(每万人平均发生的偷盗次
数)的数据见表12(时间t:
1~20表示从1964年~1983年)。
问偷盗率数据是否可看
成是平稳的?
美国某城镇年偷盗率数据
时间12345678910
偷盗率1.372.961.913.102.082.544.073.622.911.94
时间11121314151617181920
偷盗率3.964.192.713.423.023.542.664.114.253.76
程序如下:
x0=[1.372.691.913.102.082.544.073.622.911.943.964.192.713.423.023.542.664.114.253.76];
x0=x0';x0=x0(:
)';n=length(x0);
[Xsort,ind]=sort(x0)
Rt(ind)=1:
n
t=1:
n;
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;
end
rho=gama/gama0
rho=
Columns1through9
0.52270.41610.32760.25400.20310.16160.11770.09750.0822
Column10
0.0641
>>rh02=autocorr(x)
rh02=
Columns1through9
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)
t=1:
n;
bar(t,rho)
n=length(rho2)
t=1:
n;
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所述时间序列的偏相关特性
程序如下:
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;
end
rho=gama/gama0
f(1,1)=rho
(1);
fork=2:
L
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);
end
>>f(k,k)=s1/s2;
>>forj=1:
k-1
f(k,j)=f(k-1,j)-f(k,k)*f(k-1,k-j);
end
end
|
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=parcorr(x)
pcorr2=
Columns1through9
1.00000.51960.20280.07070.03820.04160.02560.0130-0.0196
Columns10through18
0.00920.0149-0.0029-0.01680.0047-0.0132-0.0035-0.00700.0018
Columns19through21
-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:
3
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);
end
end
r=input('R=?
’);
m=input('M=?
’);
spec2=garchset('R',r,'M',m,'Display','off');
[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);
end
m2=length(y);
w=diff(y);
m3=length(w);
fori=0:
3
forj=0:
s+1
spec=garchset('R',i,'M',j,'Display','off');
[coeffX,errorsX,LLFX]=garchfit(spec,w);
num=garchcount(coeffX);
[aic,bic]=aicbic(LLFX,num,m3);
fprintf('R=%d,M=%d,AIC=%f,BIC=%f\n',i,j,aic,bic);
end
end
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程序如下:
clc,clear
loadbdata
spec2=garchset('R',0,'M',13,'Display','off');
[coeffX,errorsX,LLFX]=garchfit(spec2,w);
[sigmaForecast,w_Forecast]=garchpred(coeffX,w,n);
yhat=y(m2)+cumsum(w_Forecast);
forj=1:
n
x(m1+j)=yhat(j)+x(m1+j-s);
end
x_hat=x(m1+1:
end);
x0_hat=exp(x_hat+sigmaForecast.^2/2)
x0_hat_check=exp(x_hat)
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 程序