初次接触matlab及时间序列还请给位前辈多多指教！
我写的代码如下，可以正常运行。主要功能是平稳性检验及平稳化处理，求自相关函数、偏相关函数。AIC准则定阶。
现在的问题是在AIC定阶时我不知道结果是什么，怀疑是数据预处理没有完成，还望各位前辈多多指教。代码如下(有点长)；
clc,clear;
%function()
%+++++++++++++++++++++++++++实验数据+++++++++++++++++++++++++++++++
shiyan01=[1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18;
-0.7 -0.7 -1.7 -1.7 -2.1 -3.1 -3.6 -4.6 -5.3 -5.7 -6.5 -6.0 -6.1 -6.4 -6.7 -6.8 -6.9 -7.1];
tshiyan01=1:0.5:18; %内插期数
shiyan02=[1 10 19 24 27 35 41 46 54 61 69 76;
5.38 5.15 6.69 8.30 8.12 11.05 13.04 11.9 16.26 17.56 18.61 18.03];
tshiyan02=1:3:76; %内插期数
%+++++++++++++++++++++++++++应力数据+++++++++++++++++++++++++++++
yl0201=[1 13 16 20 22 24 26 31 34 36 43 45 49 53 57 61 65 66 70 71 72 73 77 78 79 85 94 101 104 114 120 129 136 139 148 155 161 ;
0.00000000 -4.16848320 -6.94561600 0.34752440 -0.34750120 5.91107000 6.25898880 12.17715800 22.28833280 -7.98665800 9.04317440 101.04583680 127.31974400 121.27260600 127.67566680 141.93160280 163.02686400 172.70756280 173.42532440 174.50214080 176.65640000 178.45225400 180.60804440 180.60804440 180.96742400 207.62589440 181.68625280 323.99977400 194.28077880 197.52396480 205.09870680 211.23813440 210.87680600 204.37686200 215.57588480 217.38426680 214.49113400 ];
tyl=1:3:161; %内插期数
%+++++++++++++++++++++++++++数据赋值+++++++++++++++++++++++++++++++++++
shuju=yl0201;
qs=tyl;
%++++++++++++++++++++++++++插值处理+++++++++++++++++++++++++++++++++++++
q=shuju(1,:);
d=shuju(2,:);
cz=interp1(q,d,qs,'cubic');%%三次插值
figure(6);
subplot(2,1,1)
plot(q,d,'g'),title('原始数据');
subplot(2,1,2)
plot(qs,cz,'r'),title('内插后数据');
M=length(cz);
%+++++++++++++++++++++++平稳化处理及检验++++++++++++++++++++++++++++++++++
jy=cz;
jy=jy';jy=jy(:)';
l=length(jy);
[xsort,ind]=sort(jy);
Rt(ind)=1:l;
t=1:l;
Qs=1-6/(l.*(l^2-1)).*sum((t-Rt).^2);
T=Qs.*sqrt(l-2)/sqrt(1-Qs^2);
ta=tinv(0.975,l-2);
fprintf('平稳性检验结果；\n');
fprintf(' 统计量T 下QS 统计量T的临界值ta \n');
fprintf('==================================================\n');
fprintf('%11.6f%11.6f%11.6f %11.6f\n\n',T,Qs,ta);
if abs(T)>ta %若数据不平稳进行平稳化处理
figure(1);
plot(qs,cz,'r'),title('原始数据');
fprintf('\n数据是非平稳数据，需要平稳化处理。\n')
p1=input('请选择曲线拟合阶次：');
xishu=polyfit(qs,cz,p1 ); %计算拟合曲线系数
nhqx=polyval(xishu,qs,p1); %拟合曲线
pw=cz-nhqx; %求残差
%+++++++++++++++++++++++对平稳后的时间序列进行平稳性检验+++++++++++++
cc=pw; %以防在进行平稳检验是对数列有改变
cc=cc';cc=cc(:)';
l=length(cc);
[xsort,ind]=sort(cc);
Rt(ind)=1:l;
t=1:l;
Qs=1-6/(l.*(l^2-1)).*sum((t-Rt).^2);
T=Qs.*sqrt(l-2)/sqrt(1-Qs^2);
ta=tinv(0.975,l-2);
fprintf('平稳性检验结果；\n');
fprintf(' 统计量T 下QS 统计量T的临界值ta \n');
fprintf('==================================================\n');
fprintf('%11.6f%11.6f%11.6f %11.6f\n\n',T,Qs,ta);
if abs(T)>ta
fprintf('数据仍不平稳，需从新做平稳化处理\n');
else
fprintf('数据平稳。\n');
end
else
pw=cz;
fprintf('\n数据时平稳数据。\n');
end
figure(7);
plot(qs,pw,'g'),title('平稳处理后数据');
cancha=pw(1:M-5);
jmqs=qs(1:M-5);
l1=length(cancha);
%+++++++++++++++++++++++++++++++++求自相关函数+++++++++++++++++++++++++++++
zxghs=autocorr(cancha,[],2,2); %残差数据为三十期
%autocorr()函数是时间序列自相关函数
%sm : 一个时间序列数据
%[]: 表示计算这个时间序列数据的自相关函数的延迟.
%2: 表示自相关函数在>2的所有延迟的自相关系数看作为0
%2: 表示在什么范围内时间序列的数据被看作是0
%++++++++++++++++++++++++++++++++++求偏相关函数++++++++++++++++++++++++++++
pxghs=parcorr(cancha,[],2,2);
%++++++++根据自相关函数和偏相关函数的拖尾截尾情况选择模型++++++++++++++++++++++
p1=length(zxghs);
p2=[1:1:p1];
figure(2);
plot(p2,zxghs,'r',p2,pxghs,'g'),legend('自相关函数','偏相关函数'),title('模型判断');
%+++++++++++++++++++++采用FPE定阶准则定阶++++++++++++++++++++++++++++++++++
N=length(cancha);
%++++++++++++++++++++++++++++AIC准则定阶+++++++++++++++++++++++++++++++++++
e=5; %AIC准则定阶的最高阶
t=0;
for i=0:e
for j=0:e
spec=garchset('R',i,'M',j,'Display','off');
[coeffX,errorsX,LLFX]=garchfit(spec,cancha);
num=garchcount(coeffX);
[aic,bic]=aicbic(LLFX,num,10000);
t=t+1;
a(t)=aic;
b(t)=bic;
fprintf('R=%d.M=%d.AIC=%f.BIC=%f\n',i,j,aic,bic);
end
end
d=1:1:(e+1)^2;
figure(222);
plot(d,a,'r',d,b,'g'),legend('AIC','BIC'),title('AIC准则阶次判断');
fprintf('请依据AIC、BIC值选择(p、q)\n')
p=input('AR模型阶次：');
q=input('MA模型阶次：');
spec=garchset('R',p,'M',q,'Display','off');
[coeffX,errorsX,LLFX]=garchfit(spec,cancha);%计算拟合参数
%e=zeros(1,l1);
e(1)=0;
z=0;
%cancha(0)=0;
%y(1)=cancha(1);
%++++++++++++++++++++++++MA模型++++++++++++++++++++++++++++++++++
if p==0
if q>0
if q==1
for j=q+1:l1 %%%q为参与建模的期数
e(j)=cancha(j)-coeffX.MA*e(j-1);
MA(j)=e(j)-coeffX.MA*e(j-1);
end
end
if q==2
e(1:q)=0;
for j=q+1:l1 %%%q为参与建模的期数
e(j)=cancha(j)-coeffX.MA(1)*e(j-1)-coeffX.MA(2)*e(j-2);
MA(j)=e(j)-coeffX.MA(1)*e(j-1)-coeffX.MA(2)*e(j-2);
end
end
if q==3
e(1:q)=0;
for j=q+1:l1 %%%q为参与建模的期数
e(j)=cancha(j)-coeffX.MA(1)*e(j-1)-coeffX.MA(2)*e(j-2)-coeffX.MA(3)*e(j-3);
MA(j)=e(j)-coeffX.MA(1)*e(j-1)-coeffX.MA(2)*e(j-2)-coeffX.MA(3)*e(j-3);
end
end
if q==4
e(1:q)=0;
for j=q+1:l1 %%%q为参与建模的期数
e(j)=cancha(j)-coeffX.MA(1)*e(j-1)-coeffX.MA(2)*e(j-2)-coeffX.MA(3)*e(j-3)-coeffX.MA(4)*e(j-4);
MA(j)=e(j)-coeffX.MA(1)*e(j-1)-coeffX.MA(2)*e(j-2)-coeffX.MA(3)*e(j-3)-coeffX.MA(4)*e(j-4);
end
end
end
end
%++++++++++++++++++++++++AR模型++++++++++++++++++++++++++++++++++
if q==0
if p>0
for j=p+1:l1 %%%a：q中a=2，q为参与建模的期数
%y(j)=polyval(coeffX.MA,e(j-q:j-1),q)
for i=1:p
z=coeffX.AR(i)*cancha(j-i)+z;
end
AR(j)=z;
e(j)=cancha(j)-z;%计算残差
end
end
end
%++++++++++++++++++++++++ARMA模型++++++++++++++++++++++++++++++++++
if q>0
if p>0
%for j=2:l1 %%%a：q中a=2，q为参与建模的期数
%y(j)=coeffX.AR*cancha(j-1)-coeffX.MA*e(j-1);
% e(j)=cancha(j)-polyval(coeffX.AR,cancha(j-p:j-1),p)-polyval(coeffX.MA,e(j-q:j-1),q);%计算残差
% end
%++++++++++++++++AR模型部分++++++++++++++++++++++++++++++++++
for j=p+1:l1 %%%a：q中a=2，q为参与建模的期数
%y(j)=polyval(coeffX.MA,e(j-q:j-1),q)
for i=1:p
AR(j)=coeffX.AR(i)*cancha(j-i);
end
%y(j)=z(j);
e(j)=cancha(j)-AR(j);%计算残差
end
%+++++++++++++++++++MA模型部分+++++++++++++++++++++++++++
if q==1
for j=q+1:l1 %%
e(j)=cancha(j)-coeffX.MA*e(j-1);
MA(j)=e(j)-coeffX.MA*e(j-1);
end
end
if q==2
e(1:q)=0;
for j=q+1:l1 %%%
e(j)=cancha(j)-coeffX.MA(1)*e(j-1)-coeffX.MA(2)*e(j-2);
MA(j)=e(j)-coeffX.MA(1)*e(j-1)-coeffX.MA(2)*e(j-2);
end
end
if q==3
e(1:q)=0;
for j=q+1:l1 %%%
e(j)=cancha(j)-coeffX.MA(1)*e(j-1)-coeffX.MA(2)*e(j-2)-coeffX.MA(3)*e(j-3);
MA(j)=e(j)-coeffX.MA(1)*e(j-1)-coeffX.MA(2)*e(j-2)-coeffX.MA(3)*e(j-3);
end
end
if q==4
e(1:q)=0;
for j=q+1:l1 %%%
e(j)=cancha(j)-coeffX.MA(1)*e(j-1)-coeffX.MA(2)*e(j-2)-coeffX.MA(3)*e(j-3)-coeffX.MA(4)*e(j-4);
MA(j)=e(j)-coeffX.MA(1)*e(j-1)-coeffX.MA(2)*e(j-2)-coeffX.MA(3)*e(j-3)-coeffX.MA(4)*e(j-4);
end
end
%++++++++++++++++++++++++++++++合成+++++++++++++++++++++++++++++++
if p>q
a=q;
else
a=p;
end
for j=a:l1
ARMA(j)=AR(j)+MA(j)-e(j);
end
end
end
figure(11);
plot(jmqs,e,'b'),title('建模后残差');
%++++++++++++++++++++++++++++++++五步预报++++++++++++++++++++++++++++++++++
[sigmaForecast,x_Forecast]=garchpred(coeffX,cancha,5);%预报，进行5步预报值，其中第一个返回标准差，第二个返回预报值
ybqs=qs(M-4:M);
nhqx1=polyval(xishu,ybqs,p1);
ybz=x_Forecast+nhqx1;
ybz
ybc=cz(M-4:M)-ybz;
ybc
figure(55);
subplot(2,1,1)
plot(ybqs,cz(M-4:M),'g',qs(M-4:M),ybz,'r'),title('多步预报值与插值数据'),legend('后五期插值','预报值');
subplot(2,1,2)
plot(ybqs,ybc,'r'),title('预报误差');
zdwc=max(abs(ybc));%求预测最大误差
zwc=sqrt(sum(ybc.*ybc)/5);%求预测中误差
fprintf('最大预报误差为：%11.2f\n',zdwc);
fprintf('预 报中误差为：%11.2f\n\n',zwc);
3992
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
