Given independent $Z_i \sim N(\mu_i, \sigma_i ), i=1,2$ derive the density of $(Z_1, Z_1 + Z_2)$. Given independent $Z_i \sim N(\mu_i, \sigma_i ), i=1,2$ I would like to derive the distribution of $(Z_1, Z_1 + Z_2)$ and doing so through deriving a density (but not using Characteristic functions).
My idea is to show explicitly that for some function $f(x, y)$ one has that
$$
P \left [ Z_1 \in (-\infty, a], Z_1 + Z_2 \in (-\infty , b) \right ] = \int_{-\infty}^b \int_{-\infty}^a f(x, y) dx dy.
$$
This would be sufficient to show that $f$ is the density of $(Z_1, Z_1 + Z_2)$.
We may take as known that the sum $Z_1 + Z_2$ of two independent normally distributed random variables follows a $N(\mu_1 + \mu_2, \sigma_1 + \sigma_2)$-distribution, that the joint density two of independent random variables is the product of their individual densities, as well as results from Calculus.

My idea was to do as follows: Considering the following equality of sets
$$
\left \{ Z_1 < a, Z_1 + Z_2 < b \right \} = \left \{ Z_1 < a, Z_2 < b - a \right \}
$$
from which it follows that
\begin{align*}
P \left [ Z_1 < a, Z_1 + Z_2 < b \right ] &= P \left [ Z_1 < a , Z_2 < b - a \right ] \\
&= \int_{-\infty}^{b-a} \int_{-\infty}^a f_{Z_1}(x)f_{Z_2}(y) dx dy.
\end{align*}
Substituting with the transformation $\phi(x, y) = (x, y -a)$ we have that
$$
\int_{-\infty}^{b -a}\int_{-\infty}^a f_{Z_1}(x)f_{Z_2}(y)dx dy = \int_{-\infty}^a\int_{-\infty}^b  f_{Z_1}(x)f_{Z_2}(y-a) dx dy,
$$
where
$$
f_{Z_2}(y-a) = \frac{1}{\sqrt{2 \sigma_2 \pi}}\exp \left ( -\frac{1}{2}\left (\frac{x- (\mu_2 + a)}{\sqrt{\sigma_2}} \right )^2 \right )
$$
But here I get stuck and do not know how to proceed. Is it possible to do a derivation of the density in this - or a similar - manner? [I believe the general approach to prove that $(Z_1, Z_1 + Z_2)$ follows a normal distribution would be to derive the Characteristic function of the Joint normal distribution, show that the Characteristic function of $(Z_1, Z_1 + Z_2)$ is that of a joint normal distribution, and then refer to the uniqueness of Characteristic functions.]
Thanks in advance!
 A: Observe that $\begin{bmatrix} 1&0\\1&1 \end{bmatrix}\cdot \begin{bmatrix}Z_1\\Z_2\end{bmatrix}=\begin{bmatrix}Z_1\\Z_1+Z_2\end{bmatrix}$, by Lemma 1 we get that $$\begin{bmatrix}Z_1\\Z_1+Z_2\end{bmatrix}\sim\mathcal N\left( \begin{bmatrix} 1&0\\1&1 \end{bmatrix} \cdot \begin{bmatrix}\mu_1\\\mu_2\end{bmatrix}, \begin{bmatrix} 1&0\\1&1 \end{bmatrix}\cdot \begin{bmatrix} \sigma_1^2&0\\0&\sigma_2^2 \end{bmatrix}\cdot\begin{bmatrix} 1&0\\1&1 \end{bmatrix}^T \right)$$
which is exactly $\mathcal N\left( \begin{bmatrix} \mu_1\\\mu_1+\mu_2 \end{bmatrix},\begin{bmatrix} \sigma_1^2&\sigma_1^2\\\sigma_1^2&\sigma_1^2+\sigma_2^2 \end{bmatrix} \right)$.


Lemma 1: If $Z\sim \mathcal N(\mu,\Sigma)$ is a jointly gaussian vector of size $n$ and $A$ is a full rank $n\times n$ matrix, then $AZ\sim\mathcal N(A\mu,A\Sigma A^T)$.

We compute directly for a a set $B\subseteq \mathbb R^n$ that you mention ($B=\prod_{i=1}^m ]-\infty,a_i ]$, actually this also holds for any Riemann integrable set)
\begin{align*}
\mathbb P[AZ\in B]&=\mathbb P[Z\in A^{-1} B]\\
&=\int_{A^{-1}B} f_Z(z) dz\\
&=\int_{B} f_Z(A^{-1}y) \det(A^{-1}) dy
\end{align*}
By a change of variable in integration. It remains to compute the density it corresponds to
\begin{align*}
f_Z(A^{-1}y) \det(A^{-1})&=\det(A)^{-1}\det(2\pi\Sigma)^{-\frac{1}{2}}\exp\left(-\frac{1}{2}(A^{-1}y-\mu)^T\Sigma^{-1}(A^{-1}y-\mu)\right)\\
&=\det(A)^{-\frac{1}{2}} \det(2\pi\Sigma)^{-\frac{1}{2}} \det(A^T)^{-\frac{1}{2}}\exp\left(-\frac{1}{2}(y-A\mu)^TA^{-T}\Sigma^{-1}A^{-1}(y-A\mu)\right)\\
&=\det(2\pi A\Sigma A^T)^{-\frac{1}{2}}\exp\left(-\frac{1}{2}(y-A\mu)^T(A\Sigma A^T)^{-1}(y-A\mu)\right)\\
\end{align*}
which is the density of a joint Gaussian vector $\mathcal N(A\mu,A\Sigma A^T)$.
