CODEWITHSUNDEEP
rtime-seriesvectorizationrandom-walk
HAVB
HAVB

Reputation: 1888

Vectorizing a time series ensemble generator in R

This R script generates ensembles of time series. The series are derived from function f(t) = alpha * f(t-0) + epsilon, where epsilon is a random number from a normal distribution.

The final result is a list of ensembles generated from different values of alpha.

How can it be vectorized? Using base functions would be great, but solutions that require additional packages are welcome too.

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)

Upvotes: 2

Views: 155

Answers (2)

Jacob H
Jacob H

Reputation: 4513

You might want to (re)look at the Wold Theorem. The idea is that if you recursively solve your AR(1) it will be easy to vectorize. Strictly speaking the Wold Theorem only applies/makes sense when the series do not follow a stochastic trend (i.e. the alpha is less than 1), but more about that later.

This is the recursive solution to the model Yt = alpha Yt-1 + epsilon_t:

Yt = Sum alpha^i * epsilon_t-i.

A vectorized solution is now obvious.

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 
}

The above code does loop over the alphas. There are ways to avoid this loop, however, given the relative small number of alphas I feel like there is no need to. The most important thing is that I vectorized the most expensive part of your code (i.e. the part which required numerous iterations)

This approach will be faster than Julius' approach, I think, because it truly is vectorized. I believe that replicate is part of the apply family, though I could be wrong.

Finally, when alpha > 1, your model does not really make any sense. When alpha < 1, as you can see above in the equation the shocks die off and the nearest shocks are given the most weight, e.g. .5^100*.5 is essentially zero. When alpha > 1 the weight on a shock increases over time, e.g. 2^100*.5 is really really big. Said another way, your model essentially has no predictive power as your series can be pretty much anywhere after only a few steps into the future.

Last thing, you should, as Julius suggests, change the title of your question. An AR(1) follows a random walk if and only if alpha = 1.

Upvotes: 0

Julius Vainora
Julius Vainora

Reputation: 48241

You may start with using

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

for generate.rw and 

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

for ensemble. By the way, with alpha different from 1 it is not really a random walk; check autoregressive processes of order 1 (AR(1) for short) and ?arima.sim (alternative to filter).

Upvotes: 2

Related Questions