Approximating the Distribution of the Product of Two Normally Distributed Random Variables

: The distribution of the product of two normally distributed random variables has been an open problem from the early years in the XXth century. First approaches tried to determinate the mathematical and statistical properties of the distribution of such a product using different types of functions. Recently, an improvement in computational techniques has performed new approaches for calculating related integrals by using numerical integration. Another approach is to adopt any other distribution to approximate the probability density function of this product. The skew-normal distribution is a generalization of the normal distribution which considers skewness making it ﬂexible. In this work, we approximate the distribution of the product of two normally distributed random variables using a type of skew-normal distribution. The inﬂuence of the parameters of the two normal distributions on the approximation is explored. When one of the normally distributed variables has an inverse coefﬁcient of variation greater than one, our approximation performs better than when both normally distributed variables have inverse coefﬁcients of variation less than one. A graphical analysis visually shows the superiority of our approach in relation to other approaches proposed in the literature on the topic.


Introduction
The first approach for finding the distribution of the product of two random variables with normal distribution was based on the product of two standard normal distributed random variables. In this case, Bessel type functions were used to approximate such a product [1]. Moments of the distributions of powers and products of normally distributed variables were studied in [2]. Other approaches presented in [3,4] considered a Pearson type function. In [5], there are several advances regarding the distribution of the product of two normally distributed random variables, but some of these advances need integrals that do not have an analytical solution and neither are defined at zero; see also [6]. A number of authors [7][8][9][10][11] employed the improvement of the computational methods to provide approximate numerical solutions for the mentioned integrals. Other alternative approaches are very limited. For example, the Bessel function method proposed in [9] presents a good approximation, but only for the product of two standard normal distributions and, moreover, it is not defined at zero. In [12][13][14][15], other approaches are presented, including Stein's operators and distributional theory. Another approach is to use any other distribution to approximate the probability density function (PDF) of the mentioned product. For example, in [16], the Gram-Charlier series are utilized to approximate a distribution while using the cumulative function distribution (CDF) of another

Product of Two Normal Random Variables
Let X ∼ N(µ X , σ X ) and Y ∼ N(µ Y , σ Y ) be two normally distributed random variables with correlation coefficient Corr(X, Y) = ρ ∈ [−1, 1]. Subsequently, the MGF of U = XY is given by Note that expression defined in (1) can be formulated in terms of the inverse coefficients of variation of X and Y established by δ X = µ X /σ X and δ Y = µ Y /σ Y , which are considered to be small when its value is less than one. Thus, from (1) and based on ρ, δ X and δ Y , we have that: • Skewness/kurtosis ratio: . (2) From (2), note that the skewness/kurtosis ratio of the product of two normally distributed random variables is composed by three factors. By assuming ρ ∈ [−1, 1] and δ X , δ Y ∈ [−1, 1], the ranges for these factors are analyzed as follows: Afterwards, the first factor ranges between [1, Thus, the second factor ranges between [− 16,16].
(iii) Third factor 3(1 + 2δ 2 ): its minimum value is three and reached at ρ = δ X = δ Y = 0, while its maximum value is 120 at ρ = 1, δ X = δ Y = −1. Hence, the denominator of this skewness/kurtosis ratio is in the range [3,120]. Therefore, one can deduce that the skewness/kurtosis ratio of the product of two normally distributed random variables is in the range [−0.38, 0.38].
In [28,29], the authors analyzed the coefficients of skewness and kurtosis of the product of two normally distributed random variables, according to different values of parameters of the two normal distributions.
The skewness/kurtosis ratio of the product of two normally distributed random variables is explored and depicted in Figure 1. From this figure, note that the influence of the correlation on the shape of the ratio is limited. However, the inverse coefficient of variation has an important effect on this ratio. There is a saddle-point at (0, 0) and, when the signs of δ X and δ Y are the same, the ratio tends to a maximum, whereas, for different signs, the ratio tends to a minimum. Figure 2 displays the scenario for small values of the inverse coefficient of variation, which is less than one in absolute value for the two normally distributed variables.

The Extended Skew-Normal Distribution
Now, let X be a random variable with SN distribution of location parameter ξ, scale parameter ω and skewness parameter α, which is denoted by X ∼ SN(ξ, ω, α); for more details of the SN distribution, see [25]. Subsequently, the PDF of X is given by where φ and Φ are the standard normal PDF and CDF, respectively. Note that, if α = 0, then the PDF of X defined in (3) reduces to the normal PDF. The SN distribution has interesting properties, some of which are used in this paper and presented next. If X ∼ SN(ξ, ω, α), then: As mentioned, the ESN distribution is an extension of the SN distribution, as indicated below; for more details of the ESN distribution, see [25]. In this case, the notation Y ∼ ESN(ξ, ω, α, τ) is employed, where ξ, ω, α and τ are location, scale, shape, and truncation parameters, respectively.
Note that, if τ = 0, the SN PDF defined in (3) is obtained and when τ = α = 0, the normal distribution is get. The standard ESN distribution is a particular case of the ESN distribution with parameters ξ = 0 and ω = 1. In addition, observe that the shape of the ESN PDF is very similar to the normal PDF (with different levels of skewness) when τ > 0. For negative values of τ, several types of shapes are obtained. When α is large (for example, for |α|> 50), the ESN PDF is concentrated around the mean, the variance is small, and the sign of the coefficient of skewness is determined by the sign of α. Statistics as the mean, variance, and coefficients of skewness and kurtosis of the ESN distribution can be calculated from its MGF after taking derivatives. Subsequently, following [25], the ESN MGF may be computed from the expression Hence, from (4), the following expressions are obtained: • Coefficient of skewness: • Excess of kurtosis: • Skewness/kurtosis ratio: Note that the ESN skewness/kurtosis ratio given in (5) only depends on the parameters α, τ and on the standard normal CDF, Φ namely. Therefore, there is no influence of the two other parameters ξ and ω on this ratio. The sign of α determines the sign of the skewness (positive or negative). If α is positive, the ratio moves from a positive value to a negative value, while for α < 0, the ratio goes from negative values to positive values. When α = 0, the ESN excess of kurtosis is zero and then its skewness/kurtosis ratio is not defined, but no problem exists for other values of α. When τ = 1, there is a discontinuity in this ratio. For τ > 1, the ratio is always negative and, for τ < 1, it is always positive. Notice that, with the values of skewness and kurtosis, it is possible to only calculate the values of two parameters (α, τ), in whose case a numerical solution to the equations is needed, which is sensitive to the initial point. With the values of the mean and variance, the other two parameters, ξ and ω namely, can be obtained.
Next, the ESN skewness/kurtosis ratio is graphically analyzed according to different values of its parameters; see Figure 3. In this figure, the ESN skewness/kurtosis ratio is depicted as a function of the parameters α and τ. The shape of this ratio is no smooth as in the case of the skewness/kurtosis ratio of the product of two normally distributed variables. The presence of gaps and changes on the direction of the shape of this ratio may be observed in this figure. The ratio is a non-continuous function with horizontal asymptotes at zero and ±0.4. Observe that the ESN skewness/kurtosis ratio is an increasing function of τ; see Figure 4. When τ is positive, the ESN skewness/kurtosis ratio is small.

Approximation for the Product Distribution with ESN Distributions
An approximation for the distribution of the product of two normally distributed variables using the ESN distribution is considered next. By using the MGF of the product of two normally distributed variables, the values of the mean, variance, coefficient of skewness, and excess of kurtosis are obtained. Subsequently, an ESN distribution with the same values for these statistics is assumed. By employing (4), the values of the parameters ξ, ω, α, and τ for the corresponding ESN distribution are calculated.
The parameters α and τ of the ESN distribution must be obtained in order to determine the skewness and kurtosis of the product U = XY as a function of δ X , δ Y and ρ. In a first stage, a system of two equations must be solved in terms of α and τ. Afterwards, in a second stage, to determine the ESN parameters ω and ξ, the other two equations for the mean and variance of the product must be solved. The solution must be reached using the values of parameters obtained in the previous stage. The following problems are associated with the first stage: • As non-linear equations are generated, suitable numerical methods are required.

•
As the solution obtained can be sensitive to the initial point, different points can be tested to avoid the stop at an inadequate point.

•
Values of the coefficient of skewness and the excess of kurtosis are not always reached.
However, our approximation has the following advantages: • A distribution is used considering a correlation ρ = 0 differently of other approaches [8,9], where only non-correlated variables (ρ = 0) were considered. • A known distribution is employed and its PDF for any value of the product can be calculated with no need of numerical integrals. • Inference using the ESN distribution is available.
Monte Carlo simulations are utilized with 10 9 replications of the product of two normally distributed variables while using a bivariate normal distribution with correlation ρ. Thus, the values of the mean, variance, skewness, and kurtosis are obtained according to the previously described procedure.
For this simulation, the R software was used [30]. Specifically, the sn package [31] was utilized, which provides functions for the SN and ESN PDFs and their cumulants. These cumulants are employed to derive the values of skewness and kurtosis.
The product of two normally distributed variables while using the bivariate normal distribution with correlation ρ was obtained. Its mean, variance, skewness and kurtosis were calculated with the corresponding MGF. The BB package [32] was used, which contains the BBSolve function, in order to solve the system of equations described above. Table 1 reports the values of the statistics of the product of two normally distributed variables X ∼ N(1.0, 2.0) and Y ∼ N(0.5, 2.0) with ρ = 0.5, where the performance of the approximation, based on the ESN distribution, can be analyzed. The values are the same for the mean and variance, while a slight difference for the excess of kurtosis is detected. However, in the case of the skewness, the difference is large, being around 30%. Figure 5 depicts the behavior of the approximation, where both inverse coefficients of variation are small (less than one), specifically, δ X = 0.5 < 1 and δ Y = 0.25 < 1. Now, consider the case of X ∼ N(1.0, 0.5) and Y ∼ N(0.5, 0.2), with ρ = 0.5, δ X = 2.0 > 1 and δ Y = 2.5 > 1, that is, the two variables have a large value for their inverse coefficients of variation. The results are reported in Table 1, with a high similarity between the theoretical values and those simulated by the ESN distribution; see also Figure 6.
The approximation that is proposed in this paper is sensitive to the values of the inverse coefficients of variation of the normally distributed variables. For the scenarios considered, when the two inverse coefficients of variation are less than one, the approximation based on the ESN distribution does not work well. However, when the two inverse coefficients of variation are greater than one, the approximation based on the ESN distribution works very good. Moreover, in order to have a good performance, there is only the need that one of the inverse coefficients of variation is greater than one, as reported in Figure 7 and Table 1. Therefore, the inverse coefficients of variation of both normally distributed variables determine the performance of our approximation based on the ESN distribution. Observe that the expression of the ESN skewness/kurtosis ratio, defined in (5), is always less than 0.4 (in absolute value). Thus, the ESN skewness/kurtosis ratio is never equal to the ratio of skewness over kurtosis of the product of two normally distributed variables when δ X , δ Y are in the range [−1, 1].
Note that other alternative approaches are very limited. For example, as mentioned, the approach that is based on the Bessel function method presents a good approximation, but only for the product of two standard normal distributions, and, moreover, it is not defined at zero. In Figures 5-7, other approaches like Pearson type III distribution [3], or the function presented in [7] that we call the Cui et al. approach, are depicted as comparison with our approach. Observe that, clearly, the results are worse than those obtained with our ESN approach. Notice that the function provided in [7] does not converge and then it is not possible to represent it in Figure 5. Table 1. Statistics for the product of two normally distributed variables: X and Y with ρ = 0.5 and δ X , δ Y as listed. Notice that our approximation that is based on the ESN distribution works well for several examples of the product of two normally distributed variables, when δ X or δ Y are greater than one; see Table 2 for positive and negative correlations. For both examples, the statistics of the two distributions are very similar and just slight differences are detected. Therefore, the PDF of U = XY can be well approximated using the ESN distribution with the calculated parameters appropriately. Table 2. Comparison of statistics of the product of two normality distributed variables X ∼ N(µ X , σ X ) and Y ∼ N(µ Y , σ Y ) with indicated ρ and ESN approximation statistics.

Conclusions and Future Research
This paper reported the following findings: (i) An approximation was proposed for the distribution of the product of two normally distributed random variables while using the extended skew-normal distribution. (ii) The moment generating function for the product of two normally distributed variables was considered as a function of the inverse coefficients of variation of both variables and of the corresponding correlation. From this moment generating function, the mean, variance, skewness, and kurtosis of this product were calculated. (iii) The probability density function of an extended skew-normal distribution is known and, using its cumulants, expressions for the mean, variance, skewness, and kurtosis were calculated. Skewness and kurtosis depended only on two parameters, while the mean and variance are functions of the four parameters of the extended skew-normal distribution. (iv) A numerical evaluation of the proposed approach was considered, which showed its potential for applications. (v) A graphical comparison of the proposed approach with other two approaches proposed in the recent literature on the topic was carried out showing the superiority of our approach.
In summary, the presented approximation required two stages. In the first stage, the values for the shape and truncation parameters of the extended skew-normal distribution related to the skewness and kurtosis values of the product of two normally distributed variables were obtained. The second stage consisted of determining the values of the other two parameters of location and scale, using the mean and variance of the product. The mentioned stages required the use of numerical methods to solve the corresponding equations, which revealed solutions that are sensitive to the initial point, and suggestions were given in order to avoid it. The approximation to the distribution of the product of two normally distributed variables using the extended skew-normal distribution worked well in most cases. The numerical evaluation reported the good performance of the proposed approach, suggesting it is a good choice when dealing with the distribution of the product of two normally distributed random variables. Therefore, this investigation may be an important knowledge addition to the tool-kit of diverse practitioners, including engineers, statisticians, and data scientists.
Some open problems that arose from the present investigation are the following: (i) The approximation to the distribution of the product of two normally distributed variables, while using the extended skew-normal distribution, showed some problems when the inverse coefficients of variation both are in the range [−1, 1]. For this range, the adequate values of skewness and kurtosis cannot be reached and, thus, the approximation does not work well there. Some further research should be conducted in this line in order to improve our approach. (ii) Other distributions may be considered for approximating the distribution of the product of two normally distributed random variables; see, for example, [20,21,24,[33][34][35]. (iii) Applications of the new methodology proposed in this investigation can be of interest in diverse areas, where the distribution of the product of two normally distributed random variables is required; see, for example, [28,36]. (iv) If there are two random samples, the first one from the distribution of the product of two normally distributed random variables, and the second one from the fitted extended skew-normal distribution, and then some non-parametric tests can be used to compare them. Moreover, a sophisticated approach, as, for example, bootstrapping, can be applied in this context. Therefore, the proposed methodology in this investigation promotes new challenges and offers an open door for exploring other theoretical and numerical issues. Research on these and other issues are in progress and their findings will be reported in future articles.