* * SHILLERGIBBS.RPF * Gibbs sampling analysis of Shiller Smoothness Prior * open data haversample.rat calendar(m) 1947 data(format=rats) 1947:1 2007:4 fltg ftb3 set shortrate = ftb3 set longrate = fltg * * Set the number of draws. * Set the prior degrees of freedom (NU) and mean (S2) for the precision * (reciprocal of the variance) for the regression residuals. * Set the prior precision for the second difference in the coefficients. * compute nburn=100 compute ndraws=5000 compute s2=.50^2.0 compute nu=4.0 compute hb=1.0/(.03^2) * * Generate the precision matrix for the prior. For this example, the * prior is flat on the constant and zero lag coefficient, and the second * differences on the remainder of the lag polynomial are given * independent Normal priors with common precision. In general, PRIORH * will be the (generalized) inverse of your prior variance. * declare rect dummy(22,26) declare symm priorh(26,26) fmatrix(diff=2) dummy 1 3 compute priorh=hb*(tr(dummy)*dummy) compute priorb=%zeros(26,1) * * Compute the cross-moment matrix * cmom # constant shortrate{0 to 24} longrate * * Compute the OLS regression to provide an initial value for the beta * vector. * linreg(cmom) longrate # constant shortrate{0 to 24} compute beta=betaols=%beta * * Specialized code for the bookkeeping is needed in this example We are * keeping track of the following: * * 1. Full distribution on the intercept * 2. Full distribution on the zero lag coefficient * 3. Full distribution of the sum of the lag coefficients * 4. The first and second moments for all coefficients * * 1,2 and 3 each require a series with length equal to the number of * draws. To allow us to compute quickly the sum of the lag coefficients, * we generate a vector (called SUMMER) which, when dotted with a draw * for the coefficients, will give the sum of the lag coefficients. * * 4 requires simply a pair of series with length equal to the number of * coefficients, all elements initialized to zero. * dec vect summer(%nreg) ewise summer(i)=%if(i>=2,1.0,0.0) * set inter 1 ndraws = 0.0 set coeff1 1 ndraws = 0.0 set sums 1 ndraws = 0.0 * set comoment1 1 %nreg = 0.0 set comoment2 1 %nreg = 0.0 * infobox(action=define,progress,lower=-nburn,upper=ndraws) "Gibbs Sampler" do draw=-nburn,ndraws infobox(current=draw) * * Draw residual precision conditional on previous beta * compute rssplus=nu*s2+%rsscmom(%cmom,beta) compute hu =%ranchisqr(nu+%nobs)/rssplus * * Draw betas given hu * compute beta=%ranmvpostcmom(%cmom,hu,priorh,priorb) if draw<=0 next * * Do the bookkeeping here. For 1,2,3, stuff the function of the * current draw into the "DRAW" element of the result series. * compute inter(draw)=beta(1) compute coeff1(draw)=beta(2) compute sums(draw)=%dot(beta,summer) * * Accumulate the sum and sum of squares of the coefficients * set comoment1 1 %nreg = comoment1+beta(t) set comoment2 1 %nreg = comoment2+beta(t)^2 end do draw infobox(action=remove) density(smoothing=1.5) inter 1 ndraws ginter finter density(smoothing=1.5) coeff1 1 ndraws gcoeff1 fcoeff1 density(smoothing=1.5) sums 1 ndraws gsums fsums scatter(style=lines,footer="Posterior for Intercept") # ginter finter scatter(style=lines,footer="Posterior for Lag 0") # gcoeff1 fcoeff1 scatter(style=lines,footer="Posterior for Sum") # gsums fsums * set comoment1 1 %nreg = comoment1/ndraws set comoment2 1 %nreg = sqrt(comoment2/ndraws-comoment1^2) set upper 1 %nreg = comoment1+2.0*comoment2 set lower 1 %nreg = comoment1-2.0*comoment2 set ols 1 %nreg = betaols(t) * graph(number=0,footer="Graph of Lag Distribution") 4 # comoment1 2 %nreg # upper 2 %nreg 2 # lower 2 %nreg 2 # ols 2 %nreg