Simulazione di GARCH in R

Question

Sto facendo una simulazione di un modello GARCH. Il modello in sé non è troppo rilevante, quello che vorrei chiederti riguarda l'ottimizzazione della simulazione in R. Più di ogni cosa, se vedi qualsiasi spazio per la vettorizzazione, ci ho pensato ma non riesco a vederlo. Finora quello che ho è questa:

Let:

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

Quindi questo è il codice:

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

voglio fare una simulazione dei rendimenti 5 periodi d'ora, così ho verrà eseguito questo diciamo 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)) 

Penso che questo è in esecuzione ragionevolmente veloce, ma vorrei chiedere a voi se non v'è alcun modo di appr oaching questo problema in un modo migliore.

Grazie!

Answer 1

Invece di utilizzare i numeri nel loop, è possibile utilizzare i vettori della dimensione N: che rimuove il loop nascosto 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) 

Answer 2

edificio sulla risposta di Vincent, tutto quello che è stato cambiato dfining zt tutto in una volta e il passaggio ad apply(ret, 1, sum)rowSums(ret) e accelerato un po '. Ho provato sia compilato, ma nessun grande 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 

camera probabilmente ancora di miglioramento :-)