 % the following code implements the scalar
 % state space model:
 % x_t+1=x_t +r_t
 % y_t  =x_t + v_t
 %
 % E(x_0)=0, sqrt{var(x_0)}=sigma
 % E(v_t)=0, sqrt{var(v_t)}=sigma_v
 % E(r_t)=0, sqrt{var(r_t)}=sigma_r
 % x_, {v_t} and {r_t} are all independent
 clear all
 sigma=1;
 sigma_v=.05;
 sigma_r=1;
 
 T=60  % define the number of time steps
 
 x=randn*sigma; % generating the initial state
 K=[]; % to store the gain
 xHat=[0];  % to store the estimates, initial estimate is 0
 y=[]; % to store the sequence of observations
 
 
 for i=1:T
 temp=sigma(i)^2+sigma_r^2;
 temp=temp/(temp + sigma_v^2);
 K=[K,temp];  % appending new gain    
     
 % generating new state and observation
 x=[x, x(i)+randn*sigma_r];
 y=[y,x(i+1)+randn*sigma_v];
 % calculating xHat
 xhat=K(i)*(y(i)-xHat(i))+xHat(i);
 xHat=[xHat,xhat ];    
     
     
 % recursion for the error variance
 temp=sigma(i)^2+sigma_r^2;
 temp=sigma_v^2*temp/(sigma_v^2 + temp);
 sigma=[sigma,sqrt(temp)];
 
 end
 
 
 figure
plot(x(2:end))
hold on
plot(y,'xr')
plot(xHat(2:end),'k')
title('true state =blue, obs =cross, xHat=black')