*
* GARCHIMPORT.RPF
* RATS Version 8, User's Guide, Example 16.3
* Importance sampling for MC integration on a GARCH model.
*
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 draw=1,ndraws
   *
   * Draw coefficients from a multivariate t
   *
   compute betau=xbase+%ranmvt(fxx,drawdf)
   *
   * 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 %ranmvt 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
   compute bxx=bxx+weight*%outerxx(betau)
   *
   *  Save the drawn value for persist and weight.
   *
   compute persist(draw)=betau(4)+betau(5)
   compute weights(draw)=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,smoothing=1.5) persist 1 ndraws xx dx
scatter(style=line,header="Density of Persistence Measure (alpha+beta)")
# xx dx