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]] ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.