**************************************************************************************************************                                                   
* PROCEDURE:  condfcbyKF.src
*
* syntax: @condfcbyKF y ycond stpt endpt lags ndraws fcper PiRes SigmaRes ForecastRes
*
* PURPOSE:  Generating conditional forecasts using Kalman filter
*           producing draws of forecasts reading in draws of VAR coefs and error var-cov matrix
*           Kalman filter implementation follows approach outlined in Banbura, Giannone, and Lenza (2014),
*             as implemented in Tom Doan's documentation of the use of DLM for forecasting missing data
*
* INPUT:    y                  vector of series of observed endogenous variables
*           ycond              time series of vec corresponding to y filled with NAs and conditions where appropriate
*           stpt               starting point of estimation (meaning observations back to stpt-lags) serve as initial values
*           endpt              end point of estimation sample
*           lags               Lag length of the VAR.
*           ndraws             number of total simulation draws (ndraws=1 to just compute forecasts based on posterior mean coefs)
*           fcper              number of periods ahead to forecast
*           PiRes              vector of simulated draws of matrix Pi
*           SigmaRes           vector of simulated draws of Sigma
*           ForecastRes        rect array of time series of forecasts, covering posterior draws and num of variables
*
* AUTHOR:   Todd Clark, Federal Reserve Bank of Cleveland, todd.clark@clev.frb.org

***********************************
*********************************** declarations
***********************************
procedure condfcbyKF y ycond stpt endpt lags ndraws fcper PiRes SigmaRes ForecastRes
type vec[ser] y
type int stpt endpt lags ndraws *fcper
type vec[rec] PiRes 
type vec[symm] SigmaRes  
type rec[ser] *ForecastRes
type ser[vec] ycond

local int nvar draws statedim stpt i time
local model varmodel 
local ser[vec] ystsp xstsp states
local ser[symm] vstates
local ser[symm] svmat 
local symm swmat
local rec amat cmat
local vec zvec
local vec[rec] dlmstuff

****************** set up, mostly for st space form of VAR
comp nvar = %rows(y)
comp statedim = nvar*lags
dim forecastres(ndraws,nvar)

gset ystsp endpt+1 endpt+fcper = ycond                ;* meas vector (holding conditions)
gset ystsp stpt-lags endpt = %xt(y,t)
gset xstsp stpt endpt+fcper = %zeros(statedim,1)      ;* state vector (will hold actual data in history)

comp cmat = %zeros(statedim,nvar)                    ;* coef matrix in meas equation
do i = 1,nvar
 comp cmat(i,i) = 1.
end do i
comp swmat = %zeros(statedim,statedim)                     ;* var-cov matrix of innovations in state equation
comp zvec = %zeros(statedim,1)                             ;* for intercepts in state equation

/*  

this block below gives forecast confidence intervals that are way narrower than I get with DLS approach; just leaving at all zeros gives the same
narrow confidence intervals:  

gset svmat stpt endpt = %zeros(nvar,nvar)                 ;* for estimation sample, the variance of innovations to the meas equation = 0
gset svmat endpt+1 endpt+fcper = %fill(nvar,nvar,%NA)
do time = endpt+1,endpt+fcper
 do i = 1,nvar
  if %valid(ystsp(time)(i))==1
   comp svmat(time)(i,i) = 0.0   
 end do i
end do time
*/

*** need to make measurement error large for variables not observed to match results against DLS approach
gset svmat endpt+1 endpt+fcper = 1000000000000.^2*%identity(nvar) 
do time = endpt+1,endpt+fcper
 do i = 1,nvar
  if %valid(ystsp(time)(i))==1
   comp svmat(time)(i,i) = 0.0   
 end do i
end do time

** need to define VAR for use below
system(model=varmodel)
variables y                  
lags 1 to lags
det constant
end(system)

****************** calculations: loop over draws of coefs to form forecasts
do draws = 1,ndraws

 ** assign draws of coefficients to state space params
 comp %modelsetcoeffs(varmodel,pires(draws))         ;* assign coefs to defined VAR, for use below in computations
 comp amat = %modelcompanion(varmodel)               ;* companion form slope coefficients into coef matrix of state equation
 comp %psubmat(swmat,1,1,SigmaRes(draws))            ;* draw of error var-cov matrix in variance of state equation
 comp %psubmat(zvec,1,1,tr(%xrow(pires(draws),%rows(pires(draws)))))  ;* intercepts

 ** compute initial value of state vector and variance
 comp dlmstuff = %dlminit(amat,swmat,%identity(statedim),zvec)

 ** run just Kalman filter through data sample
 dlm(c=cmat,a=amat,y=ystsp,sw=swmat,z=zvec,sx0=dlmstuff(1),x0=dlmstuff(2),type=filter) stpt endpt states vstates

 ** now do simulation smoother to draw conditional forecasts
 dlm(c=cmat,a=amat,y=ystsp,sv=svmat,sw=swmat,z=zvec,sx0=vstates(endpt),x0=states(endpt),type=csimulate) endpt+1 endpt+fcper states
 
 ** now store forecasts
 do i = 1,nvar
  set forecastres(draws,i) endpt-(%calendar()(6)) endpt = y(i){0}   ;* back-fill the forecasts with actual data, for computation of 4Q or 12M averages of some variables
  set forecastres(draws,i) endpt+1 endpt+fcper = states(t)(i)
 end do i
 
end do draws

end