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

%%%%%%%%%基于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