What is the probability density function of $X+Y-Z$, where $X$, $Y$, and $Z$ are three independent and non identical exponential random variables?

Question

Suppose that X, Y and Z are three exponential random variables, what is the probability density function of their combination X+Y-Z?

I've already asked for this here but, as you can see in the comments, when I'm programming the CDF I obtain results that are not acceptable

Answer 1

Assume $P_X(x)=\lambda_xe^{-x\lambda_x},P_Y(y)=\lambda_ye^{-y\lambda_y},P_Z(z)=\lambda_ze^{-z\lambda_z}$ and $x,y,z,\lambda_x,\lambda_y,\lambda_z\ge0$. We need $F_W(w) = \text{Pr}\{W\le w\}$ where $W = X+Y-Z$. Therefore for a given $W=w$, we have to integrate the 3-dimensional triangular polyhedron formed by the plane $X+Y-Z=W$ in the non-negative orthant. We have

$$F_W(w)=\text{Pr}\{X+Y-Z\le w\}=\int_{0}^{\infty}\int_{0}^{w+z}\int_{0}^{w-y+z}P_{X,Y,Z}(x,y,z)dzdydx \\ = \int_{0}^{\infty}P_Z(z)dz\int_{0}^{w+z}P_Y(y)dy\int_{0}^{w-y+z}P_X(x)dx $$

Where $P_{X,Y,Z}(x,y,z) = P_X(x)P_Y(y)P_Z(z)$ comes from the independence.

$$\Rightarrow F_W(w) = \int_{0}^{\infty}\lambda_ze^{-z\lambda_z}dz\int_{0}^{w+z}\lambda_ye^{-y\lambda_y}dy\int_{0}^{w-y+z}\lambda_xe^{-x\lambda_x}dx \\ =\int_{0}^{\infty}\lambda_ze^{-z\lambda_z}dz\int_{0}^{w+z}\lambda_ye^{-y\lambda_y}\left(-e^{-\lambda_x(w+z-y)}+1\right)dy \\ =\int_{0}^{\infty}\lambda_ze^{-z\lambda_z}\left[\frac{\lambda_y}{\lambda_x-\lambda_y}e^{-\lambda_x(w+z)} \left( -e^{(\lambda_x-\lambda_y)(w+z)} +1\right) + \left( -e^{-\lambda_y(w+z)}+1 \right)\right]dz \\ =1-\frac{\lambda_z}{\lambda_y+\lambda_z}e^{-\lambda_yw}+\frac{\lambda_y\lambda_z}{(\lambda_x-\lambda_y)(\lambda_x+\lambda_z)}e^{-\lambda_xw} - \frac{\lambda_z\lambda_y}{(\lambda_x-\lambda_y)(\lambda_z+\lambda_y)}e^{-\lambda_yw} \\ =1+\frac{\lambda_y\lambda_z}{(\lambda_x-\lambda_y)(\lambda_x+\lambda_z)}e^{-\lambda_xw}-\frac{\lambda_z\lambda_x}{(\lambda_y+\lambda_z)(\lambda_x-\lambda_y)}e^{-\lambda_yw} $$

Finally from $f_W(w)=\frac{d}{dw}F_W(w)$ we have

$$f_W(w)=\frac{-\lambda_y\lambda_z\lambda_x}{(\lambda_x-\lambda_y)(\lambda_x+\lambda_z)}e^{-\lambda_xw}+\frac{\lambda_y\lambda_z\lambda_x}{(\lambda_y+\lambda_z)(\lambda_x-\lambda_y)}e^{-\lambda_yw}\\ $$

This holds true when $w \ge 0$. When $w < 0$ we have

$$F_W(w)=\text{Pr}\{X+Y-Z\le w\}= \int_{-w}^{\infty}P_Z(z)dz\int_{0}^{w+z}P_Y(y)dy\int_{0}^{w-y+z}P_X(x)dx $$

Which is similar to previous one and yields

$$F_W(w)=\int_{-w}^{\infty}\lambda_ze^{-z\lambda_z}\left[\frac{\lambda_y}{\lambda_x-\lambda_y}e^{-\lambda_x(w+z)} \left( -e^{(\lambda_x-\lambda_y)(w+z)} +1\right) + \left( -e^{-\lambda_y(w+z)}+1 \right)\right]dz\\ =e^{\lambda_zw}-\frac{\lambda_z}{\lambda_y+\lambda_z}e^{\lambda_zw}+\frac{\lambda_y\lambda_z}{(\lambda_x-\lambda_y)(\lambda_x+\lambda_z)}e^{\lambda_zw} - \frac{\lambda_z\lambda_y}{(\lambda_x-\lambda_y)(\lambda_z+\lambda_y)}e^{\lambda_zw}\\ =\frac{\lambda_x\lambda_y}{(\lambda_y+\lambda_z)(\lambda_x+\lambda_z)}e^{\lambda_zw} $$

Please notice that $x,y,z\ge0$ so when $w<0$ we must have $z \in [-w,\infty)$. Then we get

$$f_W(w) = \frac{\lambda_x\lambda_y\lambda_z}{(\lambda_y+\lambda_z)(\lambda_x+\lambda_z)}e^{\lambda_zw}$$

The MATLAB code to use this distribution as a function is presented below. The parameters $\lambda_x,\lambda_y,\lambda_z$ and the vector of required values for $f_W(w)$ i.e. $w$ are given and the value of distribution as a vector the same length of input $w$ will be returned.

function f_w = PDF_W(lambda_x,lambda_y,lambda_z,w)
    f_w = (w>=0).*(-lambda_x*lambda_y*lambda_z)/(lambda_x-lambda_y)/(lambda_x+lambda_z).*exp(-lambda_x.*w)+...
        (w>=0).*(lambda_x*lambda_y*lambda_z)/(lambda_x-lambda_y)/(lambda_y+lambda_z).*exp(-lambda_y.*w)+...
        (w<0).*(lambda_x*lambda_y*lambda_z)/(lambda_x+lambda_z)/(lambda_y+lambda_z).*exp(+lambda_z.*w);
end

Program for C.D.F. is as follows

function P_w = CDF_W(lambda_x,lambda_y,lambda_z,w)
    f_w = (w>=0)+(w>=0).*(*lambda_y*lambda_z)/(lambda_x-lambda_y)/(lambda_x+lambda_z).*exp(-lambda_x.*w)-...
        (w>=0).*(lambda_x*lambda_z)/(lambda_x-lambda_y)/(lambda_y+lambda_z).*exp(-lambda_y.*w)+...
        (w<0).*(lambda_x*lambda_y)/(lambda_x+lambda_z)/(lambda_y+lambda_z).*exp(+lambda_z.*w);
end