Ornstein–Uhlenbeck process: joint probability as a Gaussian The problem
Consider a stochastic process with the following three properties:

*

*The process is Markov, meaning that $p(x_n,t_n|x_{n-1},t_{n-1},\ldots x_1, t_1) = p(x_n,t_n|x_{n-1},t_{n-1}).$

*The conditional probability is
$$p(x_n, t_n|x_{n-1},t_{n-1}) = \left[ 2 \pi \sigma^2 (1 - e^{-2 \gamma (t_n - t_{n-1})})\right]^{-1/2} \exp \left[- \frac{\left(x_n - x_{n-1}e^{-\gamma(t_n - t_{n-1})} \right)^2}{2 \sigma^2 \left( 1 - e^{-2 \gamma (t_n - t_{n-1})}\right)} \right] \, .$$

*The unconditioned probability is
$$\bar{p}(x) = \left[ 2 \pi \sigma^2\right]^{-1/2} \exp \left[- \frac{x^2}{2 \sigma^2}\right] \, .$$
I'm reading Balakrishnan's Elements of Nonequilibrium Statistical Mechanics in which one of the exercises is to show that for this process, the n-point probability density is
$$
p(x_n,t_n,x_{n-1},t_{n-1},\ldots,x_1,t_1) =
(2 \pi)^{-n/2} \text{det}(\mathbf{\Sigma})^{-1/2}
\exp\left[ - \frac{1}{2} \mathbf{x}^\intercal (\mathbf{\Sigma}^{-1}) \mathbf{x}\right] \tag{$\star$}
$$
where $\mathbf{\Sigma}$ is a positive definite matrix and $\mathbf{x}$ represents a vector of the random variables $\{x_1, x_2,\ldots,x_n\}$.
I would like to know how to show this.
What I've tried
One approach could be to make the demonstration explicitly.
We can use some basic probability and the Markov property to write
\begin{align}
p(x_n, t_n, x_{n-1}, t_{n-1},\ldots x_1, t_1) &=
p(x_n,t_n|x_{n-1},t_{n-1},\ldots x_1, t_1) \times p(x_{n-1},t_{n-1},\ldots x_1, t_1)\\
\text{(use Markov)} \qquad &= p(x_n,t_n|x_{n-1},t_{n-1}) \times p(x_{n-1},t_{n-1},\ldots x_1, t_1) \\
&= p(x_n,t_n|x_{n-1},t_{n-1}) \times p(x_{n-1},t_{n-1}| x_{n-2}, t_{n-2}) \times p(x_{n-2},t_{n-2},\ldots x_1, t_1) \\
\text{(repeat)} \qquad &= \bar{p}(x_1) \prod_{i=2}^n p(x_i, t_i|x_{i-1} t_{i-1}) \, .
\end{align}
From here I tried plugging in the explicit forms of the two-point conditional probability.
It looks like there should be enough structure to make progress.
For example, the product of exponentials can be written as a single exponential of a sum.
So if we can make that sum look like the matrix contraction we're going for, then we'd be done (assuming that the product of the prefactors works out to the intended determinant).
However, I have not yet succeeded to make progress in this direction. The main problem is that the thing in the exponentials is a fraction and they have all different denominators.
Another thing I considered was to somehow just prove that the stipulated form $(\star)$ has the right properties, i.e. yields the right moments, to match the three properties we're given.
I think this might be a good approach but I may need some hints about how to do the necessary Gaussian integrals, e.g. to compute things like $\langle x_i x_j \rangle$.
 A: The (conditional) distribution is always Gaussian, so one could just state the result as obvious, given we're not, in the question as worded, interested in writing $\mathbf{\Sigma}$ in terms of the O-U parameters.
Now that's not maybe very satisfying, so instead let's get the covariance matrix as a function of the observation times, O-U parameters.
\begin{align}
p(x_n, t_n,\ldots x_1, t_1)
&= \bar{p}(x_1) \prod_{i=2}^n p(x_i, t_i|x_{i-1} t_{i-1}) \\
&= \underbrace{\left( 2 \pi \sigma^2\right)^{-1/2} \exp\left(-\frac{x_1^2}{2\sigma^2} \right)}_{\bar{p}(x_1)}
\prod_{i=2}^n\left( 2 \pi \sigma^2 (1 - e^{-2 \gamma \tau_i})\right)^{-1/2} \\
&\times \exp \left(- \frac{1}{2}\sum_{i=2}^n \frac{\left(x_i - x_{i-1}e^{-\gamma \tau_i} \right)^2}{\sigma^2 \left( 1 - e^{-2 \gamma \tau_i}\right)}\right)
\end{align}
where $\tau_i \equiv t_i - t_{i-1}$.
If we define $\tau_1 = \infty$ and $x_0 = 0$, which corresponds to the first sample point having no memory of the past (i.e. the first sample point is not conditional on any previous data), then we can simplify the expression to
\begin{align}
p(x_n, t_n,\ldots x_1, t_1)
&= \prod_{i=1}^n\left( 2 \pi \sigma^2 (1 - e^{-2 \gamma \tau_i})\right)^{-1/2} \\
&\times \exp \left(- \frac{1}{2}\sum_{i=1}^n \frac{\left(x_i - x_{i-1}e^{-\gamma \tau_i} \right)^2}{\sigma^2 \left( 1 - e^{-2 \gamma \tau_i} \right)} \right)
\end{align}
Simplifying the notation with
\begin{align}
  a_i &\equiv \frac{1}{\sigma^2 \left( 1 - e^{2 \gamma \tau_i} \right)} \\
  b_i &\equiv e^{-\gamma \tau_i} 
\end{align}
so that the sum inside the exponential is
\begin{align}
\sum_{i=1}^n a_i(x_i - b_i x_{i-1})^2
&= \sum_{i=1}^n a_i(x_i^2 - 2b_i x_ix_{i-1} + b_i^2x_{i-1}^2) \\
\end{align}
which can be expressed as a contraction on a symmetric matrix $A$, i.e. $\mathbf{x}^\intercal A \mathbf{x}$, where
$$
A = \begin{pmatrix}
a_1 + a_2 b_2^2 & -a_2 b_2 & & & & \\
-a_2 b_2 & a_2 + a_3 b_3^2 & & & & \\
& & \ddots \\
& & & a_{n-1} + a_n b_n^2 & -a_n b_n \\
& & & -a_n b_n & a_n
\end{pmatrix}
\, .
$$
Therefore we've shown that
\begin{align}
p(x_n, t_n,\ldots x_1, t_1)
&= \prod_{i=1}^n \left(2\pi \sigma^2 \left( 1 - e^{-2 \gamma \tau_i} \right) \right)^{-1/2}
\exp \left[ -\frac{1}{2} \mathbf{x}^\intercal A \mathbf{x} \right] \\
&= (2 \pi)^{-n/2} \prod_{i=1}^n \left(\sigma^2 \left( 1 - e^{-2 \gamma \tau_i} \right) \right)^{-1/2}
\exp \left[ -\frac{1}{2} \mathbf{x}^\intercal A \mathbf{x} \right] \\
&= (2 \pi)^{-n/2} \left( \prod_{i=1}^n a_i \right)^{-1/2}
\exp \left[ -\frac{1}{2} \mathbf{x}^\intercal A \mathbf{x} \right]
\, .
\end{align}
which matches the desired form with $A = \Sigma^{-1}$, as long as
$$
\text{det}A = \prod_{i=1}^n a_i
$$
which remains to be shown.
