GARCHのシミュレーションR

Question

私はGARCHモデルのシミュレーションを行っています。モデル自体はあまり関連性がありません。私はRのシミュレーションを最適化することについてお聞きしたいと思います。ベクトル化の余地があれば、何と考えても見えません。これまでのところ私が持っているものは以下です：

はしてみましょう：私は今から戻って5つの周期のシミュレーションをしたい

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

ので、私：

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

は、これはコードですこれを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)) 

私はこれがかなり速く実行されていると思いますが、私はあなたに何か方法があるかどうか質問したいと思いますこの問題をより良い方法で解き明かします。

ありがとうございます！

Answer 1

ループ内の数字を使用する代わりに、sapplyに隠されたループを削除するサイズNのベクトルを使用できます（ ）。ヴィンセントの応答に

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

建物、私が変更されたすべてを一度ztをdfiningとrowSums(ret)にapply(ret, 1, sum)を切り替え、それはかなり高速化しました。私は両方のコンパイルを試してみましたが、重大な差分：

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 

はおそらくまだ