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);