Oops, I forgot the 'cumsum' stuff. Here it is again, hopefully working this time.
> require(MASS) > fz <- function(n, t, rho){ + f <- array(dim = c(t, 2, n)) + V <- matrix(c(1, rho, rho, 1), ncol = 2) + for(i in 1 : n){ + f[,, i] <- apply(mvrnorm(n = t, mu = c(0,0), Sigma = V), + 2, cumsum) + } + return(f) + } > > set.seed(123) > out <- fz(1000, 20, rho = 0.7) > ts.plot(out[,,500]) # one realization of the two random walks > apply(out, 1, function(x) cor(t(x))[1,2]) # cor at times 1,...,20 [1] 0.7306115 0.7229866 0.7103006 0.6970370 0.7016218 [6] 0.6968639 0.6933828 0.7039743 0.7006013 0.6969842 [11] 0.6935726 0.6853286 0.6809578 0.6754606 0.6804075 [16] 0.6881239 0.6896976 0.6917677 0.6929290 0.6888701 Best, Giovanni On Tue, 2010-05-11 at 17:29 -0400, Sergio Andrés Estay Cabrera wrote: > Dear Giovanni, > > Thanks so much for your answer, but your script returns gaussian white > noise not a random walk, at least the time series generated don't have > the expected periodogram for a random walk. That's the reason why I use > cumsum, the sum of a white noise is a easy way produce a random walk > which periodogram scales at 1/f. > > finally I solve my problem with a brute-force algorithm. Searching in > the web I found some pages about random walks and neural networks where > people said that this task is very difficult and is in a exploratory > step nowadays. > > Thank you very much for your help > > Sergio > > > > > Giovanni Petris wrote: > > Hi Sergio, > > > > Your function does not estimate what you want. In fact it does not > > estimate anything useful. A random walk is not stationary; in > > particular, the variance at time t is t. Therefore, estimating variances > > based on one run, averaging over time, does not make any sense. This is > > what you are doing with the command cor(cumsum(s[,1]),cumsum(s[,2])) > > (the correlation is obtained by standardizing the covariance, dividing > > by the square root of the variances). > > > > On the other hand, if you replicate your simulation (as you do) and > > compute the correlation of the simulated random walks at any _fixed_ > > time, you can see that the result is what you expect. Here it is. > > > > > >> require(MASS) > >> fz <- function(n, t, rho){ > >> > > + f <- array(dim = c(t, 2, n)) > > + V <- matrix(c(1, rho, rho, 1), ncol = 2) > > + for(i in 1 : n){ > > + f[,, i] <- mvrnorm(n = t, mu = c(0,0), Sigma = V) > > + } > > + return(f) > > + } > > > >> > >> set.seed(123) > >> out <- fz(1000, 20, rho = 0.7) > >> apply(out, 1, function(x) cor(t(x))[1,2]) # cor at times 1,...,20 > >> > > [1] 0.7306115 0.6810394 0.6879712 0.7095613 0.7325612 0.6922090 > > 0.7007823 > > [8] 0.6828095 0.7144861 0.6932104 0.7061211 0.6781331 0.6898726 > > 0.6926878 > > [15] 0.6949656 0.6995899 0.7524582 0.6912938 0.6931928 0.6820243 > > > > Best, > > Giovanni > > > > On Mon, 2010-05-10 at 14:55 -0400, Sergio Andrés Estay Cabrera wrote: > > > >> Dear R users and specially Albyn and Giovanni, > >> > >> thanks for your answers, but in fact I supposed the same at the > >> beginning of my problem. However, when I generate the data seldom I > >> obtain the expected correlation. For example using this code: > >> > >> fz<-function(n,t,rho){ > >> f<-NULL > >> for(i in 1:n){ > >> s<-rmvnorm(n=t,mean=c(0,0),sigma=matrix(c(1,rho,rho,1),ncol=2),method='svd') > >> paso<-cor(cumsum(s[,1]),cumsum(s[,2])) > >> f<-c(f,paso)} > >> f<-f > >> } > >> > >> and then plot the histogram of the results, it is possible to observe > >> that the distribution of the values is asymmetric with most of the > >> simulations close to 1 when the value of rho is higher than 0.3 and > >> looks like a uniform distribution with values below 0.3. > >> > >> I suspect than the only possibility is using a brute force algorithm. > >> > >> > >> Any advice would be helpful > >> > >> Sergio A. Estay > >> *CASEB * > >> Departamento de Ecología > >> Universidad Catolica de Chile > >> > >> > >> > >> > >> Albyn Jones wrote: > >> > >>> Sums of correlated increments have the same correlation as the original > >>> variables... > >>> > >>> library(mvtnorm) > >>> X<- matrix(0,nrow=1000,ncol=2) > >>> for(i in 1:1000){ > >>> Y <- rmvnorm(1000,mean=mu,sigma=S) > >>> X[i,] <- apply(Y,2,sum) > >>> } > >>> cor(Y) > >>> [,1] [,2] > >>> [1,] 1.0000000 0.4909281 > >>> [2,] 0.4909281 1.0000000 > >>> > >>> So, unless you meant that you want the _sample_ correlation to be > >>> pre-specified, you are all set. > >>> > >>> albyn > >>> > >>> On Sun, May 09, 2010 at 09:20:25PM -0400, Sergio Andrés Estay Cabrera > >>> wrote: > >>> > >>> > >>>> Hi everybody, > >>>> > >>>> > >>>> I am trying to generate two random walks with an specific correlation, > >>>> for > >>>> example, two random walks of 200 time steps with a correlation 0.7. > >>>> > >>>> I built the random walks with: > >>>> > >>>> x<-cumsum(rnorm(200, mean=0,sd=1)) > >>>> y<-cumsum(rnorm(200, mean=0,sd=1)) > >>>> > >>>> but I don't know how to fix the correlation between them. > >>>> > >>>> With white noise is easy to fix the correlation using the function > >>>> rmvnorm > >>>> in the package mvtnorm > >>>> > >>>> I surfed in the web in the searchable mail archives in the R web site > >>>> but > >>>> no references appears. > >>>> > >>>> If you have some advices to solve this problems I would be very thankful. > >>>> > >>>> Thanks in advance. > >>>> > >>>> Sergio A. Estay > >>>> *CASEB * > >>>> Departamento de Ecología > >>>> Universidad Catolica de Chile > >>>> > >>>> -- > >>>> “La disciplina no tiene ningún mérito en circunstancias ideales. ” – > >>>> Habor Mallow > >>>> > >>>> ______________________________________________ > >>>> 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. > >>>> > >>>> > >>>> > >>> > >>> > >> > > > > > > > > > > > > > > ______________________________________________ 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.