I wrote some code for a Gibbs sampler for a VAR with Minnesota Prior and exogenous variables.
I used the haversample data.
Could you have a look at the code and tell me whether it is fine?
Code: Select all
************************
********* Data *********
************************
open data haversample.rat
cal(q) 1959
data(format=rats) 1959:1 2006:4 ftb3 gdph ih cbhm
set loggdp = log(gdph)
set loginv = log(ih)
set logc = log(cbhm)
***************************************
************ General Setup ************
***************************************
compute lags=4 ;*Number of lags
compute nstep=50 ;*Number of response steps
compute ndraws=1000 ;*Number of keeper draws
compute nburn=500 ;*Number of burn-in draws
compute keep = ndraws - nburn ;*Number of kept draws
compute nvari = 3 ;*Number of endogebous variables
compute nexog = 1 ;*Number of exogenous variables
compute ncoef = nvari*lags+nexog*lags+1 ;*Number of coefficients per equation, +1 for intercept +12 for each exogenous variable
***********************************
************ VAR Setup ************
***********************************
*set time = t
compute tight= 0.2
compute other= 0.5
system(model=var)
variables loggdp loginv logc
lags 1 to lags
det constant ftb3{1 to lags}
end(system)
estimate(noprint,ols,sigma) ;*for initial draw of VAR coefficients
impulse(noprint,model=var,steps=nstep,results=impulses) ;*initial set of of impulse responses
compute start = %REGSTART() ;* Starting Period for Structural Residual Estimation
compute end = %REGEND() ;* Ending Period for Structural Residual Estimation
compute nvar=%nvar ;* Number of equations estimated
declare vect olssee(nvar) ;* vector of variances, with standard errors estimated from univariat AR(1) models
compute i = 1
dofor x = loggdp loginv logc
linreg(noprint) x
# constant x{1}
compute olssee(i) = sqrt(%SEESQ)
compute i = i+1
end do
declare vect[rect] %%responses(keep) ;* Storage Vector of Rectanculars for Impulse Response Draws ;* Storage Series of vectors for Kilian IRF Projection
declare series[rect] bgibbs ;* Storage Series of Rectangulars to save the VAR Coefficient Draws
gset bgibbs 1 keep = %zeros(ncoef,nvari) ;* To dimension the storage matrix inside the series
***********************************
********* Minnesota Prior *********
***********************************
* Create the prior mean and (square root of) precision. The olssee
* values are used to scale the precision. The prior mean is 1 for the
* first own lag in each equation, and 0 otherwise. The constant is given
* a mean of zero with a standard deviation according to Meinen and Roehe: lambda4*olssee(i), lambda4=100
dec rect minnmean(ncoef,nvar) minnprec(ncoef,nvar)
do i=1,nvar
do j=1,nvar
do l=1,lags
compute minnmean((j-1)*lags+l,i)=%if(i==j.and.l==1,1,0) ;* sets first on log to 1 and everything else to zero
compute minnprec((j-1)*lags+l,i)=(olssee(j)/olssee(i))*%if(i==j,l/tight,l/(other*tight)) ;* computes the corresponding minnesota precisions
end do l
end do j
do k = 0,nexog*lags
compute minnmean(ncoef-k,i)=0.0,minnprec(ncoef-k,i)=100*olssee(i) ;* sets the minnesota mean and precision for the exogenous regressors
end do k
end do i
ewise minnprec(i,j)=minnprec(i,j)^2
compute hprior=%diag(%vec(minnprec))
compute bprior=minnmean
compute sigmad=%sigma
cmom(model=var) ;* Cmom computes the cross moment matrix %cmom, which is the X'X matrix. It is deterministic and data dependent.
*********************************
********* Gibbs Sampler *********
*********************************
infobox(action=define,progress,lower=1,upper=ndraws) "Gibbs Sampler"
do draw= 1,ndraws
infobox(current=draw)
*
* Draw b given sigma
*
compute bdraw =%ranmvkroncmom(%cmom,inv(sigmad),hprior,bprior)
compute rssmat=%sigmacmom(%cmom,bdraw)
*
* Draw sigma given b
*
compute sigmad=%ranwisharti(%decomp(inv(rssmat)),%nobs)
if draw<=nburn
next
compute %modelsetcoeffs(var,bdraw)
impulse(noprint,model=var,results=impulses,steps=nstep,flatten=%%responses(draw-nburn))
* Bookkeeping:
compute bgibbs(draw-nburn) = bdraw ;* Saves the Draws of the VAR Coefficients
end do draw
infobox(action=remove)
***********************
********* IRF *********
***********************
@mcgraphirf(model=var,center=median,percent=||0.16,0.84||)
Best
Jules