Let $X_1$,$X_2$, $X_3$ denote a random sample from the distribution having p.d.f $$f(x) = e^{-x},\,\, 0<x <\infty,$$ zero elsewhere. Show that $$y_1 = \frac{x_1}{x_1 + x_2} $$ $$y_2 = \frac{x_1 +x_2}{x_1 + x_2 +x_3} $$
$$y_3 = {x_1 + x_2 +x_3} $$
are mutually stochastically indepedent.
According to the definition of continuous independent random variables, the joint pdf of these random variables is equal to the product marginal pdf of these variables
The joint can be can be calculated using the Jacobian, but first I have to solve the linear system of equations The answer is $$x_1 = {y_1y_2y_3} $$ $$x_2 = {y_2y_3}-{y_1y_2y_3} $$ $$x_3 = {y_3}-{y_2y_3} $$
and the determinant of the Jacobian is
$${y_2}{y_3^2} $$
How do I obtain the joint pdf of these functions? How do I obtain the marginal of each variable; which limits should I take to integrate?
Using theorem 4 in the source you posted,
$$f_{Y_1Y_2Y_3}(y_1,y_2,y_3) = f_{X_1X_2X_3}(x_1,x_2,x_3)\lvert J \rvert$$
where
$$ J = \begin{vmatrix} \frac{\partial x_1}{\partial y_1} & \frac{\partial x_1}{\partial y_2} & \frac{\partial x_1}{\partial y_3}\\ \frac{\partial x_2}{\partial y_1} & \frac{\partial x_2}{\partial y_2} & \frac{\partial x_2}{\partial y_3}\\ \frac{\partial x_3}{\partial y_1} & \frac{\partial x_3}{\partial y_2} & \frac{\partial x_3}{\partial y_3} \end{vmatrix} $$
Asumming that $X_1$, $X_2$ and $X_3$ are independent,
$$f_{X_1X_2X_3}(x_1,x_2,x_3) = e^{-(x_1+x_2+x_3)} = e^{-y_3} \qquad \text{for }y_3>0$$
and as you said, $\lvert J \rvert = y_2y_3^2$. Therefore,
$$f_{Y_1Y_2Y_3}(y_1,y_2,y_3) = y_2y_3^2e^{-y_3} \qquad \text{for }y_3>0, 0\lt y_1,y_2 \lt 1$$
To better interpret this result, we can write
$$f_{Y_1Y_2Y_3}(y_1,y_2,y_3) = f_{Y_2Y_3\mid Y_1}(y_2,y_3\mid y_1)f_{Y_1}(y_1) $$
And from this we can see that
$$f_{Y_2Y_3\mid Y_1}(y_2,y_3\mid y_1) = f_{Y_2Y_3}(y_2,y_3) = y_2y_3^2e^{-y_3} \qquad \text{for }y_3>0, 0\lt y_2 \lt 1$$
and
$$f_{Y_1}(y_1) = 1\qquad \text{for }0\lt y_1 \lt 1$$
The marginal of $Y_3$ can be calculated by integrating $f_{Y_2Y_3}(y_2,y_3)$ over $y_2$,
$$f_{Y_3}(y_3) = \int_0^1y_2y_3^2e^{-y_3}dy_2 = \frac{1}{2}y_3^2e^{-y_3}\qquad \text{for }y_3>0$$
In the same way,
$$f_{Y_3}(y_3) = \int_0^\infty y_2y_3^2e^{-y_3}dy_3 = 2y_2\qquad \text{for }0\lt y_2 \lt 1$$
Informal note about the limits of integration:
For $Y_1$: By definition,
$$y_1 = \frac{x_1}{x_1 + x_2}.$$
Since $x_1,x_2 \gt 0$, $y_1 \gt 0$. In addition, the denominator is always greater than the numerator, then $y_1 \lt 1$, which give us $0 \lt y_1 \lt 1$.
For $Y_2$: Same argument than for $Y_1$
For $Y_3$: $y_3 \gt 0$ since it is the sum of positive numbers.