[R] calculating recursive sequences
Vaidotas Zemlys
mpiktas at gmail.com
Thu Dec 29 13:33:32 CET 2005
Hi,
I was trying to repeat the estimation of threshold GARCH models from
the book "Analysis of Financial Time Series" by Ruey S. Tsay, and I
was succesfull, but I had to use "for" loop, which is quite slow. The
loop is necessary, since you need to calculate recursive sequence. Is
there a faster way to do this in R, without using loops?
The model is such:
r_t = \mu + \alpha_2 r_{t-2} + a_t
a_t = \sigma_t\varepsilon_t
\sigma_t^2 =
\beta_1a_{t-1}^2+\beta_2\sigma_{t-1}^2+
1_{\{a_{t-1}>0\}}(\gamma_0+
\gamma_1a_{t-1}^2+\gamma_2\sigma^2_{t-1})
It is asummed that \varepsilon_t are iid and normal with zero mean and
variance one. The data given is r_t, and you have to estimate
variables, \mu, \alpha, \beta and \gamma. Since
\varepsilon_t=\frac{a_t}{\sqrt{sigma_t}}
using the equations we calculate a_t and \sigma_t and estimate the
variables using maximum likelihood method. a_t can be estimated
directly using first equation and rt. \sigma_t^2 depends on
sigma_{t-1}^2, so it should be calculated recursively.
The function calculating negative log-likelihood of this problem I wrote:
garchln <- function(p,rt) {
n <- length(rt)
at <- rt[4:n]-p[1]-p[2]*rt[4:n-2]
u <- as.numeric(at>0)
h <- rep(0,length(at))
# h is \sigma_t^2
for(i in 1:(length(h)-1)) {
h[i+1] <- p[3]*at[i]^2+p[4]*h[i]+u[i]*(p[5]+p[6]*at[i]^2+p[7]*h[i])
}
#Maximum likelihood function
sum(log(h[-1])+(at[-1]^2)/h[-1])/2
#list(h=h[-1],at=at[-1])
}
For fitting I used optim, with methods "Nelder-Mead" and "BFGS",
Initial parameter values from the book are
0.03 -0.03 0.10 0.60 0.10 0.05 0.10
The fitted values from the book are
0.043 -0.022 0.098 0.954 0.060 -0.052 -0.069.
The link to the data used:
http://www.gsb.uchicago.edu/fac/ruey.tsay/teaching/fts/d-ibmln99.dat
For this problem recursive sequence is linear, so it is possible to
calculate it as a linear equations solution, but it is easy to think
of the case where the recursion is non-linear. Is the speed-up
possible only by writing C or Fortran code with loops?
Vaidotas Zemlys
--
Doctorate student, http://www.mif.vu.lt/katedros/eka/katedra/zemlys.php
Vilnius University
More information about the R-help
mailing list