Векторизация генератора ансамбля временных рядов в R

Question

Этот скрипт R генерирует ансамбли временных рядов. Серия получается из функции f (t) = alpha * f (t-0) + epsilon, где epsilon является случайным числом из нормального распределения.

Конечным результатом является список ансамблей, созданных из разных значений альфа.

Как это можно проставить в векторе? Использование базовых функций было бы здорово, но решения, требующие дополнительных пакетов, тоже приветствуются.

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

Возможно, вы захотите (пере) рассмотреть теорему Вольда. Идея состоит в том, что если вы рекурсивно решаете свой AR (1), это будет легко прорисовать. Строго говоря, теорема Вольда применима/имеет смысл, когда ряд не следует стохастическому тренду (т. Е. Альфа меньше 1), но об этом позже.

Это рекурсивное решение модели Yt = альфа Yt-1 + epsilon_t:

Yt = Sum альфа^я * epsilon_t-я.

Теперь векторизованное решение становится очевидным.

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 
} 

Вышеприведенный код делает петлю над альфа. Однако есть способы избежать этого цикла, учитывая относительное небольшое количество альфов, которые, как мне кажется, не нужны. Самое главное, что я векторизовал самую дорогую часть вашего кода (т. Е. Часть, которая требовала многочисленных итераций)

Этот подход будет быстрее, чем подход Юлиуса, я думаю, потому что он действительно векторизован. Я считаю, что replicate является частью семьи apply, хотя я могу ошибаться.

Наконец, когда alpha> 1, ваша модель действительно не имеет никакого смысла. Когда альфа < 1, как вы можете видеть выше в уравнении, удары отмирают, и наибольший вес получают самые близкие удары, например. .5^100 * .5 по существу равна нулю. Когда альфа> 1, вес при ударе возрастает с течением времени, например. 2^100 * .5 действительно очень большой. С другой стороны, ваша модель, по сути, не имеет прогностической силы, так как ваша серия может быть практически любой после нескольких шагов в будущее.

Последнее, что вам следует, как предлагает Юлий, изменить название вашего вопроса. AR (1) следует случайному блужданию тогда и только тогда, когда альфа = 1.