* * GIBBSVAR.RPF * Example of Gibbs sampling on a Bayesian VAR with a "Minnesota" prior * compute lags=4 ;*Number of lags compute nstep=16 ;*Number of response steps compute nburn=500 ;*Number of burn-in draws compute ndraws=2500 ;*Number of keeper draws * 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) * * 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=varmodel) variables loggdp loginv logc ftb3 lags 1 to lags specify(type=symmetric,tight=tight) other det constant end(system) ************************************************************************************ * * Estimate with a prior * estimate * * Same system without a prior. This is needed to get the likelihood * information. * estimate(ols,sigma) compute nvar=%nvar * * 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 * ewise minnprec(i,j)=minnprec(i,j)^2 compute hprior=%diag(%vec(minnprec)) compute bprior=minnmean * compute sigmad=%sigma cmom(model=varmodel) * dec vect[series] forecast(nvar) dec vect[series] forestderr(nvar) compute fstart=%regend()+1 compute fend =fstart+nstep-1 * do i=1,nvar set forecast(i) fstart fend = 0.0 set forestderr(i) fstart fend = 0.0 end do i * infobox(action=define,progress,lower=-nburn,upper=ndraws) "Gibbs Sampler" do draw=-nburn,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<=0 next compute %modelsetcoeffs(varmodel,bdraw) simulate(model=varmodel,cv=sigmad,results=simresults,steps=nstep) do i=1,nvar set forecast(i) fstart fend = forecast(i)+simresults(i) set forestderr(i) fstart fend = forestderr(i)+simresults(i)^2 end do i end do draw infobox(action=remove) * do i=1,nvar set forecast(i) fstart fend = forecast(i)/ndraws set forestderr(i) fstart fend = sqrt(forestderr(i)/ndraws-forecast(i)^2) end do i * do i=1,nvar set lower fstart fend = forecast(i)-2.0*forestderr(i) set upper fstart fend = forecast(i)+2.0*forestderr(i) graph(header="Forecasts of "+%modellabel(varmodel,i)) 4 # forecast(i) fstart fend # lower fstart fend 2 # upper fstart fend 2 # %modeldepvars(varmodel)(i) fstart-12 fstart-1 end do i