ar_smart.m

第一次上载程序

clear;
%%%AR
% load   Eunitefhgz.mat    %%  Eunitef competition-data 2001
load  Eunitefhgz97.mat
e=0.001;                         %矩阵的系数，一般是小于1的数，但也不宜过小，可能导致计算过程发散
NN=30;                           %采样次数
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%fhysh=fhysh/100;
for h=1:1:48
    for i=1:1:250
      y(h,i)=fhysh((i-1)*48+h);  % 按每天48小时组成矩阵
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for d=1:1:48                     %每天48小时


  
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
  p=5;                           %AR模型的阶数
  start=16;                      %序列的起始点
  Y=0;                           %初始化成0，为了下一次循环进入
  canshu=0;                      %初始化
  P=0;                           %初始化
  P=1/e*eye(p,p);                %初始P矩阵                
  
  for i=1:1:p                    %初始化 canshu(:,1);
    canshu(i,1)=0;
  end
  
  for N=1:1:NN                   % 采样次数
    
    for i=1:1:p
      Y(i,N+1)=y(d,start+(N-1)-i+1);%%% 顺序是从前向后,最后一个数是114，滚动的取值，每循环一次N，多一个
    end
    
    %%%%%%%%%%%%%%%%% ＃＃求最佳参数向量fai＃＃ %%%%%%%%%%%%%%%%%%%%%%%%%%%
    ZZ=1/(1+Y(:,N+1)'*P(:,:,N)*Y(:,N+1));%% 实数
    P(:,:,N+1)=P(:,:,N)-ZZ*P(:,:,N)*Y(:,N+1)*Y(:,N+1)'*P(:,:,N);  %
    canshu(:,N+1)=canshu(:,N)+P(:,:,N+1)*Y(:,N+1)*(y(d,start+(N-1)+p)-Y(:,N+1)'*canshu(:,N));%%Yt
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    for i=1:1:p  
        fai(p,N+1,i)=canshu(i,N+1);%% 记录一下，用于以后预报
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  end%% NN
  
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ＃＃ 预报 ＃＃ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    CC=0;
    DD=0;
    EE=0;

    for i=1:1:p
      CC=CC+fai(p,NN+1,i)*y(d,NN+start-i);  %%接近预报那天的前P个数来预报，不管优化的P是多大，都是预报最大P的后面的三天，而fai是倒着取的，假如fai比P小，就取后面几个数作为fai
    end
    yb(d,NN+start)=CC; 
    for i=2:1:p
      DD=DD+fai(p,NN+1,i)*y(d,NN+start+1-i);
    end
    yb(d,NN+start+1)=DD+fai(p,NN+1,1)*yb(d,NN+start);
    for i=3:1:p
      EE=EE+fai(p,NN+1,i)*y(d,NN+start+2-i);
    end
    yb(d,NN+start+2)=EE+fai(p,NN+1,1)*yb(d,NN+start)+fai(p,NN+1,2)*yb(d,NN+start+1);
   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 
end%% d



%%%%%%%%%%%%%%%%%%%%%%%% ＃＃计算误差＃＃ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
for i=1:1:48
     wucha1(i)=(y(i,NN+start)-yb(i,NN+start))/y(i,NN+start)*100;
end
for i=1:1:48
     wucha2(i)=(y(i,NN+start+1)-yb(i,NN+start+1))/y(i,NN+start+1)*100;
end
for i=1:1:48
     wucha3(i)=(y(i,NN+start+2)-yb(i,NN+start+2))/y(i,NN+start+2)*100;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


wucha=[wucha1 wucha2 wucha3];

i=1:1:48;plot(i,wucha1(i),'.-',i,wucha2(i),'.-',i,wucha3(i),'.-');grid on;