From: McGehee, Robert <Robert.McGehee_at_geodecapital.com>

Date: Fri 05 Aug 2005 - 03:06:40 EST

}

Qp.mean

}

set.seed(1)

plot(simCus5(Qp0=4.5))

R-help@stat.math.ethz.ch mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide!

http://www.R-project.org/posting-guide.html

R-help@stat.math.ethz.ch mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Fri Aug 05 03:12:33 2005

Date: Fri 05 Aug 2005 - 03:06:40 EST

Spencer,

On the first iteration of your simulation, all of the Qp.t + z.t < 0, so
you're adding a vector of rep(4.5, 20000) to a random distribution with
mean -0.1. So one would expect on iteration 2, your mean would have
dropped by about 0.1 (which it does). This process continues until about
the 20th iteration when we start seeing that a large number of our
initial starting points are floored at zero (because of the pmax). For
points greater than zero, we continue to subtract an average of 0.1
(actually less than this), but for those points already at zero, we're
actually adding a mean of 0.348 (since we can never subtract from a zero
number in this case), which starts the trajectory moving upward towards
its asymptote.

#That is, for those paths far above 0.1, we are subtracting

> mean(rnorm(10000, mean = -0.1))

[1] -0.1059246

#And for those paths already at zero, we are adding

> mean(pmax(0, rnorm(10000, mean = -0.1)))

[1] 0.3482376

To see a simulation a bit closer to what you were expecting, replace the starting values with a random distribution with mean Qp0.

i.e. replace

> Qp.t <- rep(Qp0, nSims)

with

> Qp.t <- rnorm(nSims, Qp0, sd = 3.7)

Robert

-----Original Message-----

From: Spencer Graves [mailto:spencer.graves@pdf.com]
Sent: Thursday, August 04, 2005 12:16 PM
To: r-help@stat.math.ethz.ch

Subject: [R] Counterintuitive Simulation Results

I wonder if someone can help me understand some
counterintuitive

simulation results. Below please find 12 lines of R code that
theoretically, to the best of my understanding, should produce
essentially a flat line with no discernable pattern. Instead, I see an
initial dramatic drop followed by a slow rise to an asymptote.

The simulation computes the mean of 20,000 simulated
trajectories of

400 observations each of a one-sided Cusum of independent normal
increments with mean EZ[t] = (-0.1) and unit variance. Started with any

initial value, the mean of the Cusum should approach an asymptote as the

number of observations increases; when started at that asymptote, it should theoretically stay flat, unlike what we see here.

I would think this could be an artifact of the simulation methodology, but I've gotten essentially this image with several independently programmed simulations in S-Plus 6.1, with all six different random number generators in R 1.9.1 and 2.1.1 and with MS Excel. For modest changes in EZ[t] < 0, I get a different asymptote but

pretty much the same image.

#################################################simCus5 <- function(mu=-0.1, Qp0=4.5, maxTime=400, nSims=20000){

Qp.mean <- rep(NA, maxTime)

Qp.t <- rep(Qp0, nSims)

for(i in 1:maxTime){

z.t <- (mu + rnorm(nSims)) Qp.t <- pmax(0, Qp.t+z.t) Qp.mean[i] <- mean(Qp.t)

}

Qp.mean

}

set.seed(1)

plot(simCus5(Qp0=4.5))

################################################# Thanks for your time in reading this. Best Wishes, Spencer Graves

Spencer Graves, PhD

Senior Development Engineer

PDF Solutions, Inc.

333 West San Carlos Street Suite 700

San Jose, CA 95110, USA

spencer.graves@pdf.com

www.pdf.com <http://www.pdf.com>

Tel: 408-938-4420

Fax: 408-280-7915

R-help@stat.math.ethz.ch mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide!

http://www.R-project.org/posting-guide.html

R-help@stat.math.ethz.ch mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html Received on Fri Aug 05 03:12:33 2005

*
This archive was generated by hypermail 2.1.8
: Fri 03 Mar 2006 - 03:39:40 EST
*