[R] Random walk
Giovanni Petris
gpetris at uark.edu
Tue May 11 22:56:21 CEST 2010
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 at 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.
> >>
> >>
> >
> >
> >
>
>
More information about the R-help
mailing list