/****** Multivariate Structural Time Series Gauss Code   **********

Author   : Geoffrey Shuetrim
 London School of Economics
 & Reserve Bank of Australia
 Email: G.Shuetrim@lse.ac.uk
Date: 6 March 1996

This code is written and submitted for public, non-commercial use.
There are no performance guarantees.
Please acknowledge this code (and its author) if you find it useful.


This is an upgrade of the univariate Structural Time Series Gauss
Code (previously listed the archive) to handle a very wide class of
multivariate models.
********************************************************************/


/*****************************************************************

SSModel : This program does full maximum likelihood estimation of
Structural Times Series models.  Thus, it is useful for estimation
such models as Random walks with noise, ARIMA specifications,
seasonality and cycles and time varying parameter models.
It can also be used to estimate cointegrating systems.
Essentially, it is capable of estimating models that can be
expressed in the state space form, so long as the the errors are
Gaussian (Because the Kalman filter is used to extract the
prediction errors).
For more information on such models a very useful reference
is:

Harvey, A.C. (1989), "Forecasting Structural Times Series Models
and the Kalman Filter", Cambridge University Press.
This reference will be shortened to FSK.

This code has been developed specifically for use in estimating time
varying parameter models as a part of my PHd research.
While I am fairly confident that it is working well, bugs may remain.
If people find errors/innefficiencies or can see improvements in the
algorithms used, I would appreciate hearing of them.  Also, if anyone
can figure out why the likelihood function is misspecified severely
when I concentrate out the measurement equation error variance, I
would love to hear the explanation.  It cost me about a month of messing
around before I gave up on the concentrating out approach.

**********************************************************************/

/*********************************************************************
The source code has been designed to separate the code that users need to
write for their particular specification and the code that is generally
applicable to the univariate state space form.  The generally applicable
code should never need to be altered.

The general State Space Form (SSF) can be written as:

Measurement equation:

yt = Zt * at + dt + et

  yt is N by 1
  Zt is N by K
  at is K by 1
  dt is N by 1
  et is N by 1

  E(et) = 0
  Var(et) = h which is N by N

Transition Equation:

at = T * at-1 + ct + Rt * ut

  T is K by K
  ct is K by 1
  Rt is K by L
  ut is L by 1 where L <= K.

  E(ut) = 0
  Var(ut) = Q which is L by L.

         *******************##########*********************

The generally applicable source code (that users never alter )
consists of the following procedures:

1. Kalman(Parameters)

This procedure filters the data, extracting information about the
state in
each of the time periods (its algorithm is described in FSK p. 106).
Note that the code is capable of detecting missing observations in the
measurement variable.  It checks yt at each period to find out whether the
data element is missing values (uses ismiss(yt)) and skips the updating
equation if it is missing.

2. LL(Parameters)

This procedure constructs the Log Likelihood function for the full set of
'T' observations where the data runs from 1 to T.  Although it is possible
to do concentrated ML, the LL procedure specified here does not concentrate
out the variance of the measurement equation error.  See FSK p. 126.

3. Smoother(Parameters)

The smoother generates estimates of the state in each time period, t, using
information available from the entire sample (1 to T) rather than from the
observations up to and including t.  It implements the fixed interval
smoother described in FSK on p. 154.


These procedures are attached after the example code for estimating
a random walk with noise model.

****************#############********************

GLOBAL CONSTANTS

_Inf.

This constant gives the value to be used in construction a diffuse prior
for the initial state vector a0.  It is used in the procedure that
constructs a0 and more specifically, its covariance matrix P0.

 ******************#######********************

Now for the part that requires more participation from users.  Here is the
specification of the procedures that users are required to write themselves.
An example is included of these procedures for the Random Walk with
Noise model.

The procedures that have to be written are:

1. Initial(Parameters)

This procedure specifies the state vector and its covariance matrix at
t0.

4. SystemM(Parameters)

This procedure sets up all of the system and data matrices:
T,H,Q,R,d,c,y,Z
Note that the matrices are set up for each observation allowing them
to be time varying.

**********************************************************************/

/*********************************************************************

Following is the mainline of code and the necessary procedures to
estimate the a bivariate random walk with noise model.

Note that this system can be easily modified to capture a
simple cointegrating system.
**********************************************************************
******/

@ Open the Required Library Files @

    library cml,pgraph;

    cmlset;
    graphset;


@ Set the Global Parameters @

tempry = 0;

    _Inf = 1e5;

    obs = 100; @ used here to generate the data @

@ Simulation Parameters @

    sigmaM1 = 0.1; @ variance of the measurement equation @
    sigmaM2 = 0.1;

    sigmaU1 = 0.01; @ Variance of the stochastic level @
    sigmaU2 = 0.01;

    parms = sigmaM1|
        sigmaM2|
        sigmaU1|
        sigmaU2;

@ Generate the Data @

    @ Generate the Measurement equation Residuals @
      M1 = rndn(obs,1) * sqrt(sigmaM2);
      M2 = rndn(obs,1) * sqrt(sigmaM2);

    @ Generate the Level equation Residuals @
      U1 = rndn(obs,1) * sqrt(sigmaU2);
      U2 = rndn(obs,1) * sqrt(sigmaU2);

    u = 0~0;
    yvec = 0~0;
    uvec = 0~0;
    t = 1;
    do while t <= obs;
        u = u + (U1[t]~U2[t]);
        y = u + (M1[t]~M2[t]);
        yvec = yvec|y;
        uvec = uvec|u;
        t = t + 1;
    endo;

    yvec = trimr(yvec,1,0);
    uvec = trimr(uvec,1,0);

    data = yvec;

@ Define restrictions on the parameter space, Variances are positive @

   _cml_C = eye(4);
   _cml_D = zeros(4,1);

@ Do the Estimation and output results @

    {parms,f,g,h,retcode} = CML(data,0,&LL,parms);
    {parms,f,g,h,retcode} = CMLCLPrt(parms,f,g,h,retcode);

@ Generate the Filtered Data using the Smoother and the natural parm @

    {BKvec,VVec,FVec,PPVec,PuVec} = Kalman(parms,data);
    {BSvec,PSVec} = Smoother(parms,data);

@ Graph the smooth, filtered and actual state (2 dimensional) through
time @

    xy(seqa(1,1,rows(UVec)),UVec~BKVec~BSVec);

end;

/* NOW FOR THE PROCEDURES THAT USERS WILL HAVE TO CUSTOMISE FOR THEIR
MODEL*/

/*********************************************************************
INITIAL : Initialises the state vector and its covariance matrix

INPUT:
1.   Parms

OUTPUT:
1.   Initial State Vector at t0
2.   Initial State Covariance Matrix at t0
**********************************************************************/
proc(2) = Initial(parms);
local   Bu,Pu;
    Bu = zeros(2,1);
    Pu = eye(2) .* _Inf;
    retp(Bu,Pu);
endp;

/*********************************************************************
 SystemM : Generate the System Matrices

INPUT:
1.  Parameter Vector (K by 1)

OUTPUT:
1.  System Matrices for the Kalman Filter as described by Harvey
1.  yt and Zt
**********************************************************************/
proc(8) = SystemM(parms,Data);
local T,H,Q,R,d,c,y,Z;

@ Transition Matrix @
    T = eye(2);
@ Transition equation Covariance Matrix @
    Q = diagrv(eye(2),parms[3 4]);
@ Measurement Equation Covariance Matrix @
    H = diagrv(eye(2),parms[1 2]);
@ Transition Equation Error Coefficient Matrix @
    R = eye(2);
@ Transition Equation Interventions @
    c = zeros(2,1);
@ Measurment Equation Observations @
    y = Data';
@ Measurment Equation Matrix@
    Z = eye(2);
@ Measurement Equation Interventions @
    d = zeros(2,1);

    retp(T,h,Q,R,d,c,y,Z);
endp;

/* NOW FOR THE PROCEDURES THAT NEVER HAVE TO BE ALTERED */

/*********************************************************************
*******
KALMAN : Generates the states vector using the Kalman Filter
INPUT:
1.   Parameters
OUTPUT:
1.   Bvec is the State Vector through time.  It is a T by 1 vector.
2.   Vvec is the T by 1 vector of prediction errors.
3.   Fvec is the T by 1 vector of prediction error variances.
4.   PpVec vectorised Pp matrices for each observation 1 through to T.
5.   PuVec vectorised Pp matrices for each observation 1 through to T.
**********************************************************************
******/

proc(5) = Kalman(parms,data);
local   Bp,         @ State Vector Prediction @
        Bu,         @ Updated State Vector @
        Bvec,       @ stores the state vector for all time (T by K)
matrix @
        Pp,         @ state covariance prediction @
        PpVec,      @ Store vectorised Pp matrices  for all time @
        Pu,         @ Updated state covariance estimate @
        PuVec,      @ Store Pu matrices @
        yp,         @ predicted measurement variable, Firm Return @
        yu,         @ Updated measurement variable estimate, Firm
Return @
        y,@ Actual Dependent Variable @
        Z,@ Actual Explanatory Variables @
        V,@ Prediction error @
        Vvec,       @ Prediction Error vector @
        F,@ Prediction error variance @
        Prefix,     @ Kind of Weighting term in the updating equation
@
        Fvec,       @ Prediction Error Variance matrix F @
        DFvec,      @ Determinant of F @
        T,Q,H,R,d,c,@ System Matrices @
        W,@ the matrix allowing for missing observations @
        mask,       @ The mask setting 1 for observations that occur @
        vmask,Fmask,@ Templates for storing the v and F matrix output
@
        i;@ Counter over time periods in the data @

@ Get the i=0 State vector and Covariance Matrix @
    {Bu,Pu} = Initial(parms);

@ Filter i=1..T, using index i cos T already used for a system matrix
@

    i = 1;
    do while i <= rows(Data);

        @ Get the System Matrices and Data @
        {T,H,Q,R,d,c,y,Z} = SystemM(parms,Data[i,.]);

        @ Predict the State Vector and Covariance Matrix @
        Bp = T*Bu + c;
        Pp = T*Pu*T' + R*Q*R';

        vmask = miss(ones(rows(y),1),1);
        Fmask = miss(ones(rows(y),rows(y)),1);

        if not ismiss(packr(y)); @ Some new data exists @

       @ Flag missing data. See Harvey FSK pp. 144-145 @
  mask = missrv(y,-999999.99) ./= -999999.99;

  W = diagrv(eye(rows(y)),mask);
  W = delif(W,sumc(W) .== 0);

  Z = W * Z;
  d  = W * d;

  @ Predict the next observation on the firm's return @
  yp = Z*Bp + d;

  @ The prediction error @
  v = packr(y) - yp;

  @ The prediction error variance at t @
  F = Z*Pp*Z' + selif(selif(H,mask)',mask)';

tempry = Z;

  @ Update the state vector and Covariance Matrix where
possible @

  Prefix = Pp*Z'invpd(F);
  Bu     = Bp + Prefix*v;
  Pu     = Pp - Prefix*Z*Pp;

  vmask[1:rows(v)] = v;
  v = vmask;

  Fmask[1:rows(F),1:cols(F)] = F;
  F = Fmask;

        else;
  v = vmask;
  F = Fmask;
  Bu = Bp;
  Pu = Pp;
        endif;

        @ Store the Filtration Output @
        if i /= 1;
  Bvec  = Bvec  | Bu';
  Vvec  = VVec  | v';
  Fvec  = FVec  | F;
  PpVec = PpVec | Pp;
  PuVec = PuVec | Pu;
        else;
  Bvec  = Bu';
  Vvec  = v';
  Fvec  = F;
  PpVec = Pp;
  PuVec = Pu;
        endif;

        i = i + 1;
    endo;

    retp(Bvec,Vvec,Fvec,PpVec,PuVec);
endp;

/*********************************************************************
*******
LL : Calculates the loglikelihood Function for the entire set of
observations

INPUT
1.   Parameters.

OUTPUT
1.   Loglikelihood function value.
**********************************************************************
******/
proc(1) = LL(Parms,Data);

local   BVec, @ Kalman Filter Output @
        VVec,
        FVec,
        PpVec,
        PuVec,
   F,
        v,
   n,     @ Number of measurement equations @
        i,
        Logl; @ Log Likelihood function value @

@ Filter the data to extract the prediction error @
    {Bvec,Vvec,Fvec,PpVec,PuVec} = Kalman(Parms,Data);

@ Get the number of the Measurement Variables @
 n = cols(Vvec);

@ Iterate over Observations for Log Likelihood Computation @

   LogL = 0; @ Initialise the vector of log likelihood values @
    i = 1;
    do while i <= rows(Vvec);

        @ Get v and F dumping missing observations @
        v = packr(Vvec[i,.]');
   if not ismiss(v);
  F  = FVec[(i-1)*n+1:(i-1)*n+rows(v),1:rows(v)];

       @ Construct the full Likelihood @
 LogL = LogL | (n*ln(2*Pi*det(F)) + v'invpd(F)*v);
   endif;

        i = i + 1;
    endo;

    retp(-0.5*trimr(LogL,1,0));
endp;

/*********************************************************************
*******
Smoother : Generates the states vector using a Fixed Interval
Smoother.

INPUT:
1.   Parameters.

OUTPUT:
1.   Bvec, the State Vector.
2.   PuVec, the State Vector Covariance matrix.

See Harvey's text on Structural time series models and
the Kalman Filter (FSK), p 154 especially!
**********************************************************************
******/
proc(2) = Smoother(parms,data);
local   Pus, @ Smoothed Covariance Matrix @
        Ps,  @ PStar matrix in Harvey FSK text p. 154 @
        Bvec,
        PpVec,
        PuVec,
        VVec,
        FVec,
        Pp,
        Pu,
        T,h,Q,R,d,c,y,Z,
   k,
        i;

@ Call the Kalman Filter Procedure @
    {BVec,VVec,FVec,PpVec,PuVec} = Kalman(parms,Data);

@ Get the dimension of the state vector @
    k = cols(Bvec);

@ Retrieve the First Smoothed Covariance Matrix for observation T @
    Pus = PuVec[rows(PuVec)-k+1:rows(PuVec),.];

@ Last State Vector in period T is already smoothed as such @
    i = Rows(BVec) - 1;
    do while i >= 1;

        @ Get the System Matrices and Data @
        {T,H,Q,R,d,c,y,Z} = SystemM(parms,Data[i,.]);

        @ Get Pp for i+1/i  which is row i+1 of PpVec @
        Pp = PpVec[k*i+1:k*(i+1),.];

        @ Get Pu for i which is row i of PuVec @
        Pu = PuVec[k*(i-1)+1:k*i,.];

        @ Construct the P Star Matrix @
        Ps = Pu*T'invpd(Pp);

        @ Construct Smooth State Vector for period i @
        BVec[i,.] = (Bvec[i,.]' + Ps * (Bvec[i+1,.]' - T*BVec[i,.]'))';

        Pus = Pu + Ps * (Pus - Pp) * Ps';
        PuVec[k*(i-1)+1:k*i,.] = Pus;

        i = i - 1;
    endo;
    retp(BVec,PuVec);
endp;

