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,-100,1000,700])
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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%比較模型擬合情況(去除趨勢項和季節項)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%最小二乘法估計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
ar=zx(X,Y);%調用最小二乘法函數計算AR(p)的參數

res=zeros(len-p,1);
for i=p+1:1:len
    res(i-p)=ar'*flipud(y2(i-p:i-1,:) );%用AR模型擬合隨機項
end

errArraylist=[];%誤差列表
errSum=0;%誤差和
for i=p+1:1:len
    err=abs((y2(i)-res(i-7))/y2(i));
    errArraylist=[errArraylist,err];
    errSum = errSum + err;
end
averageError=errSum/length(errArraylist) %平均預測誤差
i=p+1:1:len;
subplot(3,2,5)
plot(i,errArraylist(i-p))%誤差走勢圖
title('修正數據相對誤差走勢圖')

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

%%%%%%%%%%%%%%%%%%%模型擬合情況%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
predictTrend=trend(p+1:len,:);%預測趨勢項
predictPeriod=period(p+1:len,:);%預測季節項
predictRandom=res;%預測隨機項
predictResult=predictRandom.*predictPeriod.*predictTrend;%預測結果

preErrArraylist=[];%誤差列表
preErrSum=0;%誤差和
for i=p+1:1:len
    preErr=abs((YOrigin(i)-predictResult(i-7))/YOrigin(i));
    preErrArraylist=[preErrArraylist,preErr];
    preErrSum = preErrSum + preErr;
end
preAverageError=preErrSum/length(preErrArraylist) %平均預測誤差
i=p+1:1:len;
subplot(3,2,6)
plot(i,preErrArraylist(i-p))%誤差走勢圖
title('原始數據相對誤差走勢圖')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%








%%%%%%%%%%%%%%%%%%%%%%%%%%%%%預測函數%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


preLen=30;%預測長度(人為規定，可改)
%%%%%%%%%%%%%%%%%%%%%%%%%建立矩陣%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%預測每周的自變量矩陣
m=[len+1:1:len+preLen]';
preWeek=[];
for i=1:1:weeklength
    for k=len+1:1:len+preLen
        m(k-len)=delta(mod(k-i,weeklength));%調用沖激函數delta賦值
    end
    preWeek=[preWeek,m];
end

%預測春節的自變量矩陣
preSpring=zeros(preLen,21);%21 = 農歷1月初一+前5天 + 后15天 
%預測五一,十一的自變量矩陣
preFive1=zeros(preLen,11);%11=五一 + 前5天 + 后5天
preTen1=zeros(preLen,11);%11=十一 + 前5天 + 后5天

%(暫時認為預測中沒有春節,五一，國慶,這幾項數據都是0)
preZ=[Z;preWeek,preSpring,preFive1,preTen1];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%計算%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
prePeriod=preZ*B;%%預測季節項
i=[1:1:len+preLen]';
X=[(i.^0),i,(i.^2)];
preTrend=X*A;%預測趨勢項


preRes=zeros(len-p+preLen,1);%預測數據
preTemp=y2;%暫存數據
for i=p+1:1:len+preLen
    if i<=len
        preRes(i-p)=ar'*flipud(y2(i-p:i-1,:) );%用AR模型擬合隨機項
    else
        preRes(i-p)=ar'*flipud(preTemp(i-p-p:i-p-1,:));%用AR模型預測隨機項
        preTemp=[preTemp;preRes(i-p)];%%%增加的新數據到暫存數據以便預測
    end
end



predictTrend=preTrend(p+1:len+preLen,:);%預測趨勢項
predictPeriod=prePeriod(p+1:len+preLen,:);%預測季節項
predictRandom=preRes;%預測隨機項
predictResult=predictRandom.*predictPeriod.*predictTrend;%預測結果


i=1:1:len;
i1=p+1:1:len;
i2=len:1:len+preLen;

figure(2)
plot(i,YOrigin(i),i1,predictResult(i1-p),'g',i2,predictResult(i2-p),'r')%預測結果和原始數據的比較
title('預測與原始數據比較圖')
legend('原始數據','擬合數據','預測數據')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%