* * 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