The RATS Software Forum

Code: Select all

**************************************************************************************************************                                                   
* PROCEDURE:  SSVAR.src
*
* syntax: @SSVAR(options) y d stpt endpt prenobs lags prelags ndraws burndraws meanfirstlag Pipriorprec shrinkage Psipriormean Psistdev fcper &
*                   PiRes PsiRes SigmaRes ForecastRes
*
*           Written by Todd Clark of the Federal Reserve Bank of Kansas City, drawing on code by Tom Doan 
*
* PURPOSE:  Using Gibbs sampling as in Villani, Mattias (2006), ``Inference in Vector Autoregressive Models with an Informative Prior 
*            on the Steady State,'' Sveriges Riksbank Research Paper No.\ 19, forthcoming in the Journal of Applied Econometrics.
*            Program generates posterior estimates (including forecasts, if the user chooses)
*            The model/estimator is similar to the normal-diffuse of Kadiyala and Karlsson (1997), "Numerical Methods for Estimation and 
*            Inference in Bayesian VAR Models," Journal of Applied Econometrics, pp 99-132, extended to allow a prior on the steady state.  
*            The model takes the form 
*            Pi(L)*[y(t)-Psi*d(t)] = eps(t),     eps(t)~iid N(0,Sigma).
*            where:  Pi(L)=I-Pi2*L-...-Pik*L^k, and the prior on Sigma is diffuse.
*
* NOTES:  (1) priors for slope coefs of VAR, Pi(L), are handled with modified Litterman approach, without options
*         (2) priors for the steady state coefs Psi, are handled as in Villani
*
* INPUT:    y                  vector of series of observed endogenous variables
*           d                  vector of series of exogenous (deterministic) variables
*           stpt               starting point of estimation (meaning observations back to stpt-lags) serve as initial values
*           endpt              end point of estimation sample
*           prenobs            number of observations in pre-sample to use for setting priors.  If the actual number available
*                                falls short of prenobs, the priors (in this case, just the prior variance on Pi) are pulled
*                                from values that need to be defined from outside the procedure.
*           lags               Lag length of the VAR.
*           prelags            lag of AR model used in pre-sample estimation that pins down some priors
*           ndraws             number of total posterior draws
*           burndraws          number of burnin draws
*           MeanFirstLag       Prior mean vector for first own lag (Pi). Only used for the Litterman prior.
*           Pipriorprec        If there are enough data in the training sample, the prior precision is estimated based on the pre-sample;
*                                  if there isn't enough data, the prior precision needs to be given (input) in Pipriorprec
*           Shrinkage          Hyperparameters in prior.  Smaller values produces priors which are more concentrated around zero.
*                              - Shrinkage(1): Baseline tightness parameter.        
*                              - Shrinkage(2): Own/others shrinkage                 
*                              - Shrinkage(3): Lag length shrinkage parameter.      
*           Psipriormean       prior mean for Psi (vector)
*           Psistdev           prior standard deviation for Psi (vector)
*           fcper              number of periods ahead to forecast (set to 0 to skip forecasting)
*
* OUTPUT:   PiRes              vector of post-burn draws of Pi
*           PsiRes             vector of post-burn draws of Psi
*           SigmaRes           vector of post-burn draws of Sigma
*           ForecastRes        rect array of time series of forecasts, covering post-burn draws and num of variables

*   -Options:
*     printinit/[noprintinit] : display initial priors/values (default is noprint)
*     printmean/[noprintmean] : display posterior means and standard deviations (default is noprint)
*
*
* FIRST:    2008-05-09
*
*

***********************************
*********************************** declarations
***********************************
procedure SSVAR y d stpt endpt prenobs lags prelags ndraws burndraws meanfirstlag Pipriorprec shrinkage Psipriormean Psistdev $
  fcper PiRes PsiRes SigmaRes ForecastRes
type vec[ser] y d
type int stpt endpt prenobs lags prelags ndraws burndraws fcper
type vec meanfirstlag shrinkage Psipriormean Psistdev
type vec[rec] *PiRes *PsiRes
type vec[symm] *SigmaRes  
type rec[ser] *ForecastRes
type symm *Pipriorprec

** start here

option switch prinit 0
option switch prmean 0

* note that what this local stuff is used for becomes clearer below
local int nvar ndet draws_burn ncoef i j l draws row col nobs rawstpt ystpt availprenobs 
local index dvarlist
local vec[ser] ydtr residuals dminus forecasts
local rec[ser] tempresiduals
local vec univar ux_psi ux Pihat_ols Pimean Pidraw Pipriormean Psidraw Psimean
local symm psipriorv invomegaprior hx Omegamean wparm hb Sigma invSigma dd xxxols
local rec minnmean minnprec uprime Pidraw_mat Psidraw_mat fx fx_Omega meanmat sdmat ydcov
local rec pimat
local model varmodel ssmodel

***********************************
*********************************** dimensioning
***********************************
***** key dimensions
comp nvar = %rows(y)
comp ndet = %rows(d)
comp ncoef = nvar*lags              ;* number of slope coefficients in each equation
comp draws_burn = ndraws - burndraws 
comp nobs = endpt-stpt+1

***** stuff stored, for all post-burn draws
dim PiRes(draws_burn) PsiRes(draws_burn) SigmaRes(draws_burn) ForecastRes(draws_burn,nvar)

***** stuff for Pi computations
dim univar(nvar)  ;* for storing univariate st devs used in setting Pi prior standard deviations
dim ux(nvar*ncoef) ;* used for drawing innovations
dim Pidraw(ncoef*nvar) ;* vector for drawing Pi
***** stuff for Pi prior (Litterman mean and variance)
comp minnmean = %zeros(ncoef,nvar) 
comp minnprec = %zeros(ncoef,nvar) ;* precision^2 = inverse of the variance

***** stuff for Psi computations
dim ux_psi(nvar*ndet) ydtr(nvar) residuals(nvar) dminus(ndet) 
comp psipriorv = %diag(Psistdev.*Psistdev)
comp invomegaprior = inv(psipriorv)
comp uprime = %zeros(nvar*ndet,(lags+1)*nvar*ndet)
comp %psubmat(uprime,1,1,%identity(nvar*ndet))
comp pimat = %zeros(nvar,nvar)

***********************************
*********************************** setting priors (Minneapolis prior on slope coefs Pi) from pre-sample estimates
*********************************** Note:  precision is only set if the available sample size (availprenobs)
*********************************** is at least as great as specified by the inputs (prenobs+prelags).
*********************************** If the sample isn't big enough, precision is left at the value provided in the
*********************************** procedure call.
comp rawstpt = 0
do i = 1,nvar
 inquire(series=y(i)) ystpt *
 comp rawstpt = %imax(rawstpt,ystpt)
end do i
comp availprenobs = stpt-1 - rawstpt + 1

if availprenobs>=(prenobs+prelags)
 {
  do i = 1,nvar
    linreg(noprint) y(i) stpt-prenobs stpt-1
    # constant y(i){1 to prelags}
    comp univar(i)=sqrt(%seesq)  ;* %rss/%nobs
  end do i
 }

***** setting mean and precision
do i=1,nvar
 do j=1,nvar
  do l=1,lags
   comp minnmean((j-1)*lags+l,i)=(i==j.and.l==1)*MeanFirstLag(i)
   if availprenobs>=(prenobs+prelags)
    comp minnprec((j-1)*lags+l,i)=univar(j)/univar(i)*%if(i==j,float(l)**shrinkage(3)/shrinkage(1),float(l)**shrinkage(3)/(shrinkage(2)*shrinkage(1)))
  end do l
 end do j
end do i

comp Pipriormean=%vec(minnmean)         ;* prior mean, in vector form
if availprenobs>=(prenobs+prelags)
 {
  ewise minnprec(i,j)=minnprec(i,j)**2
  comp Pipriorprec=%diag(%vec(minnprec))  ;* inverse of prior variance matrix
 }

************************************ setting initial values needed in MCMC
comp Psidraw = Psipriormean + %ranmvnormal(%decomp(psipriorv))   ;* initial draw of Psi
cmom(center) stpt endpt
# y
comp sigma=%ranwisharti(%decomp(inv(.9*%cmom)),nobs)    ;* initial draw of Sigma, using .9*var-cov of y as scale matrix

if prinit==1
 {
  dis ''
  dis 'initial Psi = ' Psidraw
  dis ''
  dis 'initial Sigma = ' Sigma
 }

***********************************
*********************************** more set up for MCMC
***********************************
smpl stpt endpt  ;* sample for estimation

***** set up D matrix for steady state prior, outside loop since it is just based on fixed data
comp dvarlist = %rlempty()
enter(varying) dvarlist
# dvarlist d
do j = 1,ndet
  set dminus(j) stpt-lags endpt = -1.*d(j){0}
end do j
do i = 1,lags
 do j = 1,ndet
  enter(varying) dvarlist
  # dvarlist dminus(j){i}
 end do j
end do i
cmom(matrix=dd,nocenter)
# dvarlist

***** model in just deterministic variables, used in loop to demean y variables
system(model=ssmodel)
variables y
det d
end(system)
** now need to assign coefs and run through filter to define demeaned y (ydtr) for first draw
comp Psidraw_mat = tr(%vectorect(Psidraw,nvar))
comp %modelsetcoeffs(ssmodel,Psidraw_mat)
steps(model=ssmodel,results=tempresiduals) * (endpt-(stpt-lags)+1) stpt-lags
do i = 1,nvar
  set ydtr(i) stpt-lags endpt = tempresiduals(2,i)
end do i

***** demeaned VAR used in loop, in demeaned y variables
system(model=varmodel)
variables ydtr
lags 1 to lags
kfset xxxols
end(system)

** note:  because using the option (model=varmodel) with estimate() causes MacRATS v7 to crash, I
**        don't use it.  As a result, to be sure the estimate() is estimating the right model, it has to 
**        be the case that, in the above defs, varmodel is defined after ssmodel

***********************************
*********************************** MCMC
***********************************

do draws=1,ndraws

 ********************** STEP 1:  drawing Pi (VAR slope coefs)
 ****** forming OLS estimates and moments for ydtr = y - Psi*d
 estimate(noprint) ;* ,  note:  including the model option causes a crash in v7 of RATS
 comp Pihat_ols=%vec(%modelgetcoeffs(varmodel))
 ****** drawing slope coefficients, and reassigning them to model
 comp invSigma = inv(sigma)
 comp hb=%kroneker(invsigma,inv(xxxols))
 comp hx=inv(hb+Pipriorprec)
 comp Pimean=hx*(Pipriorprec*Pipriormean+hb*Pihat_ols)
 comp fx=%decomp(hx)
 comp Pidraw = Pimean+fx*(ux=%ran(1.0))
 comp Pidraw_mat = %vectorect(Pidraw,ncoef)
 comp %modelsetcoeffs(varmodel,Pidraw_mat)

 ********************** STEP 2:  drawing Psi (steady state coefficients)
 ** first redefine ydtr to be raw data
 do i = 1,nvar
  set ydtr(i) stpt-lags endpt = y(i){0}
 end do i
 ** forming Y series
 steps(model=varmodel,results=tempresiduals) 
 do i = 1,nvar
   set residuals(i) = tempresiduals(2,i)
 end do i
 cmom(nocenter)
 # residuals dvarlist
 comp ydcov = %xsubmat(%cmom,1,nvar,nvar+1,%cols(%cmom))
 ** setting up U matrix
 do i = 1,lags
  ewise pimat(row,col) = Pidraw_mat((col-1)*lags+i,row)
  comp %psubmat(uprime,1,(ndet*nvar)*i+1,%kroneker(%identity(ndet),tr(pimat)))
 end do i

 ** posterior computations
 comp omegamean = inv(%mqform(%kroneker(dd,invsigma),tr(uprime)) + invomegaprior)
 comp Psimean = omegamean*(uprime*%vec(invsigma*ydcov) + invomegaprior*Psipriormean)
 comp fx_omega = %decomp(omegamean)
 comp Psidraw=Psimean+fx_omega*(ux_psi=%ran(1.0))
 comp Psidraw_mat = tr(%vectorect(Psidraw,nvar))
 comp %modelsetcoeffs(ssmodel,Psidraw_mat)

 ********** STEP 3: Draw Sigma
 **********  (first filter data to get residuals:  first subtract out means, and then run through demeaned VAR,
 **********    then compute residual var-cov, then take Sigma draw)
 steps(model=ssmodel,results=tempresiduals) * (endpt-(stpt-lags)+1) stpt-lags
 do i = 1,nvar
   set ydtr(i) stpt-lags endpt = tempresiduals(2,i)
 end do i
 steps(model=varmodel,results=tempresiduals) ;* filtering demeaned y's to get residuals
 do i = 1,nvar
   set residuals(i) = tempresiduals(2,i)
 end do i
 cmom(matrix=wparm)
 # residuals
 comp sigma=%ranwisharti(%decomp(inv(wparm)),nobs)
   
 ******* storage of MCMC draws, and construction of foreasts for storage
 if draws>burndraws
  {
   *** storing everything but forecasts
   comp PiRes(draws-burndraws) = Pidraw_mat
   comp PsiRes(draws-burndraws) = Psidraw
   comp SigmaRes(draws-burndraws) = Sigma
   *** construction and storing forecasts:  (1) generate forecasts of (y - steady states), and (2) add back in steady state
   if fcper>0
    {
     simulate(model=varmodel,results=forecasts) * fcper endpt+1 sigma  
     do i = 1,nvar
      set forecastRes(draws-burndraws,i) endpt+1 endpt+fcper = forecasts(i){0} + %dot(%eqncoeffs(%modeleqn(ssmodel,i)),%eqnxvector(%modeleqn(ssmodel,i),t))
     end do i
    }
  }

end do draws

if prmean==1
{
dis ''

comp meanmat = %zeros(%rows(PsiRes(1)),%cols(PsiRes(1)))
comp sdmat = %zeros(%rows(PsiRes(1)),%cols(PsiRes(1)))
do i = 1,draws_burn
 comp meanmat = meanmat + PsiRes(i)
 comp sdmat = sdmat + PsiRes(i).*PsiRes(i)
end do i
comp meanmat = meanmat*(1./draws_burn)
comp sdmat = %sqrt(sdmat*(1./draws_burn) - meanmat.*meanmat)

dis 'Psi mean ='
dis ###.### meanmat
dis 'Psi st dev ='
dis ###.### sdmat

** MCMC mean and standard deviation for Pi
comp meanmat = %zeros(%rows(PiRes(1)),%cols(PiRes(1)))
comp sdmat = %zeros(%rows(PiRes(1)),%cols(PiRes(1)))
do i = 1,draws_burn
 comp meanmat = meanmat + PiRes(i)
 comp sdmat = sdmat + PiRes(i).*PiRes(i)
end do i
comp meanmat = meanmat*(1./draws_burn)
comp sdmat = %sqrt(sdmat*(1./draws_burn) - meanmat.*meanmat)

dis 'Pi mean ='
dis ###.### meanmat
dis 'Pi st dev ='
dis ###.### sdmat

** MCMC mean and standard deviation for Sigma
comp meanmat = %zeros(%rows(SigmaRes(1)),%cols(SigmaRes(1)))
comp sdmat = %zeros(%rows(SigmaRes(1)),%cols(SigmaRes(1)))
do i = 1,draws_burn
 comp meanmat = meanmat + SigmaRes(i)
 comp sdmat = sdmat + SigmaRes(i).*SigmaRes(i)
end do i
comp meanmat = meanmat*(1./draws_burn)
comp sdmat = %sqrt(sdmat*(1./draws_burn) - meanmat.*meanmat)

dis 'Sigma mean ='
dis ###.### meanmat
dis 'Sigma st dev ='
dis ###.### sdmat

}

end