Difference between two Gaussian stochastic variables Suppose I have a sequence of time steps $t_1 < t_2 < t_3< \dots<t_n$. To each time step corresponds a random Gaussian variable $X_1,X_2,X_3,\dots,X_n$. Since $X = \{X_1,X_2,X_3,\dots,X_n\}$ is normally distributed, I can say that
$$ X \sim N \left( 0, \sigma^{2} \right) $$
where the distribution of $X$ has zero mean (assumption) and standard deviation $\sigma$. If I define some stochastic variable $Y_j$ such that
$$ Y_j = X_{j} - X_{j-1} $$
with $Y=\{Y_1,Y_2,\dots, Y_n\}$ for $2 \leq j \leq n$, what can I say about the distribution of $Y$? Naively, I'd think that
$$ Y \sim N \left( 0, \sigma^{2} \right) $$
which means that $Y$ is distributed the same way as $X$. My intuitive understanding is that since any $X$ is normally distributed, then any linear combination of $X$ is also normally distributed with the same mean and variance. Is this correct? If so, how can I show this rigorously?
 A: Let $X\sim\mathcal N(0,\sigma^2I_n)$ be a vector of $n$ IID centered normal variables with variance $\sigma^2$ and let $Y=AX$, where $I_n$ is the identity matrix and
$$
A=\begin{pmatrix}
 1 &  0 &  0 & \dots & \dots & 0 & 0 & 0\\
-1 &  1 &  0 & \dots & \dots & 0 & 0 & 0\\
 0 & -1 &  1 & \dots & \dots & 0 & 0 & 0\\
 0 &  0 & -1 & \ddots & \dots & 0 & 0 & 0\\
 \vdots & \vdots & \vdots & \ddots & \ddots & \vdots & \vdots & \vdots\\
 0 &  0 &  0 & \dots & \ddots & 1 & 0 & 0\\
 0 &  0 &  0 & \dots & \dots & -1 & 1 & 0\\  
 0 &  0 &  0 & \dots & \dots & 0 & -1 & 1
\end{pmatrix}.
$$
This gives $Y_1=X_1$ and $Y_j=X_j-X_{j-1}$ otherwise.
Using integration by substitution for $X=A^{-1}Y$ (notice that $|A|=1$) and $x^{\mathrm t}$ to denote the transpose, we have
\begin{aligned}
\mathbb P(Y\in\mathcal E)&=
\int\unicode{120793}\{Ax\in\mathcal E\}\frac{1}{\sqrt{(2\pi)^n\sigma^2}}e^{-\frac{1}{2}x^{\mathrm{t}}(\sigma^2I_n)^{-1}x}\mathrm dx\\
&=\int\unicode{120793}\{x\in\mathcal E\}\frac{|A^{-1}|}{\sqrt{(2\pi)^n\sigma^2}}e^{-\frac{1}{2}x^{\mathrm{t}}\Sigma^{-1}x}\mathrm dx,
\end{aligned}
where $\Sigma^{-1}=(A^{-1})^{\mathrm{t}}(\sigma^2I_n)^{-1}A^{-1}=(A(\sigma^2I_n)A^{\mathrm t})^{-1}=(\sigma^2AA^{\mathrm{t}})^{-1}$.
Recall the multiplicativity of the determinant, which gives $|A^{-1}|=1$, and thereby $Y\sim\mathcal N(0,\sigma^2AA^{\mathrm{t}})$. This shows that $Y$ is a multivariate normal vector. Clearly, we want to know the covariance exactly, so
$$
\Sigma=\sigma^2AA^{\mathrm{t}}=\sigma^2
\begin{pmatrix}
 1 &  -1 &  0 & \dots & \dots & 0 & 0 & 0\\
 -1 &  2 &  -1 & \dots & \dots & 0 & 0 & 0\\
 0 &  -1 &  2 & \dots & \dots & 0 & 0 & 0\\
 0 &  0 &  -1 & \ddots & \dots & 0 & 0 & 0\\
 \vdots & \vdots & \vdots & \ddots & \ddots & \vdots & \vdots & \vdots\\
 0 &  0 &  0 & \dots & \ddots & 2 & -1 & 0\\
 0 &  0 &  0 & \dots & \dots & -1 & 2 & -1\\  
 0 &  0 &  0 & \dots & \dots & 0 & -1 & 2
\end{pmatrix}.
$$
Normal variables are independent iff they are uncorrelated. Thus, for example, $Y_n$ is independent of $(Y_m)_{m\le n-2}$, but clearly, the components of $Y$ are not mutually independent (since not all off-diagonal entries are $0$).
A: Credits to @Matija for the answer, but here's my shorter take on it.
Suppose $X\sim N(0,1)$. If $Y=AX+b$, then $Y\sim N(b,A^{T}A)$
Proof:
First, denote the PDF of $X$ as $f_{X}(x)$ such that
$$
f_{X}(x) = \frac{1}{\sqrt{(2\pi)^{n}}}\text{Exp}\left(-\frac{1}{2\sigma^{2}}X^{T}X\right)
$$
One can ask, given $f_{X}(x)$, what is the PDF of $Y$? Namely $f_{Y}(y)$? I claim (without proof) that if $Y=f(X)$, then injectively, $X=f^{-1}(Y)=G(Y)$ where $f^{-1}\equiv G$. Then
$$
f_{Y}(y) = f_{X}(G(y))|J|
$$
where $J$ is the Jacobian, with matrix elements defined as
$$
J_{ij}=\frac{\partial{G_{i}}}{\partial Y_{j}}
$$
Since $Y$ is a simple linear combination of $X$, $J$ can be easily computed:
$$
Y = Ax+b \\
X = A^{-1}(Y-b) \equiv G(Y)
$$
Hence $\partial G/\partial Y = A^{-1}$ and $|J| = |A^{-1}| = |A|^{-1}$. Since we know the PDF for $X$, we can calculate $f_{X}(G(y))$:
$$
f_{X}(G(y)) = \frac{1}{\sqrt{(2\pi)^{n}}}\text{Exp}\left(-\frac{1}{2\sigma^{2}}\left(A^{-1}(Y-b)\right)^{T}\left(A^{-1}(Y-b)\right)\right) \\
=\frac{1}{\sqrt{(2\pi)^{n}}}\text{Exp}\left(-\frac{1}{2\sigma^{2}}\left((Y-b)^{T}A^{-T}\right)\left(A^{-1}(Y-b)\right)\right)\\
=\frac{1}{\sqrt{(2\pi)^{n}}}\text{Exp}\left(-\frac{1}{2\sigma^{2}}(Y-b)^{T}(A^{T}A)^{-1}(Y-b)\right)\\
=\frac{1}{\sqrt{(2\pi)^{n}}}\text{Exp}\left(-\frac{1}{2\sigma^{2}}(Y-b)^{T}(C)^{-1}(Y-b)\right)
$$
where $C=A^{T}A$. Note that
$$
|C| = |A^{T}A| = |A^{T}||A|=|A|^{2}
$$
hence
$$
|A| = \sqrt{|C|} \longrightarrow |A|^{-1} = \frac{1}{\sqrt{|C|}}
$$
Therefore, the final PDF of $Y$ gives
$$
f_{Y}(y) = f_{X}(G(y))|J|\\
=\frac{1}{\sqrt{|C|}}f_{X}(G(y))\\
=\frac{1}{\sqrt{(2\pi)^{n}|C|}}\text{Exp}\left(-\frac{1}{2\sigma^{2}}(Y-b)^{T}(C)^{-1}(Y-b)\right)
$$
But this is simply the PDF for a multivariate distribution with mean $b$ and covariance matrix $C$. Hence I can say that
$$
Y\sim N(b,C)
$$
