/*
** This is an example of the use of Gibbs sampling applied to a VAR with
** a standard Minnesota prior. Different priors can be handled by changing
** the way that bprior and hprior (the mean and precision of the prior) are
** created.
**
** Written by Tom Doan, Estima, April 2003.
*/

Between here and next line of ****
a)  Set the four variables (lags, nvar, nstep, ndraws)
b)  Read the data and perform any required transformations
c)  Replace the VARIABLES list in the system definition with your variables
d)  Change the settings of tight and other.
*/
compute lags=4			;*Number of lags
compute nvar=2			;*Number of variables
compute nstep=16		;*Number of response steps
compute ndraws=2500		;*Number of draws
*
cal 1959 1 4
allocate 1998:4
*
open data drisampl.rat
data(format=rats) / gdpq fm1
log gdpq
log fm1
*
* These are the controlling parameters for the symmetric "Minnesota" prior - the
* own lag tightness and the relative tightness on the other lags.
*
compute tight=.1
compute other=.5
*
* First, estimate the system with a prior via mixed estimation, to see if
* we come up with something similar via Gibbs sampling.
*
system(model=priormodel)
variables gdpq fm1
lags 1 to lags
specify(type=symmetric,tight=tight) other
det constant
end(system)
*
estimate
*
* Same system without a prior. This is needed to get the likelihood information
*
system(model=varmodel)
variables gdpq fm1
lags 1 to lags
kfset xxx
det constant
end(system)
*
*************************************************************************************
estimate
*
* Get the estimated standard errors from the regressions for use in scaling the prior
*
dec vect olssee(nvar)
ewise olssee(i)=sqrt(%sigma(i,i))
*
* 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 zero precision (that is, a flat
* prior).
*
compute ncoef=lags*nvar+1
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)=(i==j.and.l==1)
         compute minnprec((j-1)*lags+l,i)=olssee(j)/olssee(i)*%if(i==j,1.0/tight,1.0/(other*tight))
      end do l
   end do j
   compute minnmean(ncoef,i)=0.0,minnprec(ncoef,i)=0.0
end do i
*
* This works better with everything in vec form, rather than as an KxN matrix.
*
ewise minnprec(i,j)=minnprec(i,j)**2
compute hprior=%diag(%vec(minnprec))
compute bprior=%vec(minnmean)
*
* Except that there's one spot where the KxN form is more convenient, so we
* create bover and bdiff as overlayed copies of each other.
*
dec rect bover
dec vect bhat bmean(nvar*ncoef) bdiff(nvar*ncoef) ux(nvar*ncoef)
overlay bdiff(1) with bover(ncoef,nvar)
*
* Pull out the information from the OLS regression
*
compute sigma=%sigma
compute xx   =inv(xxx)
compute bhat=%vec(%modelgetcoeffs(varmodel))
*
* The loop below does the Gibbs sampler. All I'm doing here is accumulating the means of
* the coefficients. However, the same basic code could easily be adapted for doing
* Monte Carlo integration for impulse response functions or forecasts.
*
dec vect gibbsmean(nvar*ncoef)
compute gibbsmean=%const(0.0)
do draws=1,ndraws
*
*     Draw b given sigma
*
   compute hb=%kroneker(inv(sigma),xx)
   compute hx=inv(hb+hprior)
   compute bmean=hx*(hprior*bprior+hb*bhat)
   compute fx=%decomp(hx)
   compute bdraw=bmean+fx*(ux=%ran(1.0))
*
*     The draw for the coefficients is in bdraw as a vector
*
   compute gibbsmean=gibbsmean+bdraw
*
*     Draw sigma given b
*
   compute bdiff=bdraw-bhat
   compute wparm=%sigma+(1.0/%nobs)*%mqform(xx,bover)
   compute sigma=%mqform(%nobs*inv(%ranwishart(nvar,%nobs-ncoef)),tr(%decomp(wparm)))
end do draws
*
* Postprocessing
*
compute gibbsmean=(1.0/ndraws)*gibbsmean
disp gibbsmean