Replicating simulation results from a paper

Question

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.

Answer 1

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.