k=1:81;
X=[125,284,201,27,106,195,177,225,152,131,198,95,30,87,...
296,1,267,67,336,310,400,379,23,190,103,191,162,176,152,...
11,229,41,74,49,29,92,148,13,302,91,68,180,489,175,331,...
127,266,4,56,583,620,261,32,188,243,84,66,289,105,636,...
39,49,161,50,126,198,126,119,7,167,112,261,389,251,230,...
16,259,591,653,631,547];
m(1)=0;
%累加运算
Y(1)=X(1);
for i=2:76;
Y(i)=Y(i-1)+X(i);
end
figure(1)
plot(Y)
%二次累加
y(1)=Y(1);
for j=2:76;
y(j)=y(j-1)+Y(j);
end
% figure(2)
for i=1:75
B(i,:)=[y(i),1];
end
for i=2:76
Yn(i-1)=y(i);
end
%最小二乘
b=B'*B;
b=inv(b)*B';
b=b*Yn';
% 计算预测值y1
y1(1)=y(1);
for k=1:81
y1(k+1)=(b(1).^k)*Y(1)+(1-b(1).^k)/(1-b(1))*b(2);
end
for k=1:81
y0(1)=y(1);
y0(k+1)=y1(k+1)-y1(k);
end
% % 方法二
% y1(1)=y(1);
% for k=1:81
% y1(k+1)=b(1)*y(k)+b(2);
% end
% for k=1:81
% y0(k+1)=y1(k+1)-y1(k);
% end
figure(2)
plot(Y,'b');
hold on
plot(y0,'r')
% 计算差值
z=y0(1:76)-Y(1:76);
% AR模型拟合残差
for i=4:76
B1(i,:)=[z(i-3),z(i-2),z(i-1),z(i),1];
end
b1=B1'*B1;
b1=inv(b1)*B1';
b1=b1*z';
% 残差拟合值
z1=B1*b1;
figure(3)
plot(z,'b');
hold on;
plot(z1,'r');
