# Replicating simulation results from a paper

I’m reading R. Kreps' paper Parameter uncertainty in (log)normal distributions and trying to figure out how the simulations were done. In order to generate Figure 1, Eqn (2.41) was used. So this is what I did in R:

This is the equation they used: $Z_{eff} = v + z\sqrt{\frac{n(1+v^2)}{w}}$ where

Note: This code is for the case when $n=3$

v <- rt(10000000, 1)/sqrt(3-2)
w <- rchisq(10000000,2)
z <- rnorm(n=10000000, m=0, sd=1)
z_eff_3 <- v + z * sqrt((3*(1+v*v))/w)
plot(density(z_eff_3),ylim=c(0,1))


However, I got a completely different graph. I would truly appreciate it if someone can take a look and explain to me what I did wrong.

-
n a sample size and Zeff is defined in Eqn (2.41), but where are v, z, and w defined? Are you sure that z is a realisation from a N(1, 1) random variable? –  ocram Dec 20 '12 at 4:35

Note: In the original post, I noted that (aside from a minor error which has since been corrected) this was actually a problem with plotting the density, not with the simulation. At the time I wasn't sure what was going on, but I found that truncating the data seemed to solve the problem, so I recommended that. I have since figured out the real problem, so I decided to edit this to suggest a better solution.

With $n=3$, the distribution we are simulating from has very heavy tails. Even though the mean is zero and the inter-quartile range is about 4, in 100,000 draws (or more) we get a few on the order of 1,000,000 in absolute value.

By default, the density function estimates the density at 512 equally-spaced points that completely cover the range of the data (plus a little). In this case, we were only interested in plotting the density over a very small portion of that range--for $x \in [0,5]$. As a result, when we plotted that window, we got what looked like a flat horizontal line, which connected the estimates at points on the order of -2000 and 2000. Not very useful.

My original solution, truncating the data, gave decent looking plots, but wasn't really very satisfying. So I played around with it a bit and figured out the problem and how to solve it: the density function has arguments from and to, which give lower and upper bounds between which is estimates the density. By setting from=0, to=5 we get the following plot, without any truncation:

This is both easier and more correct (since the density is distorted slightly by the truncation) than my original solution. The complete code I used (adapted slightly from the code provided in the question) is:

v <- rt(100000, 1)/sqrt(3-2)
w <- rchisq(100000,2)
z <- rnorm(n=100000, m=0, sd=1)
z_eff_3 <- v + z * sqrt((3*(1+v*v))/w)

plot(density(z_eff_3,from=0,to=5),xlim=c(0,5))


Hope that helps.

-
Many thanks, Jonathan. Can you please share the R code? Thanks again! –  act00 Dec 20 '12 at 6:36
Here it is.. v <- rt(10000000, 1)/sqrt(3-2) w <- rchisq(10000000,2) z <- rnorm(n=10000000, m=0, sd=1) z_eff_3 <- v + z * sqrt((3*(1+v*v))/w) z_eff_3 <- z_eff_3[which((abs(z_eff_3))<10)] plot(density(z_eff_3),xlim=c(0,5),ylim=c(0,1)) –  act00 Dec 20 '12 at 9:30
I tried to produce the CDF in Fig 3 (for n=3), but again, I got a different graph. This is the R code: F_z_eff_3 <- ecdf(z_eff_3) plot(F_z_eff_3, xlim=c(5,10),ylim=c(0.85,0.99)) Can you please take a look and let me know what's wrong. Many thanks! –  act00 Dec 20 '12 at 20:38
If you're still working with the truncated data it's probably because we truncated pretty close in. You should be able to plot the ecdf of the original (untruncated) data, or, if you want to use the same data to plot the ecdf and the smoothed density estimate, truncating further away (at 100) worked for me. –  Jonathan Christensen Dec 20 '12 at 20:47
Thanks, JC! This works. –  act00 Dec 20 '12 at 22:22