R中GARCH模拟

我正在模拟一个GARCH模型。模型本身并不太相关，我想问你的是如何优化R中的模拟。如果你看到有任何向量化空间，那么最重要的是我已经考虑过了，但我看不到它。到目前为止，我已经是这样的：R中GARCH模拟

令：

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

所以这是代码：

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

我想从现在起5期做回报的模拟，所以我将运行这个让我们说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))

我认为这是运行速度相当快，但我想问你是否有任何方式appr以更好的方式解决这个问题。

谢谢！

来源

2012-04-02 aatrujillob

貌似'mu'和'没有定义omega'。你可以在循环外面移动'zt'并且一次生成所有的随机值，然后索引它们吗？你尝试过'库（编译器）'吗？ – Chase 2012-04-02 02:03:51

'library（compiler）; f1 < - cmpfun（randhelp）'是所有它需要给它一个旋转。有时它会有很大的提升，其他时间不会太多...但容易测试，所以值得一试IMHO。祝你好运：） – Chase 2012-04-02 02:26:06

除了在循环中使用数字之外，还可以使用大小为N的矢量：删除隐藏在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)

来源

2012-04-02 02:18:39

建筑，都让我改变了一次性dfining zt和开关apply(ret, 1, sum)到rowSums(ret)，它加快了不少。我都尝试编译，但是没有大的差异：

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

可能仍有提升空间:-)

来源

2012-04-02 02:44:21

回答

相关问题