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