work3.m

matlab时间序列预测

clear all
clc

%%%%%%%%求趋势项%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
test=datamain(1);
len = test(2);%原始数据个数
i=[1:1:len]';
X=[(i.^0),i,(i.^2)];
Y=[];
for tm=1:1:len
    test=datamain(tm);
    Y=[Y;test(1)];%读取原始数据
end
YOrigin=Y;%保留原始数据
figure('Position',[10,50,1000,600])
subplot(3,2,1)
plot(i,YOrigin(i))
title('原始数据走势图')
A=ZX(X,Y);%用最小二乘法求趋势项参数
trend=X*A;%求出趋势项
subplot(3,2,2)
plot(i,trend(i))
title('趋势项走势图')
y1=[1:1:len]';
for i=1:1:len;
    temp=datamain(i);
    y1(i)=temp(1)/trend(i);%y1=原始数据/趋势项数据
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%计算季节项%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%每周周期的自变量矩阵
weeklength = 7 ;%周期长度
m=[1:1:len]';
week=[];
for i=1:1:weeklength
    for k=1:1:len
        m(k)=delta(mod(k-i,weeklength));%调用冲激函数delta赋值
    end
    week=[week,m];
end

%春节的自变量矩阵
spring=zeros(len,21);%21 = 农历1月初一+前5天 + 后15天 
%根据样本的春节数据设定自变量矩阵
for i=0:1:20
    spring(17+i,1+i)=1;%样本数据的第1个春节
    week(17+i,:)=zeros(1,7);%消除week周期的叠加影响
    spring(401+i,1+i)=1;%样本数据的第2个春节
    week(401+i,:)=zeros(1,7);%消除week周期的叠加影响
    spring(755+i,1+i)=1;%样本数据的第3个春节
    week(755+i,:)=zeros(1,7);%消除week周期的叠加影响
end

%五一,十一的计算矩阵
five1=zeros(len,11);%11=五一 + 前5天 + 后5天
ten1=zeros(len,11);%11=十一 + 前5天 + 后5天
for i=0:1:10
    five1(117+i,1+i)=1;%样本数据的第1个五一
    week(117+i,:)=zeros(1,7);%消除week周期的叠加影响
    five1(482+i,1+i)=1;%样本数据的第2个五一
    week(482+i,:)=zeros(1,7);%消除week周期的叠加影响
    five1(847+i,1+i)=1;%样本数据的第3个五一
    week(847+i,:)=zeros(1,7);%消除week周期的叠加影响
    ten1(270+i,1+i)=1;%样本数据的第1个十一
    week(270+i,:)=zeros(1,7);%消除week周期的叠加影响
    ten1(635+i,1+i)=1;%样本数据的第2个十一
    week(635+i,:)=zeros(1,7);%消除week周期的叠加影响
    ten1(998+i,1+i)=1;%样本数据的第3个十一
    week(998+i,:)=zeros(1,7);%消除week周期的叠加影响
end


%求季节项
Z=[week,spring,five1,ten1];%综合季节项自变量矩阵(7周期，春节，五一，十一)
B=ZX(Z,y1);%调用最小二乘法计算综合季节项参数
period=Z*B;%计算季节项
subplot(3,2,3)
i=1:1:len;
plot(i,period(i))
title('季节项项走势图')
y2=[1:1:len]';
for i=1:1:len;
y2(i)=y1(i)/period(i);%为最后修正的数据y2=y1/季节项
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%计算AIC(p),BIC(p)%%%%%%%%%%%%%%%%%%%%%%%%%%%
startNumber=1;%开始搜索
endNumber=15;%结束搜索
AICM=[];
BICM=[];
 
for p=startNumber:1:endNumber

    %最小二乘法估计AR(p)的参数    
    Y=y2(p+1:len,:);%建立Y矩阵
    X=[];
    for k=p:-1:1
        X=[X,y2(k:len-1-p+k,:)];%建立参数的自变量矩阵(len-1-p+k)
    end
    a=zx(X,Y);%调用最小二乘法函数计算AR(p)的参数

    %白噪声方差的最小二乘估计
    S=0;
    for i=p+1:1:len
        S=S+(y2(i)-a'*flipud(y2(i-p:i-1,:) ))^2;%残差平方和(flipud)
    end
    delta2=S/(len-p);%除以残差个数=白噪声方差

    %AIC准则
    aicVar=log(delta2)+2*p/len;

    %BIC准则 
    bicVar=log(delta2)+p*log(len)/len;

    AICM=[AICM,aicVar];      
    BICM=[BICM,bicVar];
 end

p=startNumber:1:endNumber;
subplot(3,2,4)
plot(p,AICM(p),'r',p,BICM(p),'g')%AIC(p)和BIC(p)p=startNumbe到endNumber的值作图
title('AIC(p),BIC(p)值图')
legend('AIC(p)','BIC(p)')

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

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

p=7;%人为判断最合理的p值为7
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%比较模型拟合情况(去除趋势项和季节项)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%使用matlab的arx函数计算AR(7)模型参数


ar=arx(y2,7)
res=predict(y2,ar);

errArraylist=[];%误差列表
errSum=0;%误差和
for i=1+p:1:len
    err=abs((y2(i)-res(i))/y2(i));
    errArraylist=[errArraylist,err];
    errSum = errSum + err;
end
averageError=errSum/length(errArraylist) %平均预测误差
i=1:1:len-p;
subplot(3,2,5)
plot(i,errArraylist(i))%误差走势图
title('修正数据相对误差走势图')

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

%%%%%%%%%%%%%%%%%%%模型拟合情况%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
predictTrend=trend(1:len,:);%预测趋势项
predictPeriod=period(1:len,:);%预测季节项
predictRandom=res;%预测随机项
predictResult=predictRandom.*predictPeriod.*predictTrend;%预测结果

preErrArraylist=[];%误差列表
preErrSum=0;%误差和
for i=1+p:1:len
    preErr=abs((YOrigin(i)-predictResult(i))/YOrigin(i));
    preErrArraylist=[preErrArraylist,preErr];
    preErrSum = preErrSum + preErr;
end
preAverageError=preErrSum/length(preErrArraylist) %平均预测误差
i=1:1:len-p;
subplot(3,2,6)
plot(i,preErrArraylist(i))%误差走势图
title('原始数据相对误差走势图')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%