Raw Message

Message-ID: <CABSrU1J4vJ3fZxF_n36mN7ORg8dYcJNaxUFJf6F2Q=Bk0gyKtg@mail.gmail.com>
Date: 2014-11-06T14:47:12Z
From: Matteo Richiardi
Subject: speed issue in simulating a stochastic process

I wish to simulate the following stochastic process, for i = 1...N
individuals and t=1...T periods:

y_{i,t} = y_0 + lambda Ey_{t-1} + epsilon_{i,t}

where Ey_{t-1} is the average of y over the N individuals computed at time
t-1.

My solution (below) works but is incredibly slow. Is there a faster but
still clear and readable alternative?

Thanks a lot. Matteo

rm(list=ls())
library(plyr)
y0 = 0
lambda = 0.1
N = 20
T = 100
m_e = 0
sd_e = 1

# construct the data frame and initialize y
D = data.frame(
  id = rep(1:N,T),
  t = rep(1:T, each = N),
  y = rep(y0,N*T)
)

# update y
for(t in 2:T){
  ybar.L1 = mean(D[D$t==t-1,"y"])
  for(i in 1:N){
    epsilon = rnorm(1,mean=m_e,sd=sd_e)
    D[D$id==i & D$t==t,]$y = lambda*y0+(1-lambda)*ybar.L1+epsilon
  }
}

ybar <- ddply(D,~t,summarise,mean=mean(y))

plot(ybar, col = "blue", type = "l")

	[[alternative HTML version deleted]]