*
* GIBBS.RPF
* Example of Gibbs sampling for a multiple regression with an
* independent prior
*
* RATS User's Guide, Example from Section 16.7.
* 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)