*
* GARCHIMPORT.RPF
* Example of importance sampling, applied to MC integration on an GARCH model.
*
* RATS User's Guide, Example from Section 16.6.
*
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