Simulation of GARCH in R

Question

I am doing a simulation of a GARCH model. The model itself is not too relevant, what I would like to ask you is about optimizing the simulation in R. More than anything if you see any room for vectorization, I have thought about it but I cannot see it. So far what I have is this:

Let:

#    ht=cond.variance in t
#    zt= random number 
#    et = error term
#    ret= return
#    Horizon= n periods ahead

So this is the code:

randhelp= function(horizon=horizon){
    ret <- zt <- et <- rep(NA,horizon)#initialize ret and zt et
    for( j in 1:horizon){
      zt[j]= rnorm(1,0,1)
      et[j] = zt[j]*sqrt(ht[j])
      ret[j]=mu + et[j]

      ht[j+1]= omega+ alpha1*et[j]^2 + beta1*ht[j]
    }
    return(sum(ret))
  }

I want to do a simulation of the returns 5 periods from now, so I will run this let's say 10000.

#initial values of the simulation
ndraws=10000
horizon=5 #5 periods ahead
ht=rep(NA,horizon) #initialize ht
ht[1] = 0.0002
alpha1=0.027
beta1 =0.963
mu=0.001
omega=0


sumret=sapply(1:ndraws,function(x) randhelp(horizon))

I think this is running reasonably fast but I would like to ask you if there is any way of approaching this problem in a better way.

Thanks!

Gallaway

Looks like mu and omega are not defined. Can you move zt outside the loop and generate all the random values at once, then index into them? Have you tried the library(compiler)? – Salvia

library(compiler); f1 <- cmpfun(randhelp) is all it takes to give it a whirl. Sometimes it gives a big boost, othertimes not so much...but easy to test so worth a short IMHO. Good luck :) – Salvia

Answer 1

Instead of using numbers in your loop, you can use vectors of size N: that removes the loop hidden in sapply.

randhelp <- function(
  horizon=5, N=1e4, 
  h0 = 2e-4, 
  mu = 0, omega=0,
  alpha1 = 0.027,
  beta1  = 0.963
){
  ret <- zt <- et <- ht <- matrix(NA, nc=horizon, nr=N)
  ht[,1] <- h0
  for(j in 1:horizon){
    zt[,j]   <- rnorm(N,0,1)
    et[,j]   <- zt[,j]*sqrt(ht[,j])
    ret[,j]  <- mu + et[,j]
    if( j < horizon )
      ht[,j+1] <- omega+ alpha1*et[,j]^2 + beta1*ht[,j]
  }
  apply(ret, 1, sum)
}
x <- randhelp(N=1e5)

Bluing

Answer 2

building on Vincent's response, all I changed was dfining zt all at once and switching apply(ret, 1, sum) to rowSums(ret) and it sped up quite a bit. I tried both compiled, but no major diff:

randhelp2 <- function(horizon = 5, N = 1e4, h0 = 2e-4, 
                       mu = 0, omega = 0, alpha1 = 0.027, 
                       beta1  = 0.963 ){
    ret <- et <- ht <- matrix(NA, nc = horizon, nr = N)
    zt  <- matrix(rnorm(N * horizon, 0, 1), nc = horizon)
    ht[, 1] <- h0
    for(j in 1:horizon){
        et[, j]  <- zt[, j] * sqrt(ht[, j])
        ret[,j]  <- mu + et[, j]
        if( j < horizon )
            ht[, j + 1] <- omega + alpha1 * et[, j] ^ 2 + beta1 * ht[, j]
    }
    rowSums(ret)
}

system.time(replicate(10,randhelp(N=1e5)))
   user  system elapsed 
  7.413   0.044   7.468 

system.time(replicate(10,randhelp2(N=1e5)))   
   user  system elapsed 
  2.096   0.012   2.112 

likely still room for improvement :-)

Electrolyze