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
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)
Hope that helps.