hgyc.m

线性回归算法,其功能是对一组数据进行分析预测,数据分为年度,月度,季度数据

function [yh,me]=hgyc(aa,bz)
% input %aa第一列为时间，年度数据直接为相应年份；
                     %月度数据或季度数据的格式为六位整数--前四位为年份，后两位为月份或季度
        %aa第二列为统计指标的取值
        %bz 是数据类型的标志，对于月度数据、季度数据，统计指标为累计值的时候bz=1,否则bz=0
%aa=xlsread('D:\jjyj\datav\example.xls','sheet1');
nn=size(aa,1);
tls=1:nn;
ts=tls';
tt=aa(:,1);
yt0=aa(:,2);
yt=yt0;

%判定数据的频率
if tt<10000
    tnum=1;
else
    ms=tt-floor(tt/100)*100;
    tnum=max(ms);
end
if bz==1
    tt1=floor(tt/100);
    bn=min(tt1);
    en=max(tt1);
    for yr=bn:en
        ly=find(floor(tt/100)==yr);
        ny=size(ly,1);
        yt(ly(2:ny))=diff(yt(ly));
    end
end


if tnum==1  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%计算时间趋势
cc=ones(nn,1);
[b1,bint1,r1,rint1,stats1] = regress(yt,[cc ts]);
SIC1=log(sum(r1.*r1)/nn)+2*log(nn)/nn;

[b2,bint2,r2,rint2,stats2] = regress(yt,[cc ts ts.*ts]);
SIC2=log(sum(r2.*r2)/nn)+3*log(nn)/nn;

[b3,bint3,rc,rint3,stats3] = regress(log(yt),[cc ts]);
r3=yt-exp(log(yt)-rc);
SIC3=log(sum(r3.*r3)/nn)+2*log(nn)/nn;

if SIC1<SIC2&SIC1<SIC3&stats1(2)>finv(0.95,1,nn-2)
    xf0=[ones(5,1) (ts(nn)+1:ts(nn)+5)'];
    md=1;
    XX0=[cc ts];
    k=2;
end
if SIC2<SIC1&SIC2<SIC3&stats2(2)>finv(0.95,2,nn-3)
    xf0=[ones(5,1) (ts(nn)+1:ts(nn)+5)' (ts(nn)+1:ts(nn)+5)'.*(ts(nn)+1:ts(nn)+5)'];
    md=2;
    XX0=[cc ts ts.*ts];
    k=3;
end
if SIC3<SIC1&SIC3<SIC2&stats3(2)>finv(0.95,1,nn-2)
    xf0=[ones(5,1) (ts(nn)+1:ts(nn)+5)']
    md=3;
    XX0=[cc ts];
    k=2;
    yt=log(yt);
end

if stats1(2)<=finv(0.95,1,nn-2)&stats2(2)<=finv(0.95,2,nn-3)&stats3(2)<=finv(0.95,1,nn-2)
    xf0=ones(5,1);
    md=0;
    XX0=cc;
    k=1;  
end

[b1,bint1,r1,rint1,stats1] = regress(yt,XX0);
sic=log(sum(r1.*r1)/nn)+(k+1)*log(nn)/nn;
for j=1:5
    XX0=[XX0 lagmatrix(yt,j)];
    [b1,bint1,r1,rint1,stats1] = regress(yt(j+1:nn),XX0(j+1:nn,:));
    SIC1=log(sum(r1.*r1)/(nn-j))+(k+1+j)*log(nn-j)/(nn-j);
    sic=[sic SIC1];  
end
[sc,ln]=min(sic);
XX2=XX0(:,1:k+ln-1);
kk=size(XX2,2);
XXE=XX2;

[be,binte,re,rinte,statse] = regress(yt,XXE);
if md==3
    mape=mean(abs(exp(re(ln:nn))./exp(yt(ln:nn))))*100;
else
    mape=mean(abs(re(ln:nn)./yt(ln:nn)))*100;
end

yht=[];
for j=1:5
     xfe=xf0(j,:);
    if ln>1
        for m=1:ln-1
            xfe=[xfe yt(nn-m+j)];
        end
    end
    yls=xfe*be;
    yt=[yt;yls];
    if md==3
        yht=[yht;exp(yls)];
    else
        yht=[yht;yls];
    end
end
else %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
DD=eye(tnum,tnum);
ts0=(ts(nn)+1:ts(nn)+tnum)';
D0=zeros(nn,tnum);
for i=1:nn
    ls=tt(i)-floor(tt(i)/100)*100;
    D0(i,:)=DD(ls,:);
end
DD0=[];
for k=1:tnum
    tl=tt(nn)-floor(tt(nn)/100)*100+k;
    if tl>tnum
        tl=tl-tnum;
    end
    DD0=[DD0;DD(tl,:)];
end
        
%计算时间趋势
cc=ones(nn,1);
[b1,bint1,r1,rint1,stats1] = regress(yt,[cc ts]);
SIC1=log(sum(r1.*r1)/nn)+2*log(nn)/nn;

[b2,bint2,r2,rint2,stats2] = regress(yt,[cc ts ts.*ts]);
SIC2=log(sum(r2.*r2)/nn)+3*log(nn)/nn;

[b3,bint3,rc,rint3,stats3] = regress(log(yt),[cc ts]);
r3=yt-exp(log(yt)-rc);
SIC3=log(sum(r3.*r3)/nn)+2*log(nn)/nn;

if SIC1<SIC2&SIC1<SIC3&stats1(2)>finv(0.95,1,nn-2)
    xf0=[ts0 DD0];
    md=1;
    XX0=[ts D0];
    k=tnum+1;
end
if SIC2<SIC1&SIC2<SIC3&stats2(2)>finv(0.95,2,nn-3)
    xf0=[ts0 ts0.*ts0 DD0];
    md=2;
    XX0=[ts ts.*ts D0];
    k=tnum+2;
end
if SIC3<SIC1&SIC3<SIC2&stats3(2)>finv(0.95,1,nn-2)
    xf0=[ts0 DD0];
    md=3;
    XX0=[ts D0];
    k=tnum+1;
    yt=log(yt);
end

if stats1(2)<=finv(0.95,1,nn-2)&stats2(2)<=finv(0.95,2,nn-3)&stats3(2)<=finv(0.95,1,nn-2)
    xf0=DD0;
    md=0;
    XX0=D0;
    k=tnum;  
end
[b0,bint0,r0,rint0,stats0] = regress(yt,XX0);
sigma=sqrt(sum(r0.*r0)/(nn-k));
tstats=b0./(sigma*sqrt(diag(inv(XX0'*XX0))));
lst=find(abs(tstats)>tinv(0.95,nn-k));
XX1=XX0(:,lst');
xf1=xf0(:,lst');
ne=size(XX1,2);

[b1,bint1,r1,rint1,stats1] = regress(yt,XX1);
sic=log(sum(r1.*r1)/nn)+ne*log(nn)/nn;
pLB=0;
for j=1:5
   rls=lagmatrix(r1,j);
   xg=corrcoef(r1(j+1:nn),rls(j+1:nn));
   pLB=pLB+xg(1,2)*xg(1,2)/(nn-j);
end
LB=pLB*nn*(nn+2);

if LB>chi2inv(0.95,5)&nn>=3*(ne+1)%@@@@@@@@@
for j=1:5
    XX1=[XX1 lagmatrix(yt,j)];
    xf1=[xf1 yt(nn-j+1-tnum+1:nn-j+1)];
    [b1,bint1,r1,rint1,stats1] = regress(yt(j+1:nn),XX1(j+1:nn,:));
    SIC1=log(sum(r1.*r1)/(nn-j))+(ne+j)*log(nn-j)/(nn-j);
    sic=[sic SIC1];  
end
[sc,ln]=min(sic);
XX2=XX1(:,1:ne+ln-1);
xf2=xf1(:,1:ne+ln-1)
k=size(XX2,2);

[b0,bint0,r0,rint0,stats0] = regress(yt(ln:nn),XX2(ln:nn,:));
sigma=sqrt(sum(r0.*r0)/(nn-ne-ln+1));
tstats=b0./(sigma*sqrt(diag(inv(XX2(ln:nn,:)'*XX2(ln:nn,:)))));
lst=find(abs(tstats)>tinv(0.95,nn-ln+1-k))';
ls=find(lst<=ne);
xfe=xf2(:,lst);
XXE=XX2(:,lst);
if  size(ls,2)==0
    XXE=[cc XXE];
    xfe=[ones(tnum,1) xfe];
end

yht=[];
[be,binte,re,rinte,statse] = regress(yt,XXE);
if md==3
    mape=mean(abs(exp(re(ln:nn))./exp(yt(ln:nn))))*100;
else
    mape=mean(abs(re(ln:nn)./yt(ln:nn)))*100;
end
for p=1:tnum
    yls=xfe(p,:)*be;
    if md==3
        yls=exp(yls);
    end
    yht=[yht;yls];
end
else %@@@@@@@@@@
if md==3
    mape=mean(abs(exp(r1)./exp(yt)))*100;
else
    mape=mean(abs(r1./yt))*100;
end
yht=[];
for p=1:tnum
    yls=xf1(p,:)*b1;
    if md==3
        yls=exp(yls);
    end
    yht=[yht;yls];
end
    
end %@@@@@@@@@@    

end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if bz==1
    tf=[];
    t0=tt(nn);
    ty=floor(t0/100);
    tm=t0-floor(t0/100)*100;
    for t=1:tnum
        tm=tm+1;
        if tm>tnum;
            tm=tm-tnum;
            ty=ty+1;
        end
        tf=[tf;ty*100+tm];
    end
    yhh=[yt;yht];
    ttf=[tt;tf];
    ttf1=floor(ttf/100);
    bn=min(ttf1);
    en=max(ttf1);
    for yr=bn:en
        ly=find(floor(ttf/100)==yr);
        yhh(ly)=cumsum(yhh(ly));
    end
    yht=yhh(nn+1:nn+tnum);
end
yh=num2str(yht);
me=[num2str(mape) '%'];