Problem with distribution of $Y=X_1 X_2$ via known $f_{X_1,X_2}$ plus extension to arbitrary monomial function $Y=g(X_1, ... X_n)$

Question

I hope you're well.

Question 1

The first problem I have is essentially that I want to find the product distribution (CDF) of two random variables ($Y=g(X_1,X_2)=X_1 X_2$), who's joint PDF distribution, $f_{X_1,X_2}(x_1,x_2)$, is known (may be dependent/independent) and where everything is continuous.

I aim only to implement the solution numerically in a program, and wish to hopefully extend this to solve for the quotient distribution ($Y=g(X_1,X_2)=\frac{X_1}{X_2}$), and also for any monomial function of any number of random variables within a known joint PDF (e.g. $Y=g(X_1 X_2,X_3)=\frac{X_1^{3/2} X_2}{X_3}$ where only $f_{X_1,X_2,X_3}$ is known initially).

Hence, the formula which I am using is a very general one, so that I might eventually extend it: $$\int_{A}f_Y(y')\mathrm{d}y'=\int_{g^{-1}(A)}f_{\pmb{\vec{X}}}(\pmb{\vec{x}})\mathrm{d}\pmb{\vec{x}}$$ $$\implies F_Y(y)=\int_{g^{-1}((-\infty,y])}f_{\pmb{\vec{X}}}(\pmb{\vec{x}})\mathrm{d}\pmb{\vec{x}}.$$ Hence, for $Y=X_1 X_2$, $$F_Y(y)=\iint_{g^{-1}(x_1 x_2 \leq y)}f_{X_1,X_2}(x_1,x_2)\mathrm{d}x_1 \mathrm{d}x_2.$$

I have implemented this formulation into Octave, and it's worked wonders for when $X_1 \sim \mathcal{N}(10,1)$ and $X_2 \sim \mathcal{U}(1,3)$, and where both are independent. However, when I test for two independent standard normal distributions ($X_1=X_2 \sim \mathcal{N}(0,1)$) then I get problems. The reason I have deduced is that there is a hiccup with the region that the preimage, $g^{-1}(x_1 x_2 \leq y)$, defines. In the below image is the contour map of the joint PDF of the two std normal distributions just for visualisation purposes:

In Octave and Python (see below in code appendix), I integrate simply by using a mask function for the region of integration, and then do trapezoidal numerical integration over that region to get the 3-D "volume" as it were of the region. It is the mask function that defines the preimage $g^{-1}(x_1 x_2 \leq y)$, and so for example, when $y=1$, I get this plot of the mask function:

As one may expect logically, via the definition of the preimage, the 2nd and 4th quadrant all satisfies the region of the mask function when $y=1$. However, that means that when $y\leq 0$ there is a positive value for the CDF, which shouldn't be the case. Knowing that the product of two independent standard normal distributions is equal to a Chi squared distribution with one degree of freedom, and comparing, I get very different answers because of the preimage region including all of the 2nd and 4th quadrants, even when $y=0$, where $F_y(0)$ should equal $0$, but does not (it roughly is 0.5, since all of the 2nd and 4th quadrants cover half of the whole 2-D region).

If I were to look at the quotient distribution, the same issue would appear, since the region of integration would automatically include the whole of the 2nd and 4th quadrant when $y=0$. Even if I were to just ignore the 2nd and 4th quadrants in the product and quotient cases, I would only get $F_Y(\infty)=0.5$ which of course would be wrong.

And so, my question is how might I get around this problem to retrieve the correct region of integration for this specific case of the product distribution of $Y=X_1 X_2$ where $X_1=X_2 \sim \mathcal{N}(0,1)$, both being independent? (and possibly also solve the same issue with the quotient distribution $Y=X_1/X_2$)

This problem is basically for any product distribution where the joint PDF of the initially known two distributions extends into the 2nd and 4th quadrant (e.g. $Y = X_1 X_2$ where $X_1=X_2 \sim \mathcal{U}(-0.5,0.5)$)

Question 2

My final question is simply to ask, regarding my derivation, is there a way to find the CDF of $Y$, from initially only knowing the joint CDF of $X_1,...X_n$ (instead of the joint PDF)?

Hence, could there be some general formula like my initial one, where $$F_Y(y)=\mathcal{H}(F_{\pmb{\vec{X}}}(\pmb{\vec{x}}))?$$

And so for example, what I mean is that in my previous initial formula using $f_{\pmb{\vec{X}}}$, $\mathcal{H}$ was the whole integral formula. And so is there some formulation using only the CDF for what is initially only known?

The reason being that I image it may be easier to work with an initial joint CDF.

Any help is greatly appreciated, and sorry for this quite long question, although I hope it might be an interesting problem.

EDIT: I created a python script also (see below in appendix) and so these plots may also help visualisation:

---

P.S. I am a civil engineering student, and so my level of fundamental mathematical understanding for the problem I have may be lacking, and so apologies in advance if that is the case.

Code Appendix

Octave Code for the product distribution of $Y=X_1 X_2$ where $X_1=X_2 \sim \mathcal{N}(0,1)$ (can basically be copied and pasted into MATLAB also without any major errors):

clc;
clear all;
close all;
pkg load statistics

function pd = f_X1_X2(x1,x2)
  pd = normpdf(x1,0,1) .* normpdf(x2,0,1);
end

x1 = [-10:0.01:10]';
x2 = [-10:0.01:10]';

[X1,X2] = meshgrid(x1,x2);

f_X1_X2_PDF = f_X1_X2(X1,X2);

xs = {x1,x2};
Xs = {X1,X2};

function cp = F_Y(arrayx,arrayX,f,y)
  M = zeros(size(arrayX{1}));
  sel = (arrayX{1} .* arrayX{2} <= y);
  M(sel) = 1;

  fM = f_X1_X2(arrayX{1},arrayX{2}) .* M;
  I1 = trapz(arrayx{1},fM,2);
  I2 = trapz(arrayx{2},I1,1);

  cp = I2;
end

lb = 0;
dx = 0.1;
ub = 5;
F_Y_CDF = arrayfun(@(y) F_Y(xs,Xs,@f_X1_X2,y),[lb:dx:ub]);

figure(1);
plot([lb:dx:ub],F_Y_CDF,"k--","linewidth",2);
hold on;
plot([lb:dx:ub],chi2cdf([lb:dx:ub],1),"k-","linewidth",2);
legend("Derived Y=X1*X2","chi2cdf(z,1)");

## F_Y_CDF_diff = diff(F_Y_CDF);
## f_Y_pdf = F_Y_CDF_diff ./ dx;
## figure(2);
## plot([lb:dx:ub](2:end),f_Y_pdf);

EDIT: Python Code:

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.ticker import LinearLocator
from scipy import integrate
from scipy.stats import norm, uniform, chi2
import seaborn as sns
sns.set()

plt.rc('text', usetex=True)
plt.rc('font', family='serif')

# define joint pdf of standard norm variables
def f_X1_X2(x1,x2):
    return (norm.pdf(x1,0,1) * norm.pdf(x2,0,1))

# make arrays of x1 and x2 values
x1 = np.arange(-10,10.01,0.01)
x2 = np.arange(-10,10.01,0.01)

X1, X2 = np.meshgrid(x1,x2)

# make matrix of joint pdf values
f_X1_X2_PDF = f_X1_X2(X1,X2)

xs = [x1,x2]
Xs = [X1,X2]

# make function for CDF of product distribution of X1 and X2
# utlising 2D integration via trapezoid function
def F_Y(xs,Xs,f_x,y):
    # mask function
    M = np.zeros(np.shape(Xs[0]))
    sel = (Xs[0] * Xs[1] <= y)
    M[sel] = 1

    # integrate
    fM = f_X1_X2(Xs[0],Xs[1]) * M
    I1 = integrate.trapezoid(fM,xs[0],1)
    I2 = integrate.trapezoid(I1,xs[1],0)
    
    return I2

# make arrays necessary for plotting CDF
lb = 0
dx = 0.1
ub = 5
F_Y_CDF = np.array(list(map(lambda y : F_Y(xs,Xs,f_X1_X2,y),list(np.arange(lb,(ub+dx),dx)))))

# plot 3D PDF
fig, ax = plt.subplots(subplot_kw={"projection": "3d"})
ax.plot_surface(X1,X2,f_X1_X2_PDF,cmap=cm.coolwarm,linewidth=0,antialiased=False)
plt.title(r'Plot of $f_{X_1,X_2}$')

# plot derived CDF
fig1, ax1 = plt.subplots()
ax1.plot(np.arange(lb,(ub+dx),dx),F_Y_CDF,'k--',linewidth=2,label=r'Derived $Y=X_1 X_2$')
ax1.plot(np.arange(lb,(ub+dx),dx),chi2.cdf(np.arange(lb,(ub+dx),dx),1),'k-',linewidth=2,label=r'$X^2_1$ (what $Y=X_1 X_2$ should be)')
plt.title('Derived CDF and Analytical Comparison')
ax1.legend()

plt.show()