1

我必须对 AR(1) 时间序列采用蒙特卡罗方法。我必须生成 10,000 个长度为 100 的时间序列，然后我必须为每个时间序列获取第一步自相关 rho_1。我的问题是我只获得了自相关的 NA 值，并且计算需要很长时间。我对计算 AR(1) 时间序列没有任何问题。谢谢您的帮助 ：）

gen_ar <- function(a,b,length,start)
{
 z<-rep(0,length)                                 
 e<-rnorm(n=length,sd=1)  
 z[1]<-start
  for (i in 2:length)
  {
    z[i]<-a+b*z[i-1]+e[i]
  }
 z
}

mc <- matrix(c(rep(0,10000000)),nrow=10000)
 for (i in 1:10000)
 {
  mc[i,] <- gen_ar(0.99,1,100,0)
 }

 ac <- matrix(c(rep(0,10000)),nrow=1)
  for (i in 1:10000){
   for (j in 1:99){
    ac[i] <- cor(mc[i,j],mc[i,j+1])
   }
  }
4

1 回答 1

0

除了统计数据，我认为这可以实现您的目标，而且我没有得到 NA。我改变了它的完成方式，因为你说它进展缓慢。

mc <- matrix(rep(NA,1E5), nrow=100)
for(i in seq_len(100)){
    mc[,i] <- arima.sim(model=list(ar=0.99), n=100, sd=1) + 1
}

myAR <- function(x){
    cor(x[-1], x[-length(x)])
}

answer <- apply(mc, 2, myAR)

我跳过了最后一组嵌套的 for 循环并将它们替换为apply(). 它似乎更容易阅读，并且可能更快。此外，为了使用apply()，我创建了一个名为 的函数myAR，它执行与循环cor()中相同的计算。for()

现在，我做了一些统计调整。这些主要是在模拟步骤中。

首先，您模拟的 AR(1) 过程的系数等于 1，这对我来说似乎很奇怪（这不会是固定的，而且 arima.sim() 甚至不会让您模拟这种类型的过程）。

此外，您的“a”参数在每个时间步将时间序列加 1。换句话说，您的时间序列从 1 单调增加到 100，因为系数等于 1。这也会使您的时间序列非平稳，并且具有如此强的正斜率，该cor()函数可能会返回1为估计的相关性，无论模拟 AR 系数的值。我假设您希望长期均值悬停在 1 附近，因此 1 在模拟后简单地添加到整个时间序列中，而不是在每个时间步迭代。

假设您确实想通过在每个时间步添加一些常数 (a) 来生成非平稳时间序列，您可以执行以下操作：

myInnov <- function(N=100, a=1, SD=1) {a + rnorm(n=N, sd=SD)}
mc2 <- matrix(rep(NA,1E7), nrow=100)
for(i in seq_len(1E5)){
    mc2[,i] <- arima.sim(model=list(ar=0.99), n=100, innov=myInnov(a=1, N=100, SD=1)) + 1
}

我希望这个对你有用。

于 2013-06-18T23:05:35.367 回答