时间序列的分析——模型的识别与预测

依照某AR模型生成一段数据（1000），同时用另一MA模型生成一段数据（200），合成一段1200长度的数据 1）依赖于这1200个数据的前800个数据，识别这段数据背后的AR模型。 2）在1）的基础上对新数据进行预测，并通过后续的400个数据进行判别（数据模型是否匹配）或者模型的修正。

clear all;
clc;

%从MA模型X（t）=a（t）-1.3a（t-1）+0.4a（t-2）中提取200个样本点
K=200;
a1=randn（1K）;
x1=zeros（1K）;
x1（1）=1;
x1（2）=1;

for i=3:K
    x1（i）=a1（i）-1.3*a1（i-1）+0.4*a1（i-2）;
end

X1=x1;

%从AR模型X（t）+0.8X（t-1）=a（t）中提取1000个样本点
N=1000;
a2=randn（1N）;
x2=zeros（1N）;
x2（1）=0.1;
x2（2）=0.1;
x2（3）=1;

for i=3:N
    x2（i）=1.5*x2（i-1）-0.5*x2（i-2）+a2（i）;
end

%取x2的前800个数据作为模型识别的样本观测值
L=800;
X2=x2（11:L）;

%p1=autocorr（X1L/10）;
% autocorr（X1）;

%计算样本自相关函数和样本偏相关函数
r=zeros（1L）;           %样本协方差矩阵
for k=0:L-1             %求样本协方差
   rk=0;
   for i=1:L-k
      rk=rk+X2（i）*X2（i+k）;
   end
  r（k+1）=（1/L）*rk;
end

p=zeros（1L）;
for i=1:L            %求样本自相关函数
    p（i）=r（i）/r（1）;
end



%ppp=parcorr（X2L/10）;

 M=L/10;
pp=zeros（MM）;           %样本偏相关系数矩阵
pp（11）=p（2）;
p（1）=[];      

for k=1:M-1              %求样本偏相关函数
    s1=0;
    s2=0;
    for j=1:k
        s1=s1+p（k+1-j）*pp（kj）;
        s2=s2+p（j）*pp（kj）;
    end
    pp（k+1k+1）=（p（k+1）-s1）/（1-s2）;
    for j=1:k
    pp（k+1j）=pp（kj）-pp（k+1k+1）*pp（kk+1-j）;
    end
end  

pp1=diag（pp）;        %提取偏相关系数

subplot（211）
plot（p）;
title（‘自相关图‘）;
subplot（212）
plot（pp1）;
title（‘偏相关图‘）;


%识别模型，返回模型的阶数
%取M=L/10个自相关和偏相关系数

count=0;
for f1=1:M
    if abs（p（f1））<=1/sqrt（L）
        count=count+1;
    end
    if  count/M>=0.683
       break;
    end
end


count_1=0;
for f2=1:M
    if abs（pp1（f2））<=1/sqrt（L）
        count_1=count_1+1;
    end
    if （count_1）/M>=0.683
        break
    end
end


if f1==M&&f2    disp（‘建立AR模型‘）;
    
    %确定模型阶数
     for i=1:M                       
        if  abs（pp1（i））<=1/sqrt（L）
            jieshu=i-1;
            break
        end
    end
    fprintf（‘其阶数为:%d\n‘jieshu）;
    
    %fai=ar（X2jieshu）;             %AR模型求系数函数
    
    %最小二乘法求模型系数
    for n=1:jieshu                   
        Y=X2‘;
        Y（1:n）=[];
        m=L-n;X=[];
        for i=1:m
            for j=1:n
                X（ij）=X2（n+i-j）;
            end
        end
        xishu=inv（X‘*X）*（X‘*Y）;
    end
    
    %用后200个数据检验所得模型是否正确
    X3=x2（1L+1:end）;   %提取后200个样本点
    T=N-L;
    a3=zeros（1T）;        %由观测样本点计算得到的白噪声点
    ra=ones（1T+1）;      %白噪声的协方差
    
    a3（1）=X3（1）;
    a3（2）=X3（2）-xishu（1）*X3（1）;
    
   if jieshu>1     
    for i=1:T