/*******   SSM bootstrap procedure   ****************************
Author: Thiery Roncalli       Thierry.Roncalli@montesquieu.u-bordeaux.fr
Date: 10 Sep 1996
This code is provided without guarantee for public non-commercial use.

Here are two procedures for TSM users. KF_Gain could be used
to compute the Kalman filter gains K(t). SSM_bootstrap is
an implementation of the algorithm of Stoffer and Wall [1991],
Bootstrapping state-space models: Gaussian maximum likelihood
estimation and the Kalman filter, Journal of the American
Statistical Association, 86, 1024-1033

The kalman filter gains are given by Harvey's formula (3.2.4b):
   K(t) = T*P(t|t-1)*Z'*inv[F(t)]

Bootstrap for state space models is based on the following
innovations form representation :

   y(t) = Z*a(t|t-1) + d + v(t)
   a(t+1|t) = T*a(t|t-1) + c + K(t)*v(t)

with v(t) the innovation process.


***********************************************************
Here is an example using the SSM_bootstrap procedure.
The example requires the data in reinsel.asc and theta.

new;
library optmum,tsm,pgraph;
TSMset;

load reinsel[100,2] = reinsel.asc;

invest = reinsel[.,1];
invent = reinsel[.,2];
di = invest - lag1(invest);
y = packr(di~invent);

load theta = EXACT_ML;

beta = theta[1:8];
Pchol = (theta[9]~0)|(theta[10]~theta[11]);
SIGMA = Pchol*Pchol';

{Z,d,H,T,c,R,Q} = arma_to_SSM(beta,1,1,SIGMA);
call SSM_build(Z,d,H,T,c,R,Q,0);

y0 = y[1,.]';
{a0,P0} = SSM_ic;
call KFiltering(y,a0,P0);

ys = SSM_bootstrap(y0,a0);

xy(seqa(1,1,rows(y)),y[.,1]~ys[.,1]);
xy(seqa(1,1,rows(y)),y[.,2]~ys[.,2]);



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

/*
** kalman2.src - Time Series Modelling.
** (C) Copyright 1996 by Global Design/Thierry Roncalli.
** All Rights Reserved.
**
**  Format                        Purpose                                 Line
** ============================================================================
**  K = KF_Gain;                  Compute the Kalman filter gains          15
**
**  ys = SSM_bootstrap(y0,a0);    Bootstrap for state space models         30
*/

/*
** KF_Gain
**
** Purpose: Compute the Kalman filter gains
**
** Format:  K = KF_gain;
**
** Output:  K - Nobs*(N*m) matrix, the Kalman filter gains K(t).
**
** Remarks: To read the matrices K(t), use the reshape command
**                       K(t) = reshape(K[i,.]',N,m)';
**
*/

/*
** SSM_bootstrap
**
** Purpose: Bootstrap for state space models
**
** Format:  ys = SSM_bootsrap(y0,a0);
**
** Input:   y0 - N*1 vector, the initial measurement vector
**          a0 - m*1 vector, the initial state vector
**
** Output:  ys - Nobs*N matrix, the bootstrap data set
**
**
*/

/*
** KF_Gain
*/

proc (1) = KF_Gain;
  local N,m,Z,T,Pt_1,Ft_inv,Nobs;
  local i,Kt,K;

  if _SSM_parameters[1] /= 1;

    ERRORLOG "error: Not valid for a time-variant state space model.";
    end;

  endif;

  if _SSM_parameters[5] /= 1;

    ERRORLOG "error: Kalman filter not used.";
    end;

  endif;

  N = _SSM_parameters[2];
  m = _SSM_parameters[3];
  Z = _SSM_Z;
  T = _SSM_T;

  Pt_1 = KF_matrix(6);               /* Read P(t|t-1)  */
  Ft_inv = KF_matrix(8);             /* Read inv(F(t)) */
  Nobs = rows(Pt_1);

  K = miss(zeros(Nobs,m*N),0);

  Z = Z';

  i = 1;
  do until i > Nobs;

    if ismiss(Pt_1[i,.]);

      i = i + 1;
      continue;

    else;

      /* Harvey [1990], formula 3.2.4b */

      Kt = T*xpnd(Pt_1[i,.]')*Z*xpnd(Ft_inv[i,.]');

      K[i,.] = vec(Kt)';
      i = i + 1;

    endif;

  endo;

  retp(K);
endp;

/*
** SSM_bootstrap
*/

proc (1) = SSM_bootstrap(y0,a0);
  local K,Nobs,Ft,vt,N,m,e,P,ys,i,Pchol,et;
  local Z,T,c,d,as,Tstar,Rstar,cstar,Kt;

  K = KF_Gain;
  Nobs = rows(K);

  N = _SSM_parameters[2];
  m = _SSM_parameters[3];

  if ( rows(y0) /= N ) or ( cols(y0) /= 1 );
    ERRORLOG "error: Vector y0 is not conformable.";
    end;
  endif;

  if ( rows(a0) /= m ) or ( cols(a0) /= 1 );
    ERRORLOG "error: Vector a0 is not conformable.";
    end;
  endif;

  if ismiss(K);
    ERRORLOG "error: Not valid with missing values.";
    end;
  endif;

  Ft = KF_matrix(7);             /* Read F(t) */
  vt = KF_matrix(2);             /* Read v(t) */

  e = zeros(Nobs,N);
  P = zeros(Nobs,N^2);
  ys = zeros(Nobs,N);

  /* step 1: standardized innovations */

  i = 1;
  do until i > Nobs;

    Pchol = chol(xpnd(Ft[i,.]'))';
    et = inv(Pchol)*vt[i,.]';

    P[i,.] = vec(Pchol)';
    e[i,.] = et';

    i = i + 1;
  endo;

  /* step 2: bootstrap sample of standardized innovations */

  e = bootstrap(e);

  /* step 3: bootstrap data set of y(t) */

  Z = _SSM_Z; T = _SSM_T; c = _SSM_c; d = _SSM_d;

  Tstar = (zeros(N,N)~Z)|(zeros(m,N)~T);
  cstar = d|c;

  as = y0|a0;

  i = 1;
  do until i > Nobs;

    Kt = reshape(K[i,.]',N,m)';      /* Read the Kalman gain K(t) */
    Pchol = xpnd2(P,i);
    Rstar = Pchol|(Kt*Pchol);

    et = e[i,.]';

    as = Tstar*as + cstar + Rstar*et;

    ys[i,.] = as[1:N]';

    i = i + 1;
  endo;

  retp(ys);
endp;
