\documentstyle[11pt]{article}

\begin{document}

\def\bff{{\bf f}}
\def\bfk{{\bf k}}
\def\bfx{{\bf x}}
\def\bfy{{\bf y}}
\def\bfA{{\bf A}}
\def\bfB{{\bf B}}
\def\bfC{{\bf C}}
\def\bfE{{\bf E}}
\def\bfI{{\bf I}}
\def\bfGamma{{\bf \Gamma}}
\def\bfU{{\bf U}}
\def\bfL{{\bf L}}
\def\bfT{{\bf T}}

\def\tb#1{{\tilde {\bf #1}}}
\def\tbf{\tb{f}}
\def\tbk{\tb{k}}
\def\tbA{\tb{A}}
\def\tbB{\tb{B}}
\def\tbC{\tb{C}}
\def\tbE{\tb{E}}
\def\tbI{\tb{I}}

\def\tr{^T}
\def\inv{^{-1}}
\def\tchi{{\tilde \chi}}
\def\chiprh{\chi_{\rm PRH}}
\def\tQ{{\tilde Q}}
\def\Ahat{{\hat {\bf A}}}
\def\be{\begin{equation}}
\def\bel#1{\begin{equation}\label{eq:#1}}
\def\ee{\end{equation}}
\title{Notes on Gaussian Random Functions \\
with Exponential Correlation Functions\\
(Ornstein-Uhlenbeck Process)}
\author{George B. Rybicki}
\date{2 Dec 1994}

\maketitle

We discuss here the properties of a Gaussian random process $x(t)$ of
a very special type, namely, one that has
zero mean and the exponential correlation function
\bel{1}
    \Phi(\tau) = \langle x(t)x(t+\tau) \rangle = \sigma^2 \exp(-\alpha|\tau|)
\ee
for time lag $\tau$.  The constants $\sigma^2$  and $\alpha$ are, respectively,
the variance and the inverse correlation time of the
process.  (The quantity $\sigma$ itself is the standard deviation.)
This process is known as the {\em Ornstein-Uhlenbeck process}.  (ref.\ ?)

An interesting special case of this process is where
$\alpha$ approaches zero, while $\sigma$ becomes infinite in such
a way that $\alpha\sigma^2$ approaches a fixed constant.  We can
accomplish this by letting
\bel{18}
   \sigma^2 = {D \over 2\alpha}
\ee
where $D$ is a constant, called the {\em diffusivity}.  This limiting case
is often called the Gaussian {\em random walk} process. (ref.\ ?)

We note that the correlation function
(\ref{eq:1}) of the random walk is not defined,
since $\sigma \rightarrow \infty$ as $\alpha \rightarrow 0$.  However, the
{\em structure function}, defined as
\bel{21}
   \psi(\tau)= \langle [x(t)-x(t+\tau)]^2 \rangle 
           = 2[\Phi(0)-\Phi(\tau)]
          = 2\sigma^2[1-\exp(-\alpha|\tau|)]
\ee
does have meaning in the limit, namely,
\bel{22}
   \psi(\tau) = D|\tau|
\ee
This shows directly the meaning of $D$ as a diffusivity.
[Note that the structure function is often defined as 
$V(\tau)=(1/2)\psi(\tau)$.]

\section{Joint probability functions}

Let $t_i$, $1 \le i \le n$ be an ordered set of times,
$t_1 < t_2 < \ldots < t_n$, and let $x_i=x(t_i)$ be the values of the
process at those times.  Then from the definition of a Gaussian
process, we know that the joint probability distribution function of
the set of values $x_i$ is
\bel{2}
  P_n(x_n,\ldots,x_1)=\left[ \det(2\pi \bfC_n) \right]^{-1/2}
      \exp\left(-{1\over 2} \bfx\tr \bfC_n\inv \bfx \right)
\ee
where the $n$-dimensional vector $\bfx$ has components $x_i$ and the
$n\times n$ matrix $\bfC_n$ has components $C_{ij}=\Phi(t_i-t_j)$.
The superscript ``$^T$'' indicates matrix transposition.

It turns out (Rybicki \& Press 1995) that the matrix $\bfC_n$ has a simple
inverse, namely, the tridiagonal matrix
\bel{3}
     \bfT_n=\bfC_n^{-1} = \sigma^{-2}
   \left( \matrix{d_1  & -e_1 &     &&&\cr
                  -e_1  & d_2 & -e_2 &&&\cr
                       & -e_2 & d_3 & -e_3 &&\cr
                       &     &  \ddots&\ddots&\ddots &\cr
                       &     &        &-e_{n-2}& d_{n-1} &-e_{n-1}\cr
                       &     &        &      & -e_{n-1} &d_n  \cr}\right),
\ee
where the elements not indicated are zero.
Defining
\bel{4}
    r_i = \cases{0, &$i=0$, \cr
                \exp[-\alpha(t_{i+1}-t_i)],&$1\le i \le n-1$, \cr}
\ee
the quantities $e_i$ and $d_i$ are given by
\bel{5}
       e_i = \cases{ 0, &$i=0$, \cr
                     r_i/(1- r_i^2), &$ 1 \le i \le n-1$ \cr
                     0, &$i=n$, \cr}
\ee
and 
\bel{6}
       d_i = 1+r_ie_i+r_{i-1}e_{i-1}, \qquad 1 \le i \le n.
\ee

With the representation (\ref{eq:3}) for the inverse, the quadratic
form appearing in the exponent of Eq.\ (\ref{eq:2}) can be written,
\bel{7}
      Q_n \equiv \bfx\tr \bfC_n\inv \bfx 
=\sigma^{-2}\sum_{i=1}^{n} (d_ix_i^2 +2e_{i-1} x_i x_{i-1} )
\ee  
Substituting for $d_i$ gives
\begin{eqnarray}\label{eq:8}
Q_n &=&
\sigma^{-2}\sum_{i=1}^n \left[(1-r_{i-1}e_{i-1})x_i^2+2e_{i-1}x_ix_{i-1} \right]
    -\sigma^{-2}\sum_{i=1}^n r_ie_ix_{i}^2 \nonumber \\
  &=& \sigma^{-2}\sum_{i=1}^n \left[(1-r_{i-1}e_{i-1})x_i^2+2e_{i-1}x_ix_{i-1} 
                 -r_{i-1}e_{i-1}x_{i-1}^2\right] 
\end{eqnarray}
To get to the second form we have shifted the index in the second sum.
Now using the definitions for $d_i$ 
and $e_i$ from Eqs.\ (\ref{eq:5}) and (\ref{eq:6}), we obtain
\bel{9}
Q_n=\sigma^{-2}\sum_{i=1}^{n} {(x_i - r_{i-1} x_{i-1})^2 \over 1-r_{i-1}^2}
\ee
Therefore the joint probability distribution function $P_n$ can be
expressed
\bel{10}
  P_n(x_n,\ldots,x_1)=\left[ \det(2\pi \bfC_n) \right]^{-1/2}
 \exp\left(-\sum_{i=1}^{n} 
     {(x_i - r_{i-1} x_{i-1})^2 \over 2\sigma^2(1-r_{i-1}^2)}\right)
\ee
A special case is for $n=1$, where $r_0=0$ and
\bel{11}
   P_1(x_1) = (2\pi \sigma^2)^{-1/2} \exp\left( -{x_1^2 \over 2\sigma^2} \right)
\ee

\section{Conditional probability functions}

Conditional probabilities can be easily found from the
above joint probability distribution functions.  These are
useful for a number of reasons.  First, the joint
distribution functions can be expressed simply in terms of them
(Markov property).  Second, they can be used as the basis for
constructing fast data simulations via recursion.  Third,
they are necessary for discussion of the random walk process,
for which, as we shall see, the joint distribution becomes singular.

We are particularly
interested in finding the conditional distribution of one of the
$n$ variables, say $x_m$, keeping the others fixed.  
There are two fairly distinct cases, depending on whether $x_m$ is
an ``endpoint'', i.e., either $m=0$ or $m=n$, or an ``interior''
point, i.e., $1 < m < n$.  We shall treat the endpoint case
first and then the interior point case.

\subsection{Endpoint case}

We shall treat the endpoint case for a truncated set of times,
namely, $t_1,\ldots,t_m$, where $m \le n$. 
It is sufficient to find the conditional distribution of $x_m$
given $x_1,\dots,x_{m-1}$, since the other endpoint cases are
completely analogous.  This distribution is given by
\bel{12}
    P_m(x_m | x_{m-1},\ldots,x_2,x_1) =
  {P_m(x_m,x_{m-1},\ldots,x_2,x_1) \over 
     P_{m-1}(x_{m-1},x_{m-2},\ldots,x_2,x_1)}
\ee
Substituting the appropriate expressions for the joint probability
distribution functions on the right (replacing $n$ in Eq.\ (\ref{eq:10})
with $m$ and $m-1$), we find that $P_m(x_m | x_{m-1},\ldots,x_2,x_1)$
is given by
\bel{13'}
P_m(x_m | x_{m-1}) = 
    \left[{\det(2\pi \bfC_{m-1}) \over \det(2\pi \bfC_m)} \right]^{1/2}
  \exp \left[ -{(x_m-r_{m-1}x_{m-1})^2 \over 2\sigma^2 (1-r_{m-1}^2)} \right]
\ee
Since the integral of the condition distribution over $x_m$ must be unity,
we may also write
\bel{13}
P_m(x_m | x_{m-1}) = {1\over \sqrt{2\pi\sigma^2 (1-r_{m-1}^2)}}
  \exp \left[ -{(x_m-r_{m-1}x_{m-1})^2 \over 2\sigma^2 (1-r_{m-1}^2)} \right]
\ee
One sees that this conditional probability is independent of the
variables $x_{m-2}$, $x_{m-3}$, $\ldots$, $x_1$, which is now reflected
in our notation $P_m(x_m | x_{m-1})$.  This is an expression of the
Markov property of the process, so that the distribution of $x_m$ depends only
on the {\em immediately previous} value $x_{m-1}$.

From the conditional probability given by Eq.\ (\ref{eq:13}) it is clear that 
$x_m$ has a Gaussian distribution with mean $r_{m-1}x_{m-1}$
and variance $\sigma^2 (1-r_{m-1}^2)$.  Thus,
for small values of $\alpha(t_m-t_{m-1})$, the mean will be close
to the preceding value $x_{m-1}$ with small variance,
but for larger values there
will be a regression to the (zero) mean of the process, while the
variance will become close to the full variance $\sigma^2$.

Another important deduction can be made from these equations.
The equality of the normalization factors in Eqs.\
(\ref{eq:13'}) and (\ref{eq:13}) implies the relation
\bel{14}
   \det \bfC_m = \sigma^2 (1-r_{m-1}^2)\det \bfC_{m-1}
\ee
Since $\det \bfC_1 = \sigma^2$, we have the general result
\bel{15}
   \det \bfC_m = \sigma^{2m} \prod_{i=1}^{m-1} (1-r_i^2)
\ee
As a consequence, Eq.\ (\ref{eq:10}) can be written
\bel{16}
  P_n(x_n,\ldots,x_1)=\prod_{i=1}^n {1 \over \sqrt{2\pi\sigma^2 (1-r_{i-1}^2)}}
 \exp\left[
     -{(x_i - r_{i-1} x_{i-1})^2 \over 2\sigma^2(1-r_{i-1}^2)}\right]
\ee
or
\bel{17}
  P_n(x_n,\ldots,x_1)=P_1(x_1)\prod_{i=2}^n P_i(x_i | x_{i-1})
\ee
Therefore the complete joint distribution function can be expressed
simply as a product of conditional distibution functions.

\subsection{Interior point case}

Another important type of conditional probability involves not an
endpoint value (like $x_1$ or $x_n$), but an interior value, say $m$.
Then we desire
\bel{23}
   P(x_m | x_n,\ldots,x_{m+1},x_{m-1},\ldots,x_1)
   ={P(x_n,\ldots,x_{m+1},x_m,x_{m-1},\ldots,x_1) \over
     P(x_n,\ldots,x_{m+1},x_{m-1},\ldots,x_1)}
\ee
Using Eq.\ (\ref{eq:17}) the right hand side can be written
\bel{24}
  {P(x_{n} | x_{n-1}) P(x_{n-1} | x_{n-2}) \cdots  P(x_{m+1} | x_{m})
    P(x_{m}| x_{m-1}) \cdots P(x_{2} | x_{1})P(x_1) \over
    P(x_{n} | x_{n-1}) P(x_{n-1} | x_{n-2}) \cdots
   P(x_{m+1} | x_{m-1}) \cdots P(x_{2} | x_{1})P(x_1)}
\ee
which, after cancellations, becomes simply
\bel{25}
   P(x_m | x_n,\ldots,x_{m+1},x_{m-1},\ldots,x_1)
  ={P(x_{m+1} | x_{m}) P(x_{m} | x_{m-1}) \over   P(x_{m+1} | x_{m-1}) }
\ee
We now wish to use Eq.\ (\ref{eq:13}) for the
three conditional probabilities on the right.
However, first we must realize
that the meaning of the quantity  $P(x_{m+1} | x_{m-1})$  
is with respect to the values at $t_{m+1}$ and $t_{m-1}$ with no
intervening values.  Therefore it is necessary to replace the
value for $r_m$ in Eq.\ (\ref{eq:13}) 
with $r_{m-1}r_{m}=\exp[-\alpha(t_{m+1}-t_{m-1})]$ to obtain 
$P(x_{m+1} | x_{m-1})$.  With these substitutions,
\begin{eqnarray}\label{eq:26}
 &\, &  P(x_m | x_n,\ldots,x_{m+1},x_{m-1},\ldots,x_1) 
= \left[{1-r_{m-1}^2 r_{m}^2 \over 2\pi \sigma^2
     (1-r_{m-1}^2)(1-r_{m}^2)}\right]^{1/2} \times \nonumber\\
 &\,& \times \exp\left\lbrace -{1\over 2\pi\sigma^2} \left[
   {(x_{m+1}-r_{m}x_{m})^2 \over 1-r_{m}^2}
  + {(x_{m+}-r_{m-1}x_{m-1})^2 \over 1-r_{m-1}^2} 
- {(x_{m+1}-r_{m-1}r_{m}x_{m-1})^2 \over 1-r_{m-1}^2r_{m}^2} \right]
   \right\rbrace\nonumber \\
\end{eqnarray}
After lengthy, but straightforward, manipulations, this can be expressed,
\bel{27}
   P(x_m | x_{m+1},x_{m-1})
 = {1\over \sqrt{2\pi{\bar \sigma}^2}} \exp\left[-{(x_m-{\bar x})^2
    \over 2 {\bar \sigma}^2}\right]
\ee
where
\begin{eqnarray}\label{eq:28}
  {\bar x}  &=& {r_m(1-r_{m-1}^2)x_{m+1} +r_{m-1}(1-r_m^2)x_{m-1} \over
                   1-r_{m-1}^2r_{m}^2} \nonumber\\
  {\bar \sigma} &=& \sigma \left[ {(1-r_{m-1}^2)(1-r_m^2) \over 
                   1-r_{m-1}^2r_{m}^2} \right]^{1/2}
\end{eqnarray}
This ``interior'' conditional probability depends only on the neighboring
values of the process, namely, $x_{m+1}$ and $x_{m-1}$,
which is reflected in the notation in Eq.\ (\ref{eq:27}).  This is
the appropriate expression of the Markov property for an interior point.
We note that ${\bar x}$ approaches the values $x_{m+1}$ or $x_{m-1}$
as $t_m$ approaches the corresponding $t_{m+1}$ or $t_{m-1}$, but
there is also some regression to the (zero) mean of the process 
between these points.
The variance is zero at the boundary points but reaches a maximum
between the points, which may be close to the full variance
$\sigma^2$ if $\alpha(t_{m+1}-t_{m-1})$ is large.

\section{Case of a Random Walk}

The Gaussian random walk formulas can be derived from the preceding ones
by judicious use of the limit
$\alpha \rightarrow 0$ while maintaining Eq.\ (\ref{eq:18}).  This limit
yields
\bel{19}
 r_{m-1} \rightarrow 1, \qquad \sigma^2(1-r_{m-1}^2) \rightarrow D(t_m-t_{m-1})
\ee
so by Eqs.\ (\ref{eq:13}) and (\ref{eq:14}) the conditional probability 
$P_i(x_m | x_{m-1},\ldots,x_2,x_1)$ becomes
\bel{20}
  P_m(x_m | x_{m-1}) = {1 \over \sqrt{2\pi D(t_m-t_{m-1})}}
                       \exp\left[-{(x_m-x_{m-1})^2 \over 2D(t_m-t_{m-1})}\right]
\ee
Therefore, if the variables $x_{m-1}$, $\ldots$, $x_1$ are known,
then the variable $x_m$ is a Gaussian deviate with mean $x_{m-1}$
and variance $D(t_m-t_{m-1})$.  Note that there is no tendancy now
for the mean to regress to the (zero) mean of the process, since the
succeeding values are all statistically centered on $x_{m-1}$.
The fact that the standard deviation grows as the square root of the time
difference $(t_m-t_{m-1})$ is described by saying that the process
is a ``diffusive'' process, with ``diffusivity'' $D$.
This is often called the {\em random walk} process.



It is possible to find the analogous
interior point formulas for the random walk.
We again have Eq.\ (\ref{eq:27}) but with
\begin{eqnarray}\label{eq:28.1}
  {\bar x}  &=&  {(t_{m+1}-t_{m})x_{m-1} + (t_{m}-t_{m-1})x_{m+1} \over
                   t_{m+1}-t_{m-1}} \nonumber\\
  {\bar \sigma} &=& {D(t_{m+1}-t_{m})(t_{m}-t_{m-1}) \over
                   t_{m+1}-t_{m-1}}
\end{eqnarray}
Thus the mean ${\bar x}$ is given simply by linear interpolation
between the bounding values.  [This also implies that the Wiener
optimum interpolation for a noise-free random walk is simple
linear interpolation.]

\section{Fast Generation of Simulated Gaussian Random Functions with
  Exponential Correlation Functions}

Suppose one wants to simulate data from a Gaussian random process with
an exponential correlation function.
We solve this problem by recursion, generating successive values in the order
$x_1$, $x_2$, $\ldots$, $x_{n-1}$, $x_n$ (a similar method could
be given to generate values in the reverse order).  
The first value $x_1$ at the smallest time can be chosen
according to the probability distribution obtained from Eq.\ (\ref{eq:10}).
That is, $x_1$ can be generated simply by choosing a Gaussian random deviate
with zero mean and with variance $\sigma^2$. 
At any later stage of the process, 
say the $m$th, where $x_{m-1}$, $x_{m-2}$, $\ldots$, $x_2$, $x_1$ are
known, we must generate $x_m$ using the conditional
distribution function (\ref{eq:13}).
Thus one determines the simulated value $x_m$ as a Gaussian random deviate
with mean $r_{m-1}x_{m-1}$ and variance $\sigma^2 (1-r_{m-1}^2)$.
Successive values of the $x_i$ can be determined in this way, leading
to a complete simulation of the process.

For some applications it may be desirable to be able to add new
points to the process that are not necessarily endpoints.  This can
be done through the use of Eqs.\ (\ref{eq:27}) and (\ref{eq:28})
rather than Eq.\ (\ref{eq:13}).  This is slightly more complicated
procedure, but statistically is equivalent.

It is easy to simulate the random walk process using the same recursive
steps as for the exponential process, except that the initial value, say
$x_1$, is not well defined, since the variance of that value is
technically infinite.  However, if one simply takes the initial value
to be $x_1=0$, then all other values can be interpreted as
the increments from that initial value.




\end{document}
