MONTEARCH.RPF

MONTEARCH.RPF shows Monte Carlo examination of test statistic, in this case, a test for ARCH (Engle, 1982). The model simulated is

\begin{equation} y_t = \beta _0 + \beta _1 X_t + u_t \label{eq:boot_montearchmean} \end{equation}

\begin{equation} u_t \sim N\left[ {0,\sigma ^2 \left( {1 + \alpha {\kern 1pt} u_{t - 1}^2 } \right)} \right] \label{eq:boot_montearchvariance} \end{equation}

Because the (conditional) variance of \(u\) is not constant, the draws for the model must be generated using the non-linear systems solution procedure. SIMULATE only draws from a fixed distribution, so we have to add a third equation to the model which computes a standard Normal variate V. FRMLS for \(u\) and \(y\) are then “identities.”

This program does 10000 draws, generating the process over the period 2 through 150, starting with U(1) = 0. Since this starts “cold” from an arbitrary value, we use only the last 100 observations to compute the test statistics. The 50 discarded observations are known as the burn-in, a phrase is used often in simulations of time-dependent data.

compute ndraws=10000

compute burnobs=50

compute useobs=100

compute endobs=burnobs+useobs

The specifics of the simulated model are that the mean process \eqref{eq:boot_montearchmean} is \(2+3 X_t\) where \(X\) is a time trend, and the \(\alpha\) in \eqref{eq:boot_montearchvariance} is .2 (which is a rather modest value). As this is a "power" calculation, the choice of a value relatively close to the null makes sense. It's not clear how sensitive the results will be to the choice of the mean of the process since the efficiency of estimates of the trend may be affected by the ARCH effect in the residuals.

This sets up the X variable and the model. Note that unlike the choice of the mean function, the choice of SIGMA should be irrelevant, since all it does is to scale the generated data and a scale multiple which will wash out of the test statistic.

set x 1 endobs = t

compute sigma=1.0,alpha =.20

frml(variance=sigma^2) vdef v = 0.0

frml(identity) udef u = v*sqrt(1+alpha*u{1}^2)

frml(identity) ydef y = 2.0 + 3.0*x + u

group archmod vdef udef ydef>>y

You might think that you could simplify this three-equation setup by using a %RAN function directly in the UDEF FRML. This, however, will draw a new random number at every pass through the solver, so the model will never converge. SIMULATE draws only a single set of values and solves the model given those.

The test is the LM form, which uses the \(TR^2\) from a regression of the squared residuals on lagged squared residuals. This will have an asymptotic \(\chi^2\) distribution with 1 degree of freedom. These are saved in entries of the series TESTSTAT, which runs from 1 to NDRAWS.

Note that we have to set the "pre-sample" value of U to 0. If we don't, the UDEF applied at entry 2 will have a missing value for U{1}, which will produce a data series of all NA's. Note also that we need to give the explicit range on the first LINREG, since we don't want to use the generated Y series until entry BURNOBS+1.

set u 1 1 = 0.0

set teststat 1 ndraws = 0.0

do draw=1,ndraws

simulate(model=archmod,from=2,to=endobs)

linreg(noprint) y burnobs+1 endobs

# constant x

set usquared = %resids^2

linreg(noprint) usquared

# constant usquared{1}

compute teststat(draw)=%trsquared

end do draw

To evaluate how well the test works on this, we compare the test statistics with the .05 and .01 critical values of the test statistic under the null. This is done with SSTATS with the MEAN option. (SSTATS on a logical expression counts the number of times that the logical expression is true; SSTATS with MEAN will convert that to a probability).

compute crit05=%chisqr(.05,1),crit01=%chisqr(.01,1)

sstats(mean) 1 ndraws (teststat>crit05)>>reject05 (teststat>crit01)>>reject01

?"P(reject at .05)" #.### reject05

?"P(reject at .01)" #.### reject01

Full Program