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!
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)$.
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)$.