Simulación de GARCH en R

Question

Estoy haciendo una simulación de un modelo GARCH. El modelo en sí no es demasiado relevante, lo que me gustaría preguntarle es sobre la optimización de la simulación en R. Más que nada si ve algún espacio para la vectorización, lo he pensado pero no puedo verlo. Hasta ahora lo que tengo es la siguiente:

Let:

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

Así que este es el código:

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)) 
    } 

Quiero hacer una simulación de los rendimientos 5 periodos a partir de ahora, por lo que se desarrollará este digamos 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)) 

Creo que esto está funcionando razonablemente rápido, pero me gustaría preguntarle si hay alguna manera de apr abordar este problema de una mejor manera.

Gracias!

Answer 1

En lugar de usar números en su ciclo, puede usar vectores de tamaño N: que eliminen el bucle oculto en 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) 

Answer 2

edificio sobre la respuesta de Vicente, lo único que cambió fue dfining zt todos a la vez y cambiar apply(ret, 1, sum) a rowSums(ret) y se aceleró un poco. He intentado tanto compilado, pero sin diff importante:

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 

probable que todavía margen de mejora :-)