/****************************************************************************
              FILE: KALMAN.OPT
              Programs for maximizing Gaussian likelihood functions
              NOTE: IF the program is run with non-stable date or with 
              poor starting values it will crash. If the problem is
              the initial values then conditioning on the initial values
              (rather than using the unconditional distribution for the 
              first Pt will usually take care of the problem (set Pinit=0
              before running the kalman filter). 
              When you get convergence with the conditional likelihood
              then maybe use the estimated
              parametervalues as starting values for a second run with
              no conditioning (Pinit=1).  
***************************************************************************/
        library optmum;
/***************************************************************************
                 CREATE GLOBALS
****************************************************************************/
clear pinit ;
       load gnp5991 ;         @DATASET@
       datatemp = ln(gnp5991) ;
       nobstemp = rows(datatemp) ;
       dlgnp82 = datatemp[2:nobstemp]-datatemp[1:nobstemp-1] ;
       data = dlgnp82 ;
       Nobs= rows(data) ;
       clear datatemp,nobstemp,dlgnp82 ;                        

       "no of obs: ";; Nobs ;
       meandata = sumc(data)/Nobs ;
       data = data - meandata ;
/*       _opstmth = "bfgs stepbt";       /* starting algorithm */
@       _opstmth = "steep one";       /* starting algorithm */  @
       _opmdmth = "bfgs stepbt";     /* take-over algorithm */  */
       _opditer = 10 ;
       _opstep=2;   /* Optmum uses STEPBT procedure                        */
       _opalgr=2;   /* Optmum uses BFGS algorithm          */
 /*      _opalgr=1,2,3,4,5   Optmum uses Steep,bfgs,Sbfgs,DFP,Newton    */

/**************************************************************************/
@CHOOSE DIMENSIONS OF ARMA MODEL:@
kdim = 1 ;  ldim = 1 ;
m = maxc(kdim|ldim+1) ;
if kdim > 0 ;
  a=zeros(kdim,1) ;
endif ;
if ldim > 0 ;
  b=zeros(ldim,1) ;
endif ;

/****************** USE VAR instead ***************
@starting values@
a[1,1] = .5 ;
a[2,1] = .1 ;
b[1,1] = .2 ;
sigma  = 25 ;                        *************/


proc (2) = VECVAR(nlags) ;
local Z,Y,UHAT,sigmahat,jj,B,p ;
p = cols(data) ;
@p is larger than one in the case of VECTOR autoregression@ ;
Y = data[nlags+1:Nobs,.]' ;
Z = zeros(nlags,Nobs-nlags) ;
jj = 1;
do while (jj < nlags+1)  ;
jj ;
Z[nlags-jj+1,.] = data[jj:(Nobs-nlags-1+jj),.]' ;
jj = jj + 1 ;
endo ;
B =  Y*Z'inv(Z*Z') ;
UHAT = Y - B*Z ;
sigmahat = 1 ;
retp(B,sigmahat ) ;
endp ;

proc kalman(psi) ;
local KK,Z,R,Q,Pt,Ptpred,Ft,FTinv,alpha,alphapr,nu,LL,t,ypred ;
KK = zeros(m,m) ;
if m > 1 ;
   Z = 1~zeros(1,m-1) ;
else;
   Z =1 ;
endif ;
R = zeros(m,1) ;
Pt = zeros(m,m) ;

@transfer parameters to matrices@
if kdim > 0 ;
  KK[1:kdim,1] = psi[1:kdim] ;
endif ;
if m > 1 ;
  KK[1:m-1,2:m] = eye(m-1) ;
endif ;
if ldim > 0 ;
  R[1:1+ldim] = 1|psi[kdim+1:kdim+ldim] ;
else ;
  R[1,1] = 1 ;
endif ;
Q = psi[kdim+ldim+1]^2 ;

@initialize the Kalman Filter@
if Pinit == 1 ;
  Pt = reshape(inv(eye(m*m)-KK.*.KK)*vec(R*Q*R'),m,m)  ;
else ;
  Pt = 0 ;
endif ;
@"pt: ";; pt ;@
alpha = zeros(m,1) ;
LL =0 ;
t = 0 ;

/********the main loop **************************************/
do until t==(Nobs-1) ;
  t= t+1 ;

  @updating equations@
  alphapr = KK*alpha ;
  ypred = Z*alphapr ;
  Ptpred = KK*Pt*KK' + R*Q*R' ;
  @find the conditional likelihood for x_t@
  nu = data[t] - ypred ;
  Ft = Z*Ptpred*Z' ;    @note - H is zero here@
  Ftinv = inv(Ft) ;
@  "Ft: " ; Ft ; @
  LL = LL - .5*ln(Ft) - .5*nu'Ftinv*nu ;

  @updating equations@
  alpha = alphapr + Ptpred*Z'*Ftinv*nu ;
  Pt = Ptpred - Ptpred*Z'*Ftinv*Z*Ptpred' ;
endo ;

retp(-LL) ;
endp;

/**********************************************************************/
/************ EXECUTION ***********************************************/
/**********************************************************************/

output file = ec266a.out reset;
"----------------------------------------------------------------------";
"                     ec266a.OUT                                       ";
"----------------------------------------------------------------------";
"Date: " datestr(date) ;
"Time: " timestr(time) ;
"----------------------------------------------------------------------";

if kdim>0 ;
  "*** VAR ***" ;
  {ahat,sigma2} = vecvar(kdim) ;
  "ahat: ";; ahat ;
  sigma = sqrt(sigma2) ;
  "sigma: " ;; sigma ;
  pause(2) ;
  " " ;
else; 
  b = .1 ; sigma = 1 ;
endif ;


" ************* ML **************** " ;
" " ;
if Pinit == 1 ;
   "Initial conditions: Stationary distribution. " ;
else ;
   "Initial conditions: initial value given. " ;
endif ;
"" ;
if (kdim > 0) AND (ldim > 0) ;
  start = ahat'|b|sigma ;
elseif kdim > 0 ;
  start = ahat'|sigma ;
elseif ldim > 0 ;
  start = b|sigma ;
endif ;

  Pinit = 0 ; @0 for first Pt = 0 , 1 for stationary distribution@
   {parm,vof,g,ret}=optmum(&kalman,start) ;
if ret ==0 ;
  "NORMAL CONVERGENCE" ;
else ;
  "CONVERGENCE PROBLEM: ERROR CODE: ";; ret ;
endif ;

  Pinit = 1 ; @0 for first Pt = 0 , 1 for stationary distribution@
   {parm,vof,g,ret}=optmum(&kalman,parm) ;
if ret ==0 ;
  "NORMAL CONVERGENCE" ;
else ;
  "CONVERGENCE PROBLEM: ERROR CODE: ";; ret ;
endif ;

/************** standard deviations ************************************/
stda = -100000 ; stdb = -100000 ;
/* Put in the right formula's here */

format /mb1 /ros 16,8 ;

if kdim > 0 ;
"ahat:     ";; parm[1:kdim]'  ;
"std dev.: ";; stda ;
endif ;
if ldim > 0 ;
"bhat:     ";; parm[kdim+1:kdim+ldim]'  ;
"std dev.: ";; stdb ;

endif ;
"sigmahat: ";; parm[kdim+ldim+1]  ;

 output off;
"----------------------------------------------------------------------";
  end;