Distribution of a MA(1) process

Question

Suppose I have this MA(1) model:

$y_t = \mu + \epsilon_t + \theta \epsilon_{t-1}$ with $\epsilon_t \sim \mathcal{N}(0,\sigma^2)$

The marginal distribution of $y_t$ for all $t$ is $\mathcal{N}(\mu,\sigma^2(1 + \theta^2))$ but the classical theory says that since the $y_t$ variables are not independent, the likelihood of the sample is not the product of the marginal densities. One has to use conditional densities instead (given $\epsilon_0 = 0$):

$y_1 \sim \mathcal{N}(\mu,\sigma^2)$
$y_2 \sim \mathcal{N}(\mu + \theta \epsilon_1,\sigma^2)$
$...$
$y_t \sim \mathcal{N}(\mu + \theta \epsilon_{t-1},\sigma^2)$

And the likelihood is the product of these conditional densities, right?

Now, I wrote a small program in R simulating a MA(1) process and displaying the distribution of the generated $y_t$ variables:

n = 100000
mu = 5
theta = 3
sigma = 2

y = rep(1, n)
eps.1 = 0

for (t in 1:n)
{
    eps = rnorm(1, mean = 0, sd = sigma)
    y[t] = mu + eps + theta * eps.1
    eps.1 = eps
}

hist(y, probability = TRUE)
curve(dnorm(x, mean = mu, sd = sigma * sqrt(1 + theta^2)), col = "blue", add = TRUE)

And I find out that the joint distribution of $y_t$ is $\mathcal{N}(\mu,\sigma^2(1 + \theta^2))$. How come?