clear
t=500;
Q=0.001;
R=0.01;
F=1;
mu=0;
e1=randn(t,1)*sqrt(R);
e2=randn(t,1)*sqrt(Q);
Beta=zeros(t,1);
Y=zeros(t,1);
X=randn(t,1);
for j=2:t
Beta(j,:)=Beta(j-1,:)+e2(j,:);
Y(j)=X(j,:)*Beta(j,:)'+e1(j);
end
beta0=zeros(1,1);
p00=1;
beta_tt=[];
ptt=zeros(t,1,1);
beta11=beta0;
p11=p00;
for i=1:t
x=X(i);
beta10=mu+beta11*F';
p10=F*p11*F'+Q;
yhat=(x*(beta10)')';
eta=Y(i,:)-yhat;
feta=(x*p10*x')+R;
K=(p10*x')*inv(feta);
beta11=(beta10'+K*eta')';
p11=p10-K*(x*p10);
ptt(i,:,:)=p11;
beta_tt=[beta_tt;beta11];
end
beta2 = zeros(t,1);
wa=randn(t,1);
i=t;
p00=squeeze(ptt(i,:,:));
beta2(i,:)=beta_tt(i:i,:)+(wa(i:i,:)*chol(p00));
for i=t-1:-1:1
pt=squeeze(ptt(i,:,:));
bm=beta_tt(i:i,:)+(pt*F'*inv(F*pt*F'+Q)*(beta2(i+1:i+1,:)-mu-beta_tt(i,:)*F')')';
pm=pt-pt*F'*inv(F*pt*F'+Q)*F*pt;
beta2(i:i,:)=bm+(wa(i:i,:)*chol(pm));
end
plot([beta_tt beta2 Beta])
axis tight
legend('Kalman filter estimated \beta_{t}','Draw from H(\beta_{t})','true \beta_{t}');