R中时间序列集成发生器的矢量化

Question

这个R脚本生成时间序列的集合。该级数由函数f(t) =α* f(t-0) + epsilon导出，其中epsilon是正态分布的随机数。

最后的结果是从alpha的不同值生成的集合列表。

如何将其矢量化？使用基本函数将是很棒的，但需要额外软件包的解决方案也是受欢迎的。

steps <- 1000 # number of times each "pseudo random walk" is iterated
N <- 5 # number of walks in each ensemble

mu <- 0 # normal distribution mean
mysd <- 1 # normal distribution standard deviation

alphas <- c(0, 0.5, 0.7, 1, 1.5, 2) # set of different alphas to generate ensembles with


# Pseudo random walk generator
generate.rw <- function(steps, alpha, my, mysd) {
  epsilons <- rnorm(n = steps, mean = mu, sd = mysd)
  rw <- vector(,steps)
  rw[1] <- epsilons[1]
  for (i in 2:steps) rw[i] <- alpha * rw[i-1] + epsilons[i]
  return(rw)
}

# Ensemble generator
ensemble <- function(N, steps, alpha, mu, mysd) {
  result <- matrix(,N,steps)
  for (i in 1:N) result[i,] <- generate.rw(steps, alpha, my, mysd)
  return(result)
}

# Get a list of ensembles with different values of alpha
ensembles <- lapply(alphas, ensemble, steps = steps, N = N, mu = mu, mysd = mysd)

Answer 1

你可以从使用

filter(rnorm(steps, mu, mysd), alpha, "recursive")

对于generate.rw和

replicate(generate.rw(steps, alpha, mu, mysd), n = N)

为了ensemble。顺便说一句，与1不同的α不是真正的随机游动；检查1阶的自回归过程(简称AR(1))和?arima.sim (filter的替代方案)。

Answer 2

你可能想看一下Wold定理。这样做的想法是，如果递归地求解AR(1)，就很容易将AR(1)矢量化。严格地说，Wold定理只在序列不遵循随机趋势(即α小于1)时才适用/有意义，但以后更多地说明这一点。

这是模型Yt = alpha Yt-1 +epsilon_t的递归解决方案：

Yt =和α^i* epsilon_t-i.

一个矢量化的解决方案现在是显而易见的。

res = rep(list(NA),length(alpha))

for (i in 1:length(alpha)){

  epsilon = rnorm(n = steps, mu, mysd)

  alpha_power = alpha[i]^seq(0,(steps-1))

  res[[i]] = alpha_power%*%epsilon

  #or if you want to save each Yt, alpha_power*epsilon 
}

上面的代码确实在alpha上循环。有一些方法可以避免这个循环，但是，考虑到相对较少的阿尔法，我觉得没有必要这样做。最重要的是，我将代码中最昂贵的部分(即需要多次迭代的部分)向量化。

我认为，这种方法比Julius的方法更快，因为它确实是矢量化的。我相信replicate是apply家族的一员，尽管我可能错了。

最后，当alpha >1时，您的模型实际上没有任何意义。当α<1时，如上文所示，冲击消失，最近的冲击被赋予最大的权重，例如.5^100*.5基本上为零。当α>1时，冲击的重量随着时间的增加而增加，例如2^100*5真的很大。换句话说，你的模型本质上是没有预测能力的，因为你的系列在未来的几步之后几乎可以在任何地方出现。

最后一件事，你应该，正如朱利叶斯所建议的，改变你问题的标题。AR( 1 )跟随随机游动当且仅当alpha =1。