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

Question

此 R 脚本生成时间序列的集合。该系列源自函数 f(t) = alpha * 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。顺便说一下，如果 alpha 不同于 1，那么它并不是真正的随机游走；检查 1 阶自回归过程（简称 AR(1)）和 ?arima.sim（替代 filter）。

Answer 2

您可能想（重新）看一下 Wold 定理。这个想法是，如果你递归地解决你的 AR(1) 它会很容易矢量化。严格来说，当序列不遵循随机趋势（即 alpha 小于 1）时，Wold 定理仅 applies/makes 有意义，但稍后会详细介绍。

这是模型 Yt = alpha Yt-1 + epsilon_t:

的递归解

Yt = Sum alpha^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。有一些方法可以避免这种循环，但是，考虑到相对较少的 alpha，我觉得没有必要这样做。最重要的是我向量化了你代码中最昂贵的部分（即需要多次迭代的部分）

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

最后，当 alpha > 1 时，您的模型实际上没有任何意义。当 alpha < 1 时，正如您在上面的等式中看到的那样，冲击消失并且最近的冲击被赋予最大的权重，例如.5^100*.5 本质上是零。当 alpha > 1 时，冲击的权重会随着时间的推移而增加，例如2^100*.5 真的很大。换句话说，您的模型基本上没有预测能力，因为您的系列在未来仅需几步之后就可以出现在几乎任何地方。

最后，您应该按照 Julius 的建议更改问题的标题。当且仅当 alpha = 1 时，AR(1) 遵循随机游走。

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

Vectorizing a time series ensemble generator in R

r

time-series

vectorization

random-walk