Finding the distribution of $Z$, which is the product of two independent Pareto distributed random variables Consider two independent random variables, $X$ and $Y$. Both follow a Pareto distribution with minimum value $1$ and shape parameters $\alpha_X$ and $\alpha_Y$, respectively. Hence, the PDFs are given by
\begin{align*}
f_X(x) &= \alpha_X ~ x^{-\alpha_X-1} \\
f_Y(y) &= \alpha_Y ~ y^{-\alpha_Y-1}.
\end{align*}
I need to find the distribution of $Z$, which is the product of $X$ and $Y$: $Z \equiv XY$.

What I've tried: The distribution of $Z$ seems to be a product distribution and according to the linked page the PDF of $Z$ should be
\begin{align*}
f_{Z}(z)=\int _{{-\infty }}^{{\infty }}f_{X}\left(x\right)f_{Y}\left(z/x\right){\frac  {1}{|x|}}\,dx.
\end{align*}
However, trying to use this formula turned out to be problematic due to the following reasons (or due to my inexperience):
First, integrating from $-\infty$ to $\infty$ seems problematic if my distribution has a minimum value of $1$. I therefore replaced the lower bound of the integral by $1$ (should I?). Besides, if the lower bound is $1$, then using the absolute value of $x$ seems superfluous. Therefore, after plugging in my PDFs from above and integrating, I get:
\begin{align*}
f_Z(z) = \alpha_X \alpha_Y z^{-\alpha_Y-1} \left[ \frac{1}{-\alpha_X+\alpha_Y} x^{-\alpha_X+\alpha_Y} \right]_1^\infty.
\end{align*}
This raises the second issue: I expected the terms where I plug in $\infty$ to vanish. This however only happens if $\alpha_X > \alpha_Y$ (negative exponent of $x$). I do not see any plausible reason why such a constraint should be necessary (after all, $Z= XY$; both random variables enter $Z$ "the same way", therefore I expect a "symmetric" result). Anyway, going with this assumption I get:
\begin{align*}
f_Z(z) = \frac{\alpha_X \alpha_Y}{\alpha_X - \alpha_Y} z^{-\alpha_Y-1}
\end{align*}
I suspect this to be wrong due to my aforementioned concerns about the limit of the integral and the assumption on the shape parameters. Moreover, there is this very related question where a commentator points out that $\frac{z}{x}$ must be larger than the distribution's minimum value ($1$ in my case). This is clearly violated for sufficiently large $x$.
On the other hand, I think that the general approach is valid. The formula an the Wikipedia page is apparently correct (I even checked the original source (p. 140)). The formula suggested in this question leads me to the same result (after differentiating the resulting CDF – but subject to the same limitations as my initial approach).
Therefore, my question is:

What is the distribution of $Z$ and what is wrong with my initial approach?

 A: As requested I write out a full derivation here.
Given the pdf of two independent random variables $X, Y$ are
$$ f_X(x) = \alpha_X x^{-\alpha_X-1}, ~~ x > 1$$
$$ f_Y(y) = \alpha_Y y^{-\alpha_Y-1}, ~~ y > 1$$
and equals to $0$ otherwise. The pdf of the product $Z = XY$ is given by
$$ f_Z(z)=\int_{-\infty}^{\infty}f_X(x)f_Y\left(\frac {z} {x}\right)\frac {1} {|x|}dx$$
Since the support of $X,Y$ are $(1, +\infty)$, the integrand $\displaystyle f_X(x)f_Y\left(\frac {z} {x}\right)\frac {1} {|x|}$ is non-zero if and only if
$$ f_X(x) > 0 \text{  and  } f_Y\left(\frac {z} {x}\right) > 0 
\iff x > 1 \text{  and  } \frac {z} {x} > 1 
\iff 1 < x < z$$
From this we also verified the obvious fact that the support of $Z$ is also $(1, +\infty)$ i.e. $f_Z(z) > 0 \iff z > 1$. So for the coming calculation, we are only interested in $z > 1$ only. 
Now we can plug in the pdfs of $X, Y$, and only need to calculate when $x \in (1, z)$ as the integrand will vanish outside. The above integral becomes
$$ \begin{align} f_Z(z) 
&= \int_1^z \alpha_X x^{-\alpha_X-1} 
\alpha_Y \left(\frac {z} {x}\right)^{-\alpha_Y-1} \frac {1} {x} dx \\
&= \alpha_X\alpha_Y z^{-\alpha_Y-1} \int_1^z x^{-\alpha_X+\alpha_Y-1}dx \\
&=\begin{cases}
\displaystyle 
\alpha_X\alpha_Y z^{-\alpha_Y-1} \left. \frac {x^{\alpha_Y-\alpha_X}} {\alpha_Y-\alpha_X}\right|_{x=1}^{x=z} 
& \text{when  } \alpha_X \neq \alpha_Y \\
\displaystyle 
\alpha_X\alpha_Y z^{-\alpha_Y-1} \left. \ln x \right|_{x=1}^{x=z} 
&\text{when  } \alpha_X =\alpha_Y \triangleq \alpha \\
\end{cases} \\
&= \begin{cases}
\displaystyle 
\frac {\alpha_X\alpha_Y} {\alpha_Y - \alpha_X} z^{-\alpha_Y-1}(z^{\alpha_Y-\alpha_X}-1)  
& \text{when  } \alpha_X \neq \alpha_Y \\
\displaystyle 
\alpha^2 z^{-\alpha-1} \ln z
&\text{when  } \alpha_X =\alpha_Y \triangleq \alpha \\
\end{cases} \\
&= \begin{cases}
\displaystyle 
\frac {\alpha_X\alpha_Y} {\alpha_Y - \alpha_X} (z^{-\alpha_X-1}-z^{-\alpha_Y-1})  
& \text{when  } \alpha_X \neq \alpha_Y \\
\displaystyle 
\alpha^2 z^{-\alpha-1} \ln z
&\text{when  } \alpha_X =\alpha_Y \triangleq \alpha \\
\end{cases}
\end{align}$$
and $0$ otherwise. As a careful check, note that when $z > 1$, $\ln z > 0$ and 
$$ z^{-\alpha_X-1}-z^{-\alpha_Y-1} > 0 \iff \alpha_Y - \alpha_X > 0$$
So the calculated pdf is indeed positive for $z > 1$.
