kalman_loglik_initx_free.m

function [negloglik, xtt, Vtt] = kalman_loglik_initx_free(v,y,V1, demoit) 
%demoit is a flag to plot steps
%v is vector of parameters B Q R (mu s2 s2np); y is data; initx is just that; V1 is var of initx

%STEP 1 set the initial estimates
%Initialize by getting D-H estimates
T=length(y);
B = v(1); %mu
Q = exp(v(2)); %s2
R = exp(v(3)); %s2np
initx = v(4);

%STEP 2 using these initial estimates generate an estimate of x(t)
%using a forward pass of the Kalman filter.  This gives
%you the ML estimate of x(t) given y(1:T) and the parameter estimates.

%initialize
xtt=zeros(1,T);  Vtt=zeros(1,T); 
xtt1=zeros(1,T); Vtt1=zeros(1,T);    
Ft=zeros(1,T); vt=zeros(1,T);
At = zeros(1,T);

%set up for missing values
At(y ~= -99) = 1;
At(y == -99) = 0;
y(y == -99) = 0;

%forward pass gets you E[x(t) given y(1:t)]
x10=initx;  
V10=V1; 
for(t=1:T)
   if(t==1)
   	xtt1(1) = initx; %denotes x_1^0
      Vtt1(1) = V1; %denotes V_1^0
   else
      xtt1(t) = xtt(t-1) + B; %xtt1 denotes x_t^(t-1); Harvey 3.2.2a
      Vtt1(t) = Vtt(t-1) + Q; %Harvey 3.2.2b
   end
   Kt = Vtt1(t)*At(t)/(At(t)*Vtt1(t)*At(t)+R);
   Ft(t) = At(t)*Vtt1(t)*At(t)+R;
   vt(t) = y(t)-At(t)*xtt1(t);
   xtt(t) = xtt1(t) + Kt*( y(t) - At(t)*xtt1(t) ); %Harvey 3.2.3a
   if(demoit==1 & t>1) plot([t-1 t], [xtt(t-1) xtt(t)], '-r'); pause; end %show step if demoit is true
   Vtt(t) = Vtt1(t)*(1 - Kt*At(t)); %Harvey 3.3.3b
end
KT = Kt;

t=1:T;
if(demoit == 2) 
   plot(t, xtt, '-b');
   hold on;
   plot(t(y ~= 0), y(y ~= 0), '.r');
   hold off;
end %show whole line

%Calculate negative log likelihood
negloglik =  (1/2)*sum(vt.^2./Ft) + (1/2)*sum(At.*log(abs(Ft))) + (T/2)*log(2*pi);