Consider the model

\[d X_t = \alpha \left(\theta - X_t\right)dt + \sigma dB_t \]

By Ito Isometry, we have:

\[\langle X \rangle_t = E_t(X) = \theta \left(1 - e^{-\alpha t} \right) + X_0 e^{-\alpha t}\]

and

\[\langle X^2 \rangle_t - \langle X \rangle_t^2 = V_t(X) = \frac{\sigma^2}{2 \alpha }\left(1 - e^{-2 \alpha t} \right)\]

If we assume discrete, uniform sampling of spacing \(\tau\), then we have,

\[P(\theta | X; \alpha, \sigma) = \sqrt{\frac{1}{2\pi V_{\tau} }}^{T-1} \exp\left(-\frac{\sum_t^{T-1} \left(X_{t+1} - E\right)^2 }{2V}\right) \]

or more explicitly (taking \(\tau_i = 1\) without loss of generality),

\[P(\theta | X; \alpha, \sigma) = \sqrt{\frac{\tfrac{\alpha}{\sigma}}{\pi (1- e^{-\alpha t}) }}^{T-1} \exp\left( - \frac{\sum_t^{T-1} \left(X_{t+1} - ( \theta (1 - e^{-\alpha}) + X_0 e^{-\alpha}) \right)^2 }{\tfrac{\sigma^2}{\alpha}(1-e^{-2\alpha})}\right) \]

To integrate out \(\theta\), \(P(X | \alpha, \sigma) = \int P(X | \theta, \alpha, \sigma ) P(\theta) d\theta\), we’ll make this look like a Gaussian in \(\theta\) by completing the square. First, let us introduce a more compact notation to manipulate terms independent of \(\theta\):

\[A_t := X_{t+1} - X_0e^{-\alpha}\] \[ B := 1-e^{-\alpha} \]

Our sum in the exponent can then be written more succinctly:

\[ \sum_t^{T-1} (A_t - \theta B)^2 \]

Squaring out inside the sum and distributing the summation operator (by linearity) and extracting or summing terms constant in \(t\), we have:

\[ \sum A_t^2 - 2 \theta B \sum A_t + \theta^2 B^2 (T-1) \]

For which we complete the square in \(\theta\),

\[ B^2 (T-1) \left( \left(\theta - \frac{\sum_t A_t}{B(T-1)}\right)^2 + \frac{\sum A_t ^2 - (\sum A_t)^2}{B^2 (T-1)} \right) \]

Which lets us write \(P(\theta | X; \alpha, \sigma)\) as a normal distribution in \(\theta\) times some constant terms:

\[ P(\theta | X; \alpha, \sigma) = \exp \left(\frac{-\sum A_t^2 + (\sum A_t)^2 }{V_{\tau} }\right) \sqrt{\frac{1}{2\pi V_{\tau} }}^{T-1} \int d\theta\exp\left( -\frac{ \left(\theta - \frac{\sum_t A_t}{B(T-1)} \right)^2 }{\tfrac{2V_{\tau}}{B^2 (T-1)}}\right) \]

In which we recognize the integral as Gaussian with mean $= $ and variance \(\nu = V_{\tau}/(B^2(T-1))\), and thus can replace the integral with \(\sqrt{2 \pi \nu} = \sqrt{\frac{2 \pi V_{\tau} }{B^2(T-1)}}\),

\[ P(\theta | X; \alpha, \sigma) = \exp \left(\frac{-\sum A_t^2 + (\sum A_t)^2 }{V_{\tau} }\right) \sqrt{\frac{1}{2\pi V_{\tau} }}^{T-1} \sqrt{\frac{2 \pi V_{\tau} }{B^2(T-1)}}\]

Collecting common terms in \(V_t\),

\[ P(\theta | X; \alpha, \sigma) = \exp \left(\frac{f(A_t)}{V_{\tau}}\right) \sqrt{\frac{1}{2\pi}}^{T-1} \sqrt{\frac{2 \pi }{B^2(T-1)}} V_{\tau}^{\tfrac{1}{2}} V_{\tau}^{\tfrac{1-T}{2}} \]

\[ = \exp \left(\frac{f(A_t)}{V_{\tau}}\right) \sqrt{\frac{(2 \pi)^{2-T}}{B ^ 2(T-1)}} V_{\tau}^{1 - \tfrac{T}{2}} \]

where

\[ f(A_t) = - \sum A_t^2 + \left( \sum A_t \right)^2 \]

Then substituting in \(V_t = \tfrac{\sigma^2}{2\alpha} \left(1 - e^{-2\alpha}\right)\)

\[ = \exp \left(\frac{\tfrac{2\alpha f(A)}{1 - e^{-2\alpha}} }{\sigma^2}\right) \sqrt{\frac{(2 \pi)^{2-T}}{B ^ 2(T-1)}} \left(\tfrac{1}{2\alpha} \left(1 - e^{-2\alpha}\right)\right)^{1 - \tfrac{T}{2}}\sigma^{2(1 - T/2)} \]

Which can be written to look like an Inverse Gamma integral in \(\sigma^2\)

\[ = \mathcal{N} e^{\tfrac{-\beta}{\sigma^2}} (\sigma^2)^{-\gamma - 1} \]

where

\[ \gamma = T/2 -2,\]

\[ \beta = \frac{2 \alpha \left( \left( \sum A_t \right)^2- \sum A_t^2 \right) }{1 - e^{-2\alpha}} \] and

\[ \mathcal{N} = \sqrt{\frac{(2 \pi)^{2-T}}{B ^ 2(T-1)}} \left(\tfrac{1}{2\alpha} \left(1 - e^{-2\alpha}\right)\right)^{1 - \tfrac{T}{2}} \]

Which we can integrate as Inverse Gamma (again assuming uniform prior, now in \(\sigma^2\)),

\[ = \mathcal{N} \int e^{\tfrac{-\beta}{\sigma^2}} (\sigma^2)^{-\gamma - 1} d\sigma^2 = \mathcal{N} \frac{\beta^{\gamma}}{\Gamma(\gamma)} \]

with \(\beta\), \(\gamma\), and \(\mathcal{N}\) as defined above.