If random variables $X$, $Y$ are dependent, and $X>0$ $Y>0$, we have $P(XY\leq z)=P(X\leq z/Y)=\int_{0}^{\infty} \int_{0}^{z/y} f_{XY}(x,y) dx dy$.
How about the case of $n$ dependent random variables? i.e., how to use the joint PDF $f_{X_1,X_2,\ldots,X_n}(x_1,x_2,\ldots,x_n)$ to obtain the CDF $P(X_1 X_2\cdots X_n\leq z)$?
I have an idea, but not sure if it is right.
$P(X_1 X_2\cdots X_n\leq z)\\=E(P(X_1 \leq \frac{z}{X_2\cdots X_n}\mid X_2=x_2,\ldots,X_n=x_n))\\=\int_0^\infty \cdots \int_0^\infty dx_2\cdots dx_n \int_0^{\frac{z}{x_2\cdots x_n}}f_{X_1,X_2,\ldots,X_n}(x_1,x_2,\ldots,x_n)dx_1$.
Let $X_1,....,X_n$ be random variables with joint density function $f(x_1,...,x_n)$.
Consider the transformation $g(x_1,...,x_n)=(x_1...x_n,x_2,...,x_n)$. The inverse transformation is $g^{-1}(z_1,...,z_n)=(\frac{z_1}{z_2...z_n},z_2,...,z_n)$ and it follows that the jacobian of the inverse transformation is $|J|=\frac{1}{z_2...z_n}$.
Then by the transformation theorem the density function of the vector $(z_1,...,z_n)$ is $f_Z(z_1,...,z_n)=f(\frac{z_1}{z_2...z_n},z_2,...,z_n)|\frac{1}{z_2...z_n}|$ and the marginal density function of $Z_1$ can be written as:
$f_{Z_1}(z_1)=\int_0^\infty...\int_0^\infty f(\frac{z_1}{z_2...z_n},z_2,...,z_n)|\frac{1}{z_2...z_n}| dz_2...dz_n$
Note that the solution above can be obtained by applying the Leibniz rule to the solution you proposed. Take $F(\frac{z}{x_2...x_n},...,x_n)-F(0,x_2,...,x_n)=\int_0^\frac{z}{x_2,...,x_n} f_{X_1,...,X_n}(x_1,...,x_n) dx_1$ and try to calculate $\frac{d}{dz}\int_0^\infty...\int_0^\infty F(\frac{z}{x_2...x_n},...,x_n)-F(0,x_2,...,x_n) dx_2...dx_n.$