I used the inversion method to find an algorithm that can generate pareto-distributed pseudo-random numbers, by first acknowledging that the pareto distribution is $$ f(x)=\frac{\alpha x_0^\alpha}{x^{\alpha+1}} $$ where $x_0>0, \alpha >0$ are the scale and shape parameters respectively.
The algorithm is basically to 1) Draw $u \sim U[0,1]$, and 2) $x=F^{-1}(u)$, where we get $F^{-1}(u)$ from
$$ F(x) = \int_{x_0}^x f(y) dy=1-\Big( \frac{x_0}{x} \Big)^\alpha \implies x=F^{-1}(u)=\frac{x_0}{(1-u)^{1/\alpha}} $$
Step 1 of the algorithm is very basic, and step 2 is now finished since we have derived a formula for $x$, depending on $u$. But how do I prove that the output of the algorithm is pareto distributed, as desired? I thought of inserting the $x$ into the expression for $f(x)=f(\frac{x_0}{(1-u)^{1/\alpha}})$, but that didn't really get me anywhere.
A simple proof (for any cumulative distribution function $F$) can be found here. In your case you don't need to take the infinum in the definition of $F^{-1}$, since your density is strictly positive and thereby $F$ is bijective (if defined on $[x_0,\infty)$), so things get even simpler.
The proof shows that your constructed random variable $F^{-1}(u)$ has the correct cumulative distribution function and thereby the correct density.