*
* Example 12.1 from pp 556-558
*
calendar(w) 1962:1:5
all 1999:9:10
open data wgs1yr.dat
data(format=free,org=columns) / year t1
open data wgs3yr.dat
data(format=free,org=columns) / year t3
*
set dt1 = t1-t1{1}
set dt3 = t3-t3{1}
*
* ML estimation of the REG-AR2 model
*
boxjenk(reg,ar=2,maxl) dt3
# constant dt1
*
* Initial values for parameters
* Step one: Linear regression for the beta parameters
*
linreg dt3
# constant dt1
compute b0=%beta(1),b1=%beta(2)
*
* Step two: AR(2) regression on the residuals to get the phi
* and sigmasq values
*
linreg %resids
# %resids{1 2}
compute rstart=%regstart(),rend=%regend()
compute phi1=%beta(1),phi2=%beta(2)
compute sigmasq=%sigmasq
*
* Prior precision (inverse of covariance matrix) for the beta's
*
compute [symm] hbeta0=%diag(||.25,.25||)
*
* Prior precision for the phi's
*
compute [symm] hphi0=%diag(||1.0/.25,1.0/.16||)
*
* Prior for the sigma squared. pscale is the "prior mean" for sigma squared
* (actually the reciprocal of mean of the density for 1/sigma-squared) and pdf is
* the prior degrees of freedom for the chi-square.
*
compute pscale=.1,pdf=10.0
*
dec symm vbeta vphi
dec vect betam betadraw phim phidraw
*
compute ndraws=2100
dec vect[series] stats(5)
do i=1,5
   set stats(i) 1 ndraws = 0.0
end do i
labels stats
# "b0" "b1" "phi1" "phi2" "sigma"
*
do draws=1,2100
   *
   * Draw beta's given phi's and sigma's.
   *
   * Step one - filter the data by the autoregressive polynomial, and compute
   * the cross product matrix of the filtered data.
   *
   set dt3f = dt3-phi1*dt3{1}-phi2*dt3{2}
   set dt1f = dt1-phi1*dt1{1}-phi2*dt1{2}
   set cf   = 1-phi1-phi2
   cmom(noprint)
   # cf dt1f dt3f
   *
   *  Step two - compute the inverse of the sums of the precision matrices. The
   *  precision for the data is sigma**-2 X'X (inverse of the standard OLS
   *  covariance matrix).
   *
   compute vbeta=inv(hbeta0+%xsubmat(%cmom,1,2,1,2)/sigmasq)
   *
   *  Step three - compute the mean of the posterior. Ordinarily, this would be
   *  vbeta*(hbeta0*beta0+%xsubmat...), but since the prior mean is zero, we
   *  can skip the first term.
   *
   compute betam=vbeta*%xsubmat(%cmom,1,2,3,3)/sigmasq
   *
   *  Make a random draw from the Normal distribution for the beta's.
   *
   compute betadraw=betam+%ranmvnormal(%decomp(vbeta))
   compute beta1=betadraw(1),beta2=betadraw(2)
   *
   * Draw phi's given beta's and sigma's. This is basically the same process as above
   *
   set %resids = dt3-beta1-beta2*dt1
   cmom(noprint)
   # %resids{1 2 0}
   compute vphi=inv(hphi0+%xsubmat(%cmom,1,2,1,2)/sigmasq)
   compute phim=vphi*%xsubmat(%cmom,1,2,3,3)/sigmasq
   compute phidraw=phim+%ranmvnormal(%decomp(vphi))
   compute phi1=phidraw(1),phi2=phidraw(2)
   *
   * Draw sigma's given beta's and phi's. Compute the sum of squares of the residuals,
   * and draw from the appropriate inverse chi-squared density.
   *
   set u = %resids-phi1*%resids{1}-phi2*%resids{2}
   sstats rstart rend u**2>>%rss
   compute sigmasq=(pscale*pdf+%rss)/%ranchisqr(%nobs+pdf)
   *
   * Save the statistics of interest
   *
   compute stats(1)(draws)=beta1
   compute stats(2)(draws)=beta2
   compute stats(3)(draws)=phi1
   compute stats(4)(draws)=phi2
   compute stats(5)(draws)=sqrt(sigmasq)
end do draws
*
report(action=define)
report(atrow=1,atcol=1) "Parameter" "b0" "b1" "phi1" "phi2" "sigma"
report(atrow=2,atcol=1) "Mean"
report(atrow=3,atcol=1) "St. Error"
do i=1,5
   stats(noprint) stats(i) 101 2100
   report(atrow=2,atcol=i+1) %mean
   report(atrow=3,atcol=i+1) sqrt(%variance)
   density(type=hist,maxgrid=20) stats(i) 101 2100 xx dx
   scatter(style=bar,vmin=0.0,hlabel=%l(stats(i)))
   # xx dx
end do i
report(action=format,picture="*.###")
report(action=show)