The problem

Consider a stochastic process with the following three properties:

  1. 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}).$
  2. 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] \, .$$
  3. 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$.

  • $\begingroup$ This is the multipoint probability density for the ornstein-uhlenbeck process. So one could take shortcuts, calculate average and correlation function explucitly and plug them into the expression. But your approach is correct as well... in fact, I don't see how manipulating a quadratic form can pose a challenge, apart from just having to write a lot (no offense, just trying to encourage). $\endgroup$
    – Roger V.
    Oct 23, 2021 at 4:43
  • $\begingroup$ @RogerVadim Thanks for the encouragement. The problem I ran into with manipulating the quadratic form is that the denominators of each factor are different, and I didn't find a way to go to a common denominator without creating chaos. $\endgroup$
    – DanielSank
    Oct 23, 2021 at 21:43
  • $\begingroup$ I don't think that it is worthwhile bringing them to common denominator, especially if the piints are not equidistant: matrix $\Sigma^{-1}$ is bi- or tridiagonal (symmetric form is needed for using gaussian integral identity); you need to find its elements $\endgroup$
    – Roger V.
    Oct 24, 2021 at 6:22

1 Answer 1


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.

  • $\begingroup$ Is it possible that in the sum over both $i$ and $j$, the lower bound on $j$ is 1 instead of 2? $\endgroup$
    – DanielSank
    Nov 4, 2021 at 6:28
  • $\begingroup$ We can clean this up somewhat by defining $\tau_i = t_i - t_{i-1}$ for $i \geq 2$, $v_0=0$ and $\tau_1 = \infty$. Would you like me to edit for you? $\endgroup$
    – DanielSank
    Nov 4, 2021 at 22:29
  • $\begingroup$ I was wondering whether there would be an expression for the probability density of branching times under a Ornstein–Uhlenbeck process, without knowing the spatial position (or state value $x_i$) of the branching point. I just posted a specific question here: physics.stackexchange.com/questions/739347/… $\endgroup$ Dec 2, 2022 at 15:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.