*
* GARCHIMPORT.PRG
* Manual Example 13.3
*
open data haversample.rat
calendar(m) 1957
data(format=rats) 1957:1 2006:12 ftbs3
*
graph(header="US 3-Month Treasury Bill Rate")
# ftbs3
*
* Estimate a linear regression to get its residual variance for use as the
* pre-sample value.
*
linreg ftbs3
# constant ftbs3{1}
compute h0=%sigmasq
*
* Estimate a GARCH with Normally distributed errors
*
garch(p=1,q=1,reg,presample=h0,distrib=normal) / ftbs3
# constant ftbs3{1}
compute normallogl=%logl
*
* Set up for importance sampling.
*  FXX is the factor of the covariance matrix of coefficients
*  XBASE is the estimated coefficient vector
*  FBASE is the final log likelihood
*
compute fxx  =%decomp(%xx)
compute xbase=%beta
compute fbase=%logl
*
* This is the number of draws and the degrees of freedom of the multivariate t
* distribution from which we're drawing.
*
compute ndraws=10000
compute drawdf=5.0
*
* Initialize vectors for the first and second moments of the coefficients
*
compute sumwt=0.0,sumwt2=0.0
compute [vect] b=%zeros(%nreg,1)
compute [symm] bxx=%zeros(%nreg,%nreg)
dec vect u(%nreg) betau(%nreg)
*
* This is used for an estimated density of the persistence measure (alpha+beta).
* We need to keep the values of the draws and the observation weights in order to
* use the density instruction.
*
set persist 1 ndraws = 0.0
set weights 1 ndraws = 0.0
*
do draws=1,ndraws
   *
   *  Draw a %NREG vector from a t with drawdf degrees of freedom, mean 0 and
   *  identity scale matrix.
   *
   compute u=%rant(drawdf)
   *
   *  Premultiply by fxx to correct the scale matrix and add the mean to get
   *  the draw for the coefficients
   *
   compute betau=xbase+fxx*u
   *
   *  Compute the density of the draw relative to the density at the mode. The
   *  log of this is returned by the value of %ranlogkernel() from the %rant function.
   *
   compute gx=%ranlogkernel()
   *
   *  Compute the log likelihood at the drawn value of the coefficients
   *
   garch(p=1,q=1,reg,presample=h0,method=eval,initial=betau) / ftbs3
   # constant ftbs3{1}
   *
   *  Compute the difference between the log likelihood at the draw and the
   *  log likelihood at the mode.
   *
   compute fx=%logl-fbase
   *
   *  It's quite possible for %logl to be NA in the GARCH model (if the GARCH parameters
   *  go sufficiently explosive). If it is, make the weight zero. Otherwise exp(...) the
   *  difference in the log densities.
   *
   compute weight=%if(%valid(%logl),exp(fx-gx),0.0)
   *
   *  Update the accumulators. These are weighted sums.
   *
   compute sumwt=sumwt+weight,sumwt2=sumwt2+weight**2
   compute b=b+weight*betau,bxx=bxx+weight*%outerxx(betau)
   *
   *  Save the drawn value for persist and weight.
   *
   compute persist(draws)=betau(4)+betau(5)
   compute weights(draws)=weight
end do draws
*
* The efficacy of importance sampling depends upon function being estimated, but
* the following is a simple estimate of the number of effective draws.
*
disp "Effective Sample Size" sumwt**2/sumwt2
*
* Transform the accumulated first and second moments into mean and variance
*
compute b=b/sumwt,bxx=bxx/sumwt-%outerxx(b)
*
* Create a table with the Monte Carlo means and standard deviations
*
report(action=define)
report(atrow=1,atcol=1,fillby=cols) "a" "b" "gamma" "alpha" "beta"
report(atrow=1,atcol=2,fillby=cols) b
report(atrow=1,atcol=3,fillby=cols) %sqrt(%xdiag(bxx))
report(action=show)
*
* Estimate the density of the persistence measure and graph it
*
density(weights=weights,smpl=weights>0.00,bandwidth=.05) persist 1 draws xx dx
scatter(style=line,header="Density of Persistence Measure (alpha+beta)")
# xx dx