* Lag.prg is a RATS program to estimate a regression with a 
* distributed lag using smoothness priors (Robert J. Shiller,
* Econometrica, vol. 41, no. 4, July 1973, p. 775--88) implemented by 
* Gibbs sampling (A. E. Gelfand and A. F. M. Smith, JASA, vol. 85, 1990,
* pp. 398--409).  The code, apart from some annotation and graphics output,
* is taken from gibbs.prg, which  comes with RATS version 5 and is discussed
* in the User's Guide, example 13.2.  The program uses monthly data for 
* Jan. 1947 - March 1999, stored in the file rates.rat.
*
*
CAL 47 1 12
ALLOCATE 0 99:3
OPEN DATA rates.rat
DATA(FORMAT=RATS) / SHORTRATE LONGRATE
*
*  Compute OLS estimates to start
*
LINREG LONGRATE
# CONSTANT SHORTRATE{0 TO 24}
summarize
# shortrate{0 to 24}
summarize
# shortrate{1 to 24}
summarize
# shortrate{1 to 23}
*
*  Save the OLS coefficients, residual sum of squares, X'X and
*  X'y
*
COMPUTE RSSOLS =%RSS
COMPUTE BETAOLS=%BETA
COMPUTE XXOLS  =INV(%XX)
COMPUTE XYOLS  =XXOLS*%BETA
SET OLS 1 %NREG = BETAOLS(T)
*
open plot gibbs1.rgf
GRAPH(NUMBER=0,HEADER='OLS Lag Distribution') 1
# OLS 2 %NREG
*
* Try a third degree PDL with the 25th lag constrained to zero
dec rect r
dim r(3,25)
ewise r(i,j)=(j-26)*j**(i-1)
encode r / enc1 enc2 enc3
# shortrate{0 to 24}
linreg(unravel) longrate
# constant enc1 enc2 enc3
compute betapdl = %beta
set pdl3 1 %nreg = betapdl(t)
open plot pdl.rgf
GRAPH(NUMBER=0,HEADER='Third Order PDL') 1
# pdl3 2 %nreg
* 
*
*  Set the number of draws
*  Set the prior degrees of freedom (NU) and mean (S2) for
*  sigma**2
*  Set the prior precision for the second difference in the
*  coefficients
*
COMPUTE NDRAWS=1000
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 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)
*
*  Our initial precision and coefficients come from OLS.
*  Compute the posterior degrees of freedom
*
COMPUTE HU  =1.0/%SEESQ
COMPUTE BETA=BETAOLS
COMPUTE SDOF=NU+%NOBS
*
DEC VECT U(%NREG)
*
*  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=1,UPPER=NDRAWS) 'Gibbs Sampler'
DO DRAW=1,NDRAWS
*
*   Draw precision conditional on previous BETA
*
   COMPUTE CENTER=NU*S2+RSSOLS+%QFORM(XXOLS,BETA-BETAOLS)
   COMPUTE HU    =2.0*%RANGAMMA(SDOF*.5)/CENTER
*
*     Draw beta conditional on precision of u. In general,
*     the formula for BETA would take the form:
*       COMPUTE BETA=VBETA*(HU*XYOLS+PRIORH*prior mean)+UBETA
*
*     In this example, the prior mean is 0 so the second
*     term in parentheses is omitted.
*
   COMPUTE HBETA=HU*XXOLS+PRIORH
   COMPUTE VBETA=INV(HBETA)
   COMPUTE UBETA=%DECOMP(VBETA)*(U=%RAN(1.0))
   COMPUTE BETA =VBETA*(HU*XYOLS)+UBETA
*
*     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
*
   INFOBOX(CURRENT=DRAW)
END DO DRAW
*
INFOBOX(ACTION=REMOVE)

DENSITY INTER 1 NDRAWS GINTER FINTER
DENSITY COEFF1 1 NDRAWS GCOEFF1 FCOEFF1
DENSITY SUMS 1 NDRAWS GSUMS FSUMS
*

SET COMOMENT1 1 %NREG = COMOMENT1/NDRAWS
SET COMOMENT2 1 %NREG = SQRT(COMOMENT2/NDRAWS-COMOMENT1**2)
SET UPPER = COMOMENT1+2.0*COMOMENT2
SET LOWER = COMOMENT1-2.0*COMOMENT2

open plot gibbs2.rgf
GRAPH(NUMBER=0,HEADER='OLS and Posterrior Mean Lag Distributions') 2
# COMOMENT1 2 %NREG
# OLS 2 %NREG

open plot gibbs3.rgf
GRAPH(NUMBER=0,HEADER='Posterior Mean Plus or Minus 2se') 3
# COMOMENT1 2 %NREG
# UPPER 2 %NREG 2
# LOWER 2 %NREG 2

open plot gibbs4.rgf
SCATTER(STYLE=LINES,HEADER='Posterior Density for Lag 0')
# GCOEFF1 FCOEFF1
open plot gibbs5.rgf
SCATTER(STYLE=LINES,HEADER='Posterior Density for Sum')
# GSUMS FSUMS