AR自回归模型matlab预测程序

%AR时间序列模型估计

clear
tic
%%%%%%%%%%%%%%%第一步，模拟一个AR模型并绘制ACF，PACF图%%%%%%%%%%%%%
%s首先设定AR模型的多项式系数。AR模型中只有多项式A（q）和C（q），
a1 = 0.6;
a2 = 0.8;
a3 = 0;
a4 = 0;
c1 = 0;
c2 = 0;
c3 = 0;
c4 = 0; 
obv = 700;                        %obv是模拟的观测数目。 
A = [1 a1 a2 a3 a4];
B = [];                           %因为AR模型没有输入，因此多项式B是空的。
C = [1 c1 c2 c3 c4];
D = [];                           %把D也设为空的。
F = [];                           %AR模型里的F多项式也是空的。
m = idpoly（ABCDF11）    %这样就生成了AR模型，把它存储在m中。NoiseVariance被设定为1，1也是默认值。抽样间隔Ts设为1。 
error = randn（obv1）;           %生成一个obv*1的正态随机序列。准备用作模型的误差项。
e = iddata（[]error1）;        %用randn函数生成一个噪声序列。存储在e中。抽样间隔是1秒。                       。
y = sim（me）; 
get（y）                            %使用get函数来查看动态系统的所有性质。
Y=y.OutputData; %把y.OutputData的全部值赋给变量r，r就是一个obv*1的向量。
dlmwrite（‘data.txt‘Y） %将生成的700个数据保存在data.txt中
Y=Y‘
figure（1）
plot（Y）                 %绘出y随时间变化的曲线。 
title（‘时间序列图像‘）
%%%%%%%%%%%%%%%%-检验预处理后的数据是否符合AR建模要求，计算自相关和偏相关系数%%%%%%%%%%%
R0=0;
for i=1:500
  R0=Y（i）^2/500+R0;
end
R0
for k=1:20
  R（k）=0;
  for i=k+1:500
  R（k）=Y（i）*Y（i-k）/500+R（k）;
  end
  R %自协方差函数R   
end
x=R/R0 %自相关系数x
figure;
plot（x）
title（‘自相关系数分析图‘）;  
%%%%%%%%%%%%%%%%%%%%%%%%%%计算自相关系数%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%解Y-W方程，其系数矩阵是Toeplit矩阵。求得偏相关函数X%%%%%%%%%%%
X1=x（1）;
X11=x（1）/1;
B=[x（1） x（2）]‘;
x2=[1 x（1）];
A=toeplitz（x2）;   
X2=A\B
X22=X2（2）

B=[x（1） x（2） x（3）]‘;
x3=[1 x（1） x（2）];
A=toeplitz（x3）;   
X3=A\B
X33=X3（3）

B=[x（1） x（2） x（3） x（4）]‘;
x4=[1 x（1） x（2） x（3）];
A=toeplitz（x4）;   
X4=A\B
X44=X4（4）

B=[x（1） x（2） x（3） x（4） x（5）]‘;
x5=[1 x（1） x（2） x（3） x（4）];
A=toeplitz（x5）;   
X5=A\B
X55=X5（5）
X=[X11 X22 X33 X44 X55 ]  
%%%%%%%%%%%%%%%%%%%%%%%解Y-W方程，得偏相关函数X%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure;  
plot（X）;
grid on
title（‘偏相关函数图‘）;

%%%%%%%%%%%%%%%根据偏相关函数截尾性，初判模型阶次。%%%%%%%%%%%%%%%%%%%
%用最小二乘法估计参数计算5阶以内的模型残差方差和AIC值，应用AIC准则为模型定阶%
  S=[R0 R（1）];
  G=toeplitz（S）;
  W=inv（G）*[R（1:2）]‘           
    
  K=0;
  for t=3:498   
  r=0;  
  for i=1:2
  r=W（i）*Y（t-i）+r;
  end
  at= Y（t）-r;
  K=（at）^2+K;  
  end
  U（2）=K/（500-2） % 2阶模型残差方差 

    
K=0;T=X1;
for t=2:498
  at=Y（t）-T（1）*Y（t-1）;
  K=（at）^2+K;  
end   
  U（1）=K/（4