Product of Gaussian with a 4-level discrete random variable

Question

Similar to this question, but a small generalization of it, Imagine that I have $X \sim N(0,\sigma^2)$ and

\begin{align*} P(Y=y) = \begin{cases} \frac{1}{4} & y = -a_2\\ \frac{1}{4} & y = -a_1\\ \frac{1}{4} & y = +a_1\\ \frac{1}{4} & y = +a_2\\ \end{cases} \end{align*} where $a_2>a_1>0$. What is the distribution of $Z = XY$ ?

I followed similar steps like

\begin{align} P(Z<z) &=P(XY <z)\\ &=\sum_i \Big( P(a_iX<z|Y=a_i)P(Y=a_i)+P(-a_iX<z|Y=a_i)P(Y=a_i) \Big) \\ &=\frac{1}{4} \sum_i \Big( P(a_iX<z)+P(a_iX > -z) \Big) \\ &=\frac{2}{4} \sum_i \Big( P(a_iX<z) \Big) \\ \end{align} where the last step is due to symmetry. I know that $ P(a_iX<z)$ is $N(0,(a_i\sigma)^2)$. But what is the sum ? I thought it may be $N(0,(a_1\sigma)^2+(a_2\sigma)^2)$ but I think it is wrong. Thanks a lot

Answer 1

Most of your argument looks reasonable, though I would say $$P(Z \le z) = \frac12 P\left(X \le \frac{z}{a_1}\right) +\frac12 P\left(X \le \frac{z}{a_2}\right) = \frac12\Phi\left(\frac{z}{a_1\sigma} \right)+ \frac12\Phi\left(\frac{z}{a_2\sigma} \right)$$ as the conclusion, using $\Phi$ as the CDF of a normal distribution.

The variance of $Z$ is $\frac{\sigma^2}{2}(a_1^2+a_2^2)$, half of what you suggested.

$Z$ is a mixture distribution of two normals with different variances so does not have a normal distribution itself. Here is an illustration using R, with $\sigma=3,a_1=2,a_2=10$, where a cumulative distribution function estimated from a sample is shown in black, the mixture cdf in red (a close fit taking into account simulation noise) and a normal distribution with the same variance in blue (not so close a fit), followed by the same for the densities.

set.seed(2023)
cases <- 10^5
sigma <- 3
a1 <- 2
a2 <- 10
X <- rnorm(cases, mean=0, sd=sigma)
Y <- sample(c(-a2,-a1,a1,a2), cases, replace=TRUE)
Z <- X*Y
var(Z)
# 470.2348
1/2*((a1*sigma)^2+(a2*sigma)^2)
# 468

plot.ecdf(Z)
curve(1/2*(pnorm(x,0,a1*sigma)+pnorm(x,0,a2*sigma)),
   from=min(Z),to=max(Z), add=TRUE, col="red")
curve(pnorm(x,0,sqrt(1/2*((a1*sigma)^2+(a2*sigma)^2))),
   from=min(Z),to=max(Z), add=TRUE, col="blue")

plot(density(Z))
curve(1/2*(dnorm(x,0,a1*sigma)+dnorm(x,0,a2*sigma)),
   from=min(Z),to=max(Z), add=TRUE, col="red") 
curve(dnorm(x,0,sqrt(1/2*((a1*sigma)^2+(a2*sigma)^2))),
   from=min(Z),to=max(Z), add=TRUE, col="blue")