A Multivariate Distribution for Linear Combinations of Independent Exponential Random Variables Consider a random vector $\mathbf{X} \in \mathbb{R}^r$ whose components $X_j$ are independent exponential variables with different scale parameters $\beta_j$, $j=1,\dots,r$. Suppose I have a general $p \times r$ ($p < r$) matrix $\mathsf{M}$ with non-negative entries $M_{ij} \geqslant 0$.
What would be the multivariate probability density function for this random vector?
\begin{equation}
\mathbf{Y} = \mathsf{M} \mathbf{X} \in \mathbb{R}^p
\end{equation}
[It is essentially $p$ positive linear combinations of the $r$ exponential variables $X_j$.]
I know that each component $Y_i$ ($i=1,\dots,p$) would follow a hypo-exponential distribution, but I have no idea what a multivariate hypo-exponential distribution is like. If we are willing to approximate each hypo-exponential component $Y_i$ as a gamma random variable, then I am looking for a version of the multivariate gamma distribution given the covariance matrix (constructed from $\mathsf{M}$).
I have searched for many multivariate gamma distributions in the literature, but none of them seems to model this case. Any help is much appreciated!
Edits. I have difficulty directly deriving the density for $\mathbf{Y}$ from $\mathbf{X}$ because writing down the density in terms of $\mathbf{Y}$ requires the inversion of a non-square matrix $\mathsf{M}$. Padding $\mathsf{M}$ into a square matrix by introducing $(r-p)$ helper coordinates in $\mathbf{Y}$ requires marginalising these components out which I want to avoid doing.
 A: I will do it first for the case $r=p$ and then we can use that case to solve for the case of main interest, $p < r$. To simplify notation I will use the rate parameters $\lambda_j = 1/\beta_j$. Also, I will write $M^T$ (transpose) for your $M$. So we have $Y=M^T X$ where $X$ is a $r\times 1$ random vector and $M$ is a $r\times p$ matrix (with $m_{ij}\ge 0$, and of full rank, and for the moment, $p=r$). I will write the density as a differential form, simplifying the change of variables. This is the notation used in Muirhead: "Aspects of Multivariate Statistical Theory".  So we write the exterior product of the $r$ components of $d x$ as $\left( d x\right) = \wedge_{j=1}^r d x_j$.
By independence we can write the density of $X$ as 
$$
   f(x) =\prod^r \lambda_j \prod^r e^{-\lambda_j x_j} \, dx_j \\
= \prod^r \lambda_j e^{-\lambda^T x} \left( dx \right), \quad x_j>0
$$
where $\lambda=(\lambda_1, \dotsc,\lambda_r)$ and $x=(x_1, \dotsc, x_r)$. Then we substitute $x=M^{-T}y$ to get the density of $y$:
$$
   g(y)=\prod^r \lambda_j e^{-y^T M^{-1}\lambda} \left( d(M^{-T}y\right) \\
= \prod^r \lambda_j e^{-y^T M^{-1}\lambda} \det(M^{-1})\left( dy \right) \\
= \frac{\prod^r \lambda_j}{\det(M)} e^{-y^T M^{-1}\lambda} \left( dy \right),
\quad y_j>0
$$
Then the case $p<r$. We still assume that $M$ is of full rank. Let $m=r-p$. We will extend the $p$-vector $y$ with $m$ new coordinates, so we can use the square case. It will be useful to use the QR-decomposition of $M$, $M=QR=[Q_1 \colon Q_2] \begin{bmatrix} R_1 \\ 0\end{bmatrix}$ and then define 
$M^* =[Q_1 R_1\colon Q_2]= [Q_1 \colon Q_2]\begin{bmatrix} R_1 & 0 \\ 0& I_m  \end{bmatrix}$. Since the first factor is an orthogonal matrix with determinant $\pm 1$, it follows that the absolute value of the determinant of $M^*$ is the determinant of $R_1$, which is simply the product of its diagonal values. And since the inverse of an orthogonal matrix is simply it transpose, we get that
$$
   M^{*-1}= \begin{bmatrix} R_1 & 0 \\ 0& I_m  \end{bmatrix}^{-1}[Q_1 \colon Q_2]^{-1} = \begin{bmatrix} R_1^{-1} & 0 \\ 0& I_m  \end{bmatrix}\begin{bmatrix}   Q_1^T\\ Q_2^T \end{bmatrix}= \begin{bmatrix} R_1^{-1}Q_1^T\\Q_2^T  \end{bmatrix}
$$
Write the new $y$-vector with the $m$ new components as $y^*=[y\colon \bar{y}]$ with $\bar{y}=Q_2^T x$. Then we must write the density of $Y^*$ as a function of $y$ and $\bar{y}$, and integrate out $\bar{y}$:
$$
 g(y^*) = \frac{\prod^r \lambda_j}{\det(M^*)} e^{-(y,\bar{y})^T M^{*-1}\lambda} \left( dy^* \right)  \\
= \frac{\prod^r \lambda_j}{\det(R_1)} e^{-y^T R_1^{-1} Q_1^T \lambda} \left( dy \right) \int_0^\infty \dotsm \int_0^\infty e^{-\bar{y}^T Q_2^T \lambda }\left( d \bar{y}\right)
$$
The last $m$-fold integral is simply a product of $m$ exponential integrals. Write $\theta_j$ for the $j$th component of $Q_2^T \lambda$, and the last integral becomes $1/\prod^m \theta_j$. Then we have finally the density of $Y$ as 
$$
  g(y) = \frac{\prod^r \lambda_j}{\det(R_1) \prod^m \theta_j } e^{-y^T R_1^{-1} Q_1^T \lambda} \left( dy \right) , \quad y_j>0 
$$
(really hope I got everything right!)
