%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%基于Kalman滤波算法的自适应AR模型%%%%%%%%%%%%%%%%%%%%%%
%%优点:算法收敛速度快;缺点:算法较复杂。
clc;
clear all;
close all;
load data3;
z=A3;
n=length(z);
N=8108;
M=18;%预报步数
M1=30;%最大阶数
n=18;%预报步数,调整n值可得对应步数下的性能指标值
c=1;%c为小的正数
q=1;
avrg=sum(z(1:N))/N;
xg=zeros(N+M,1);
s=zeros(M,1);
a=zeros(M1,1);
for i=1:N;
xg(i)=z(i)-avrg;
end
%定阶
A=zeros(M1,M1);%对应每列保存1:p(AIC准则下的阶数)对应的参数fai
for p=1:M1;
fai=zeros(p,N);
v1=randn(1,p)*0.001;
v2=randn(1,p)*0.001;
I=eye(p,p);
k=repmat(I,[1,1,N+1]);
%k(t,t-1)=k(:,:,t);
g=zeros(p,N);
k_k=repmat(I,[1,1,N+1]);
fai(:,p)=zeros(p,1);
k(:,:,p+1)=c*I;
xi=zeros(p,N+M);
arfa=zeros(N,1);
for i=p:N;
xi(:,i)=xg(i:-1:i-p+1);
end
%%%预先假设v1、v2的方差,v2方差为(0.001~0.01)*(x(t)的方差)
%Jmin=0.001*var(z(p+1:N)); %这里消除认为设定
for i=p:N-1;
%g(:,i+1)=k(:,:,i+1)*xi(:,i)*inv(xi(:,i)'*k(:,:,i+1)*xi(:,i)+Jmin);
g(:,i+1)=k(:,:,i+1)*xi(:,i)*inv(xi(:,i)'*k(:,:,i+1)*xi(:,i));
arfa(i+1)=xg(i+1)-xi(:,i)'*fai(:,i);
fai(:,i+1)=fai(:,i)+g(:,i+1)*arfa(i+1);
k_k(:,:,i+1)=k(:,:,i+1)-g(:,i+1)*xi(:,i)'*k(:,:,i+1);
k(:,:,i+2)=k_k(:,:,i+1)+10^-2*I;
end
fai(:,N);
for i=1:p;
A(i,p)=fai(i,N);
end
x_yb=zeros(N-p,1);
x_y=zeros(N-p,1);
for i=p+1:N;
x_y(i)=fai(:,N)'*xi(:,i-1);
end
for i=p+1:N;
x_yb(i)=x_y(i)+avrg;
end
g=zeros(N-p,1);
for i=p+1:N;
g(i)=z(i)-x_yb(i);
end
h=sum(g.^2);
%a(p)=log(h/(N-2*p-1))+2*p/N;
a(p)=log(h/(N-2*p-1))+2*(p+1)*log(N)/N;
end
[C,p]=min(a)
Fai=zeros(p,1);
for i=1:p;
Fai(i)=A(i,p);
end
%预报
x_y1=zeros(N+M,1);
xg(N+1)=Fai'*xg(N:-1:N-p+1);
for i=N+1:N+M;
xg(i+1)=Fai'*xg(i:-1:i-p+1);
end
for i=1:N+M;
x_y1(i)=xg(i)+avrg;
end
%绘图
figure(1);
k=N-10:N+M;
k1=N-10:N+18;
plot(k1,z(k1),'r*-');
hold on;
k=N:N+M;
plot(k,x_y1(k),'b*-');
xlabel('t步长');
ylabel('角');
legend('GBPUSD磅美真实值','预报值');
%性能指标
%y=max(abs(x(N+1:N+n)));
%z=zeros(n,1);
%for i=N+1:N+n
%z(i)=x(i)-x_y1(i);
%end
%S=sum(z.^2);
%S1=sqrt(S/n);
%yit=S1/y