* * 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)