$X$ is $\mathcal N(0,4)$, $Y$ is $\mathcal N(0,5)$, $Z = X + Y$
I need to simulate 1000 values for each of these variables, $X$,$Y$,$Z$.
I have simulated 1000 values for both $X$ and 1000 values for $Y$.
When simulating 1000 values for $Z$, should I use the values already simulated for $X$ and $Y$?
Or should I simulate new values for $X$ and $Y$?
I feel like the other answers so far have not been completely clear on why you must reuse the same samples that you obtained for $X$ and $Y$: this is necessary to obtain a sample of $(X,Y,Z)$ which has the correct joint distribution.
If $X$ and $Y$ are independent (that wasn't explicit in the question), then: $$\text{Cov}(X,Z) = \text{Cov}(X,X+Y)=\text{Cov}(X,X)=\text{Var}(X)=4$$
We then have $\text{Cor}(X,Z) = 2/3$, so $X$ and $Z$ are definitely not independent. Then, intuitively, the method you use to sample $Z$ must use your existing samples for $X$ and $Y$ in one way or another.
If we compare both approaches (using the same $X$ and $Y$, or new $X$ and $Y$):
set.seed(123)
x <- rnorm(1000, 0, sqrt(4))
y <- rnorm(1000, 0, sqrt(5))
z <- x + y
x_new <- rnorm(1000, 0, sqrt(4))
y_new <- rnorm(1000, 0, sqrt(5))
z_new <- x_new + y_new
par(mfrow=c(1,2))
plot(x,z, main = paste0("Same sample: sample correlation = ", format(cor(x,z),digits=3)))
plot(x,z_new, main = paste0("New sample: sample correlation = ", format(cor(x,z_new),digits=3)))
par(mfrow=c(1,1))
We get the following plots:
On the left, reusing the same sample of $X$ and $Y$ yields a correlation between $X$ and $Z$ which is roughly what we expected (2/3). It's not shown here but it also has the correct joint distribution of $Y$ and $Z$, and the full $(X,Y,Z)$.
On the right, the correlation is roughly zero. Using new samples of $X$ and $Y$ completely destroys the dependence structure with $Z$.
If $X \sim \mathsf{Norm}(\mu=0,\sigma=2)$ and $Y \sim \mathsf{Norm}(\mu=0,\sigma=\sqrt{5}),$ then $Z = X+Y \sim \mathsf{Norm}(\mu=0,\sigma=\sqrt{4+5} = 3).$
If you simulate many values of $X,$ then the a histogram of those values imitates the density function of $\mathsf{Norm}(\mu=0,\sigma=2).$ Also, the sample mean of the many values $X_i$ will be $\bar X \approx 0$ and the their standard deviation will be $S_X = \frac{1}{n-1}\sum_{i=1}^n (X_i - \bar X)^2 \approx 2.$ Using R for the simulating, summarizing, and graphing, we have:
set.seed(2020)
x = rnorm(10^5, 0, 2
summary(x); sd(x)
Min. 1st Qu. Median Mean 3rd Qu. Max.
-8.359695 -1.349914 -0.014586 -0.005547 1.339879 8.318011
[1] 1.994003 # aprx sigma.x = 2
hist(x, prob=T, col="skyblue2",
main="Histogram of Sample from NORM(0, 2) with Density")
curve(dnorm(x, 0, 2), add=T, col="red")
Similarly, for $Y,$ but without the graph, we get:
summary(y); sd(y)
Min. 1st Qu. Median Mean 3rd Qu. Max.
-9.607853 -1.505329 -0.006414 0.001022 1.500837 10.620498
[1] 2.235846 # aprx sigma.y = sqrt(5) = 2.236
Finally, if you add simulated vectors of values x and y, you will simulate the distribution of
$$Z = X+Y \sim
\mathsf{Norm}(\mu=0,\sigma=\sqrt{4+5} = 3),$$
Thus, a data summary of z will show $\bar Z \approx 0$ and $S_Y \approx 3.$
z = x + y
summary(z); sd(z)
Min. 1st Qu. Median Mean 3rd Qu. Max.
-13.299918 -2.027446 -0.003853 -0.004525 2.010361 12.703835
[1] 2.994405 # aprx sigma.z = 3
hist(z, prob=T, col="skyblue2", ylim=c(0,.13),
main="Histogram of Sample from NORM(0, 3) with Density")
curve(dnorm(x, 0, 3), add=T, col="red")
Notes: In case you are interested in details of R, here are a few.
(1) Normal probability functions in R, (such as rnorm to simulate and density dnorm, use the standard deviation $\sigma$ instead of the variance $\sigma^2.$
(2) If you are going to superimpose a density function on a histogram, then the histogram needs to be a 'density histogram'. This requires parameter prob=T in hist.
(3) The R procedure curve to plot a function requires the argument of the function always to be written as x--- regardless of the context.
(4) The parameter ylim=c(0,.13) of the procedure hist makes the plotting window tall enough to
accommodate the density curve without cutting off
its mode.
(5) When you write x + y in R for two vectors of the same length (here length $10^5 = 100,000)$ then R takes the
element-wise sum, yielding a vector of the same length.
The table below shows the first six out of 100,000 such additions.
head(cbind(x, y, z))
x y z
[1,] 0.7539442 1.0037551 1.7576993
[2,] 0.6030967 0.2978469 0.9009436
[3,] -2.1960463 3.9858896 1.7898433
[4,] -2.2608118 -3.0263022 -5.2871140
[5,] -5.5930686 3.8786439 -1.7144247
[6,] 1.4411470 0.3790924 1.8202394
If $Z=X+Y$, and you have $X$ and $Y$, then just sum them, that's the simulation for $Z$.
X, one forY, one forZ. As a guess, you probably don't want the samples to be independent, but how could we actually know that? You never said what you were trying to do. $\endgroup$