Code: Select all
**************************************************************************************************************
* PROCEDURE: GibbsVAR.src
*
* syntax: @GibbsVAR(options) y stpt endpt prenobs lags prelags ndraws burndraws meanfirstlag Pipriorprec shrinkage fcper &
* PiRes SigmaRes ForecastRes
*
* Written by Todd Clark of the Federal Reserve Bank of Kansas City, drawing on code by Tom Doan
*
* PURPOSE: Using Gibbs sampling to compute posterior distributions of the parameters of a VAR, along with forecasts
* The model/estimator is 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
*
* NOTES: (1) priors for slope coefs of VAR, Pi(L), are handled with modified Litterman approach, without options
*
* INPUT: y vector of series of observed endogenous 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 models 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), as in Litterman/Minneapolis convention
* 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.
* - note: putting loose prior on constant
* fcper number of periods ahead to forecast (set to 0 to skip forecasting)
*
* OUTPUT: PiRes vector of post-burn draws of Pi
* 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:
* printmean/[noprintmean] : display posterior means and standard deviations (default is noprint)
*
* FIRST: 2008-05-09
*
*
*
***********************************
*********************************** declarations
***********************************
procedure GibbsVAR y stpt endpt prenobs lags prelags ndraws burndraws meanfirstlag Pipriorprec shrinkage fcper PiRes SigmaRes ForecastRes
type vec[ser] y
type int stpt endpt prenobs lags prelags ndraws burndraws fcper
type vec meanfirstlag shrinkage
type vec[rec] *PiRes
type vec[symm] *SigmaRes
type rec[ser] *ForecastRes
type symm *Pipriorprec
** start here
option switch prmean 0
* note that what this local stuff is used for becomes clearer below
local int nvar draws_burn ncoef i j l draws nobs rawstpt ystpt availprenobs
local vec[ser] forecasts
local vec univar ux Pihat_ols Pimean Pidraw Pipriormean Pidiff
local symm hx wparm hb Sigma invSigma xxxols xx olssigma
local rec minnmean minnprec Pidraw_mat fx meanmat sdmat Piover
local model varmodel
***********************************
*********************************** dimensioning
***********************************
***** key dimensions
comp nvar = %rows(y)
comp ncoef = nvar*lags+1 ;* number of coefficients in each equation
comp draws_burn = ndraws - burndraws
***** stuff stored, for all post-burn draws
dim PiRes(draws_burn) SigmaRes(draws_burn) ForecastRes(draws_burn,nvar)
***** stuff for Pi computations
dim Pidraw(ncoef*nvar) ;* vector for drawing Pi
dim univar(nvar) ;* for storing univariate st devs used in setting Pi prior standard deviations
dim Pimean(nvar*ncoef) Pidiff(nvar*ncoef) ux(nvar*ncoef) ;* used for drawing innovations
overlay Pidiff(1) with Piover(ncoef,nvar)
***** stuff for Pi prior
comp minnmean = %zeros(ncoef,nvar)
comp minnprec = %zeros(ncoef,nvar) ;* precision^2 = inverse of the variance
***********************************
*********************************** 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
comp minnmean(ncoef,i)=0.0 ;* flat prior on intercept
if availprenobs>=(prenobs+prelags)
comp minnprec(ncoef,i)=1./(1000.*univar(i))
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
}
***********************************
*********************************** MCMC
***********************************
smpl stpt endpt ;* sample for estimation
***** setting up VAR used in loop, and obtaining OLS estimates of VAR and related moments used in Gibbs sampler
system(model=varmodel)
variables y
lags 1 to lags
det constant
kfset xxxols
end(system)
estimate(noprint)
comp Pihat_ols=%vec(%modelgetcoeffs(varmodel))
comp nobs = %nobs
comp olssigma = nobs*%sigma
comp xx = inv(xxxols)
************************************ setting initial values needed in MCMC
comp sigma=%ranwisharti(%decomp(inv(olssigma)),nobs)
do draws=1,ndraws
********************** STEP 1: drawing Pi (VAR slope coefs)
****** drawing slope coefficients, and reassigning them to model
comp invSigma = inv(sigma)
comp hb=%kroneker(invsigma,xx)
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)
********** STEP 2: Draw Sigma
comp Pidiff=Pidraw-Pihat_ols
comp wparm=olssigma+%mqform(xx,Piover) ;* remember that wparm computed this way (as in K-K) is the same as it would be computed as the sum of squared 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 SigmaRes(draws-burndraws) = Sigma
*** construction and storing forecasts
if fcper>0
{
comp %modelsetcoeffs(varmodel,Pidraw_mat)
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}
end do i
}
}
end do draws
if prmean==1
{
dis ''
** 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