周期项之和为0

代码：

%时间序列的典型分析式
%数据来源网络
x=[9007,8106,8928,9137,10017,10826,11317,10744,9713,9938,9161,8927,...
  7750,6981,8038,8422,8714,9512,10120,9823,8743,9129,8710,8680,...
  8162,7306,8124,7870,9387,9556,10093,9620,8285,8433,8160,8034,...
  7717,7461,7776,7925,8634,8945,10078,9179,8037,8488,7874,8647,...
  7792,6957,7726,8106,8890,9299,10625,9302,8314,8850,8265,8796,...
  7836,6892,7791,8129,9115, 9434,10484,9827,9110,9070,8633,9240];
D=[9007,8106,8928,9137,10017,10826,11317,10744,9713,9938,9161,8927;
   7750,6981,8038,8422,8714,9512,10120,9823,8743,9129,8710,8680;
   8162,7306,8124,7870,9387,9556,10093,9620,8285,8433,8160,8034;
   7717,7461,7776,7925,8634, 8945,10078,9179,8037,8488,7874 8647;
   7792,6957,7726,8106,8890,9299,10625,9302,8314,8850,8265,8796;
   7836,6892,7791,8129,9115,9434,10484,9827,9110,9070,8633,9240];%6行12列
st=zeros(6,12);
s=zeros(1,12);
x1=mean(D'); %求每一行平均值
for i=1:6
    for j=1:12
        st(i,j)=D(i,j)-x1(i);
    end
end
q=sum(st,1); %求每一列的和
for j=1:12
    s(j)=q(j)/6;  %周期项
end
y=zeros(72,1);
for i=1:6
    for j=1:12
        k=(i-1)*12+j
        y(k)=D(i,j)-s(j);  %y是消除周期项之后的数据
    end
end

%对消除周期项之后的数据进行拟合
a=zeros(72,2);
for i=1:72       %回归拟合
    a(i,1)=1;
    a(i,2)=i;
end
xi=inv(a'*a)*a'*y;   %两个系数
pre=zeros(1,84);   %直线
for i=1:84
    pre(i)=xi(1)+xi(2)*i;
end
subplot(2,1,1);
plot(1:72,y,1:72,pre(1:72));
jieguo=zeros(1,84);
for i=1:7
    for j=1:12
        k=(i-1)*12+j;
        jieguo(k)=pre(k)+s(j);
    end
end
subplot(2,1,2);
plot(1:72,x,'<',1:84,jieguo);