Matteo, Ah — OK, N=20, I did not catch that. You have nested for loops, which R is known to be exceedingly slow at handling — if you can reorganize the code to eliminate the loops, your performance will increase significantly.
Tom On Thu, Nov 6, 2014 at 7:47 AM, Matteo Richiardi <matteo.richia...@gmail.com > wrote: > 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. > [[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.