*
* GIBBS.RPF
* RATS Version 8, User's Guide, Example 16.4
* Multiple linear regression with independent prior
* Adapted from Koop, Bayesian Econometrics, 4.2.7 from pp 73-77.
*
open data hprice.prn
data(format=prn,org=columns) 1 546 price lot_siz bed bath stor drive $
 recroom bment gas aircond gar desireloc
*
* This is the equation we're analyzing
*
linreg price
# constant lot_siz bed bath stor
*
* Prior for variance.
*
compute s2prior=5000.0^2
compute nuprior=5.0
*
* Prior mean and variance. Convert the variance to precision by
* inverting.
*
compute [vector] bprior=||0.0,10.0,5000.0,10000.0,10000.0||
compute [symm]   vprior=$
   %diag(||10000.0^2,5.0^2,2500.0^2,5000.0^2,5000.0^2||)
compute [symm]   hprior =inv(vprior)
*
* Compute the cross product matrix from the regression. This will
* include the dependent variable as the final row. The cross product
* matrix has all the sufficient statistics for the likelihood (except
* the number of observations, which is %nobs).
*
cmom(lastreg)
*
* The two obvious places to start the Gibbs sampler are the OLS
* estimates and the prior. We'll use OLS.
*
compute bdraw =%beta
compute s2draw=%seesq
*
* We'll use the following for all MCMC estimators. This ends up taking
* nburn+1 burn-ins.
*
compute nburn =1000
compute ndraws=10000
*
dec series[vect] bgibbs
dec series       hgibbs
gset bgibbs 1 ndraws = %zeros(%nreg,1)
set  hgibbs 1 ndraws = 0.0
*
do draw=-nburn,ndraws
   *
   * Draw residual precision conditional on previous beta
   *
   compute rssplus=nuprior*s2prior+%rsscmom(%cmom,bdraw)
   compute hdraw  =%ranchisqr(nuprior+%nobs)/rssplus
   *
   * Draw betas given hdraw
   *
   compute bdraw  =%ranmvpostcmom(%cmom,hdraw,hprior,bprior)
   if draw<=0
      next
   *
   * Do the bookkeeping here.
   *
   compute bgibbs(draw)=bdraw
   compute hgibbs(draw)=hdraw
end do draw
*
@mcmcpostproc(ndraws=ndraws,mean=bmean,$
    stderrs=bstderrs,cd=bcd,nse=bnse) bgibbs
report(action=define)
report(atrow=1,atcol=1,align=center) "Variable" "Coeff" $
    "Std Error" "NSE" "CD"
do i=1,%nreg
   report(row=new,atcol=1) %eqnreglabels(0)(i) bmean(i) $
     bstderrs(i) bnse(i) bcd(i)
end do i
report(action=format,atcol=2,tocol=3,picture="*.###")
report(action=format,atcol=4,picture="*.##")
report(action=show)