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