@------------------------------------------------------@
@ Kalman filter program, sampling error in pid data    @
@ w/o intercept using mean-centered data               @
@ *** program allows for unevenly spaced observations  @
@------------------------------------------------------@

library maxlik;
#include maxlik.ext;
maxset;
open f1=pidsim.dat;
x=readr(f1,500);
dataset="pidsim";
let dep= obsdev;
__title="Dynamic Shock Errors";
_max_ParNames="PHI"|"Process"|"Measrmnt";
_max_Algorithm=2;
_max_GradTol=1e-3;
_max_MaxIters=150;
_max_LineSearch=3;
_max_FinalHess;
_max_CovPar=2;
vars=dep;
B ={.9,1.1,5};  @ starting values @
yt=x[.,4];
days=x[.,5];
plinit = 1 ;
xlinit = 0 ;
proc li(b,x);

local ptgt,xtgt,i,xtgl,ptgl,inn,h,r,hi,q,
      gmma,phi,phip,likei,elapse,prevday,qp,qinit,qsum,j,v,e1;
i=1;
xtgt=xlinit;
ptgt=plinit;
qinit=0;
prevday=1;
r = b[3,1];
phip = b[1,1];
qp = b[2,1];
qsum=0;j=0;
likei = zeros(rows(x),1);
format /rd 8,4;
do while i <= rows(x);
        elapse = days[i,1]-prevday;
        phi  = phip^elapse;

   do while j < elapse;
   qsum=qsum+phip^(2*j);
   j=j+1;
   endo;

        q    = qsum*qp;

        xtgl = phi*xtgt;
        ptgl = ptgt*phi^2+q;
        inn = yt[i,.]'-xtgl;
        h = ptgl+r;
        hi = invpd(h);
        xtgt = xtgl+ptgl*hi*inn;
        ptgt = ptgl-ptgl*hi*ptgl;
@print days[i,1]~elapse~phi~q~qsum~j;@ @ uncomment to chk updating of q @
        q=qinit;j=qinit;qsum=qinit;
        prevday=days[i,1];
        likei[i,1] = -.5*(ln(det(h))+inn'*hi*inn);
        i=i+1;
        endo;
retp(likei);
endp;

   {b,ll,g,vc,ret}=maxprt(maxlik(dataset,vars,&li,b));

print ret;
f1=close(f1);