/* DYMIMIC, a GAUSS language program to estimate a Watson-Engle Kalman  **
** filter DYMIMIC model.                                                **
**************************************************************************
** The program has three main parts.                                    **
** (1) A keword procedure "dymimic" sets up a command line of the form  **
**     dymimic <list of ind. vars> on <list of dep. vars> from <file>   **
**     Entering "dymimic" alone will produce a syntax prompt            **
** (2) The kf procedure sets up the call to MAXLIK.  It reads the data, **
**     and constructs starting values from a series of regressions.  It **
**     uses the first dependent variable as a proxy for the latent      **
**     state variable to produce starting values that are close to      **
**     ML estimates if this variable is a good proxy.                   **
** (3) The li procedure calculates log likelihood by implementing the   **
**     KF updating algorithm.  It is based on Beck's (1990) published   **
**     likelihood procedure.                                            **
**************************************************************************
** The first executable line: use \gauss\maxli0; may need to be changed **
** to name the version of maxlik that you own.                          **
**************************************************************************
** Priors: PLINIT=1 and XLINIT=mean of Y1, an assumption that is        **
** appropriate for stationary Y series, but not for unit root series.   **
** It amounts to the assumption that x(0)=mean of x(t).  The latent     **
** dependent variable x is expressed in the metric of Y1, the first     **
** dependent variable, a result that is produced by constraining        **
** ALPHA(1)=1 and BETA(1)=0.                       J. Stimson 11/19/96  **
**************************************************************************/
use \gauss\maxli0;
closeall;
f1=0; @initialize@
let X=0;
xlinit=0;
plinit=1;
B=0;S=0;
let dvindex=0;
let ivindex=0;
nind=0;ndep=0;nparm=0;dataset="";vars="";

keyword dymimic(s);   @ proc translates command line input into variable lists@
 local dep,ind,tok,s;
 dep = 0;ind = 0;tok = "";
 do until (upper(tok) $== "ON");
    { tok,s } = token(s); @token returns the 1st item from string "s"@
    if s $== "";goto usage;endif;
    ind = ind|tok;    @vector of indepdent variable names@
    endo;
 ind = trimr(ind,1,1); @remove word ON from list@

 do until (upper(tok) $== "FROM");
    { tok,s } = token(s);
    if s $== "";goto usage;endif;
    dep = dep|tok;    @vector of dependent variable names@
    endo;
 dep = trimr(dep,1,1);
 { tok,s } = token(s);
 dataset=tok;
 cls;
 call kf(dataset,ind,dep);
 retp;
Usage: @user documentation produced whenever command line parsing fails@
 print;
 print "***************************************************************";
 print "*                    SYNTAX FOR DYMIMIC                       *";
 print "***************************************************************";
 print "*                                                             *";
 print "*  \175 DYMIMIC <indvar list> ON <depvar list> FROM <dataset>    *";
 print "*                                                             *";
 print "*  Example:                                                   *";
 print "*                                                             *";
 print "*  \175 DYMIMIC x1 x2 on y1 y2 y3 from datafile                  *";
 print "*                                                             *";
 print "***************************************************************";
endp;

/* begin li proc *************************************************** */
proc li(B,X);
  local ptgt,xtgt,i,xtgl,ptgl,INN,H,R,ALPHA,BETA,HI,q,GMMA,phi,LIKEI,intcept,
    YT,ZT,d1,d2;
  YT=X[.,nind+1:nind+ndep];  @note that X passed by MAXLIK now has vars in@
  ZT=X[.,1:nind];            @order determined by vars vector@
  S=zeros(rows(X),1);
  i=1;
  xtgt=xlinit;
  ptgt=plinit;
  d1=1;   @d1 and d2 are beginning and ending rows in B vector@
  intcept = B[d1];  @ set d1=0 to suppress intercept @
  d1=d1+1;
  phi = B[d1];
  d1=d1+1;d2=d1+nind-1;
  GMMA = B[d1:d2];
  d1=d2+1;d2=d1+ndep-2;
  BETA = 0|B[d1:d2];
  d1=d2+1;d2=d1+ndep-2;
  ALPHA = 1|B[d1:d2];
  d1=d2+1;d2=d1+ndep-1;
  R = diagrv(zeros(ndep,ndep),B[d1:d2]^2);
  d1=d2+1;
  q = B[d1]^2;
  LIKEI = zeros(rows(X),1);
  do while i <= rows(X);
        xtgl = intcept+phi*xtgt+ZT[i,.]*GMMA;
        ptgl = ptgt*phi^2+q;
        INN = YT[i,.]'-ALPHA*xtgl-BETA;
        H = ALPHA*ptgl*ALPHA'+R;
        HI = invpd(H);
        xtgt = xtgl+ptgl*ALPHA'*HI*INN;
        S[i] = xtgt;
        ptgt = ptgl-ptgl*ALPHA'*HI*ALPHA*ptgl;
        plinit = ((i-1)/i)*plinit+(1/i)*ptgt;
        LIKEI[i,1] = -.5*(ln(det(H))+INN'*HI*INN); @log likelihood(i)@
        i=i+1;
        endo;
  retp(LIKEI);
endp;

/* begin kf proc *************************************************** */
proc (0) = kf(dataset,ind,dep);
  local l,d2,plinit,t,ll,g,vc,ret,n,rvar;
  local Y,dv,iv,XCOL,X2,BY,ap,bp,SV,yss,lastr,name;
  external proc indices;
  vars=ind|dep;
  nind=rows(ind); @# of ind vars@
  ndep=rows(dep); @# of dep vars@
  {ind,ivindex} = indices(dataset,ind); @gets column numbers for indexing@
  {dep,dvindex} = indices(dataset,dep);
  nparm=2+nind+2*(ndep-1)+ndep+1; @total number of parameters to be estimated@
  B=zeros(nparm,1);
  open f1=^dataset;
  n=rowsf(f1);  @N is number of rows in file@
  X=readr(f1,n);
  __title="ML Estimation of Kalman Filter DYMIMIC Specification";
  _mlalgr=2;
  _mlparnm="G0"|"Phi";
  l=1;name=vars; @build a vector of parameter names, first gammas"@
  do while l<=nind;   @assign nind gamma names@
     name[l]="G" $+ftocv(l,1,0);  @ftocv translates the number l into a string@
     _mlparnm=_mlparnm|name[l]; @add on to existing name vector@
     l=l+1;
     endo;
  l=2;
  do while l<=ndep;     @assign ndep-1 beta names@
     name[l]="B" $+ftocv(l,1,0);
     _mlparnm=_mlparnm|name[l];
     l=l+1;
     endo;
  l=2;
  do while l<=ndep;    @assign ndep-1 alpha names@
     name[l]="A" $+ftocv(l,1,0);
     _mlparnm=_mlparnm|name[l];
     l=l+1;
     endo;
  l=1;
  do while l<=ndep;  @assign ndep R names@
     name[l]="R" $+ftocv(l,1,0);
     _mlparnm=_mlparnm|name[l];
     l=l+1;
     endo;
  _mlparnm=_mlparnm|"Q";

/* set up regressions: Y(i) = beta + alphaY(1) to determine alpha,beta
   for dependent variables 2 through p */
  XCOL=ones(n,1)~X[.,dvindex[1]]; @ begin est of beta and alpha start values@
  dv=2;
  do while dv<=ndep;
     bp=2+nind+dv-1;          @ beta pointer @
     ap=2+nind+dv-1+(ndep-1); @ alpha pointer @
     Y=X[.,dvindex[dv]];
     {BY,rvar}=reg(XCOL,Y);  @reg is a mini regression proc that returns@
     B[bp]=BY[1];            @B and residual variance@
     B[ap]=BY[2];
     dv=dv+1;
     endo;

/* form a sum (SV) across all dependent indicators, and then divide by ndep*/
  dv=2;
  SV=X[.,dvindex[1]]; @ 1st dv begins state var est. @
  do while dv<=ndep;
     bp=2+nind+dv-1;          @ beta pointer @
     ap=2+nind+dv-1+(ndep-1); @ alpha pointer @
     SV=SV+(X[.,dvindex[dv]]-B[bp])/B[ap]; @add each dvar modified by parms@
     dv=dv+1;
     endo;
  SV=SV/ndep;  @now divide by ndep to get average@

/* SV is now an average of the dep vars, expressed in the metric of Y1,
   now use it to estimate transition model parameters phi and gamma */
  Y=SV;  @ begin estimation of con, phi, gamma start values @
  X2=ones(n,1)~lag1(Y);
  iv=1;
  do while iv<=nind;
     X2=X2~X[.,ivindex[iv]]; @X2 is | 1.0 y(t-1) x1(t) x2(t) ...|@
     iv=iv+1;
     endo;
  Y=trimr(Y,1,0);   @ both X2 and Y need to have 1st row removed from lagging@
  X2=trimr(X2,1,0);
  {BY,rvar}=reg(X2,Y);
  B[1:nind+2]=BY;    @now assign estimates to start value vector @
  clear X2;

/* do regressions: Y(i) = B(0) + phi*Y(i,t-1) + XB to estimate residual
   variance of each Y, "R" */
  X2=ones(n,1)~lag1(SV);  @ first, build X2 from trans model @
  iv=1;
  do while iv<=nind;
     X2=X2~X[.,ivindex[iv]];
     iv=iv+1;
     endo;
  X2=trimr(X2,1,0);  @ trim 1st row for lagging @
  l=ap+1;dv=0;      @ points to R1 @
  lastr=ap+ndep;
  do while l<=lastr;
     dv=dv+1;
     Y=X[.,dvindex[dv]];  @ once thru for each Y indicator @
     Y=trimr(Y,1,0);     @ trim to match X2 @
     {BY,rvar}=reg(X2,Y);
     B[l]=sqrt(rvar);    @now assign R estimates to start value vector @
     l=l+1;
     endo;

   @ now resid variance on transition model @
  Y=trimr(SV,1,0);
  {BY,rvar}=reg(X2,Y);
  B[nparm]=sqrt(rvar);    @this is Q est.@

  plinit = 1;
  xlinit = meanc(X[.,dvindex[1]]); @x(0) is mean of Y1@

  t=dataset$+".OUT";  @print out starting values for inspection@
  output file=^t on;
  "Parm  Starting Values Variable";
  l=1;do while l<=nparm;
  format /rd 4,0;
  print l $_mlparnm[l];;
  format /rd 8,4;
  print B[l];;
  if l<3;print "tr model";endif;
  if l>2 and l<=(nind+2);print $ind[l-2];endif;
  if l>(nind+2) and l<=(nind+ndep+1);print $dep[l-(nind+2)+1];endif;
  if l>(nind+ndep+1) and l<=(nind+2*ndep);print $dep[l-(nind+ndep+1)+1];endif;
  if l>(nind+2*ndep) and l<nparm;print $dep[l-(nind+2*ndep)];endif;
  if l==nparm;print "tr model";;endif;
  l=l+1;endo;
  "     Any key to continue";;wait;print;print;

  clear Y,X2,BY,SV;
  {B,ll,g,vc,ret}=maxprt(maxlik(dataset,vars,&li,B));  @heavy lifting@
  l=1;print;"Correlations and Communalities: State Variable and Y Indicators:";
  local cor,smean,sss,rss;
  smean=meanc(S);
  S=S-smean;
  sss=S'S;
  do while l<=ndep;
     X2=ones(n,1)~X[.,dvindex[l]];
     {BY,rvar}=reg(X2,S);
     rss=rvar*(n-2);
     cor=sqrt(1-(rss/sss));
     print $dep[l];;cor cor^2;
     l=l+1;
     endo;
  output off;
endp;

proc (2) = reg(x,y);  @ mini regression procedure @
  local b,res,rvar;
  b=invpd(x'x)*x'y;
  res=y-x*b;
  rvar=res'res/(rows(x)-cols(x));
  retp(b,rvar);
endp;
end;