# Kelvin Waves, Klein-Kramers and Kolmogorov Equations, Path-Dependent
Financial Instruments: Survey and New Results

A. Lipton
ADIA, Abu Dhabi, UAE
ADIA Lab, Abu Dhabi, UAE
The Hebrew University of Jerusalem, Jerusalem, Israel
Khalifa University, Abu Dhabi, UAE
MIT, Cambridge, MA, USA Global Head, Research and Development, Strategy and
Planning DepartmentAdvisory Board MemberVisiting Professor and Dean’s Fellow,
Business SchoolProfessor of Practice, Department of MathematicsConnection
Science Fellow, Media Lab

###### Abstract

We discover several surprising relationships between large classes of
seemingly unrelated foundational problems of financial engineering and
fundamental problems of hydrodynamics and molecular physics. Solutions in all
these domains can be reduced to solving affine differential equations commonly
used in various mathematical and scientific disciplines to model dynamic
systems. We have identified connections in these seemingly disparate areas as
we link together small wave-like perturbations of linear flows in ideal and
viscous fluids described in hydrodynamics by Kevin waves to motions of free
and harmonically bound particles described in molecular physics by Klein-
Kramers and Kolmogorov equations to Gaussian and non-Gaussian affine
processes, e.g., Ornstein-Uhlenbeck and Feller, arising in financial
engineering. To further emphasize the parallels between these diverse fields,
we build a coherent mathematical framework using Kevin waves to construct
transition probability density functions for problems in hydrodynamics,
molecular physics, and financial engineering. As one of the outcomes of our
analysis, we discover that the original solution of the Kolmogorov equation
contains an error, which we subsequently correct. We apply our
interdisciplinary approach to advance the understanding of various financial
engineering topics, such as pricing of Asian options, volatility and variance
swaps, options on stocks with path-dependent volatility, bonds, and bond
options. We also discuss further applications to other exciting problems of
financial engineering.

Keywords: affine processes, Kelvin waves, Klein-Kramer equation, Kolmogorov
equation, stochastic volatility, path-dependent options

MSC2020 Classification: 34A34, 35A22, 42A38, 60H10, 76E99, 91G20

###### Contents

  1. 1 Introduction
    1. 1.1 Background
    2. 1.2 Main Results and Paper Structure
  2. 2 Kelvin Waves
    1. 2.1 Euler Equations
    2. 2.2 Navier-Stokes Equations
  3. 3 Klein-Kramers Equation
    1. 3.1 Langevin Equation
    2. 3.2 Klein-Kramers Equation
    3. 3.3 Chandrasekhar’s Solutions
  4. 4 Kolmogorov Equation
    1. 4.1 Summary of Kolmogorov’s paper
    2. 4.2 Challenge and Response
    3. 4.3 Direct Verification
    4. 4.4 Solution via Kelvin Waves
  5. 5 Gaussian Processes
    1. 5.1 The General Case
    2. 5.2 Example: Kolmogorov Process
    3. 5.3 Example: OU Process
      1. 5.3.1 OU Process
      2. 5.3.2 Degenerate Augmented OU Process
    4. 5.4 Example: Evolution of Free and Harmonically Bound Particles
    5. 5.5 Example: Evolution of Vorticity of Two-Dimensional Flows
  6. 6 Non-Gaussian Processes
    1. 6.1 The General Case
    2. 6.2 Example: Feller Process
      1. 6.2.1 Feller Process
      2. 6.2.2 Degenerate Augmented Feller Process
      3. 6.2.3 Non-Degenerate Augmented Feller Process
    3. 6.3 Example: Augmented OU Process
      1. 6.3.1 Degenerate Augmented OU Process, I
      2. 6.3.2 Degenerate Augmented OU Process, II
      3. 6.3.3 Non-Degenerate Augmented OU Process
  7. 7 Financial Instruments
    1. 7.1 Arithmetic and Geometric Brownian Motions and Their Generalizations
    2. 7.2 Asian Options with Geometric Averaging
    3. 7.3 Volatility and Variance Swaps and Swaptions
      1. 7.3.1 Volatility Swaps and Swaptions
      2. 7.3.2 Variance Swaps and Swaptions
    4. 7.4 Path-Dependent Volatility Model
    5. 7.5 Bonds and Bond Options
      1. 7.5.1 Vasicek Model
      2. 7.5.2 CIR Model
    6. 7.6 Further Extensions
  8. 8 Conclusions
  9. A Killed Gaussian Process
  10. B Killed Non-Gaussian Process

##  1 Introduction

###  1.1 Background

The Fourier transform is the most powerful tool of an applied mathematician,
who can use it to solve linear partial differential equations (PDEs) with
spatially constant coefficients, analyze time series and perform many other
vital tasks. The principal building block of the n-dimensional Fourier method
is the wave function:

$\begin{array}[]{c}H\left(t,\mathbf{x,k}\right)=a\left(t\right)\exp\left(i\mathbf{k}\cdot\mathbf{x}\right),\end{array}$
(1)

where $\mathbf{x}$ and $\mathbf{k}$ are $n$-dimensional vectors, $\cdot$
denote the scalar product, $a\left(t\right)$ is the amplitude, and
$\mathbf{k}\cdot\mathbf{x}$ is the phase. Depending on the circumstances, the
amplitude can be a scalar or a vector. Substituting $H$ into a PDE with
spatially constant coefficients, we reduce the problem of interest to a system
of ordinary differential equations (ODEs) or a single ODE when
$a\left(t\right)$ is scalar.

In this paper, we study PDEs with coefficients linearly dependent on $x$,
which we call affine PDEs. Hence, we must use a more general approach and
consider the wave function with a time-dependent wave vector:

$\begin{array}[]{c}H\left(t,\mathbf{x,k}\right)=a\left(t\right)\exp\left(i\mathbf{k}\left(t\right)\cdot\mathbf{x}\right),\end{array}$
(2)

To the best of our knowledge, Kelvin (1887) and Orr (1907) were the first to
use such building blocks to analyze the stability of the steady motions of a
perfect fluid. Subsequently and independently, affine PDEs and the associated
wave functions were used by many researchers in various areas, including the
theory of stochastic processes, physics, biology, and mathematical finance, to
mention a few. In this respect, the Ornstein-Uhlenbeck (OU) and Feller
processes are archetypal; see Uhlenbeck & Ornstein (1930); Chandresekhar
(1943); Feller (1951, 1952). For financial applications of affine processes
see Duffie & Kan (1996); Duffie et al. (2000); Dai & Singleton (2000); Duffie
et al. (2003); Filipovic (2009), among others.

We use Kelvin waves of the form (2) in order to study transition probability
density functions (t.p.d.fs) for affine stochastic processes, which can be
either degenerate, i.e., having more independent components than the sources
of uncertainty, or non-degenerate, when every component has its source of
uncertainty. Our main objective is to link various problems of financial
engineering with their counterparts occurring in fluid dynamics and molecular
physics and showcase the interdisciplinary nature of quantitative finance and
financial modeling. Finding such connections allows us to understand better
how to model, price, and risk-manage various financial instruments and provide
additional intuition regarding their salient features. This work continues our
previous efforts in this direction; see Lipton & Sepp (2008); Lipton (2018),
Chapter 12.

There are several approaches one can use to solve affine equations
efficiently. For instance, Lie symmetries are a powerful tool for studying
certain classes of affine equations. Numerous authors describe general
techniques based on Lie symmetries; see, e.g., Ovsiannikov (1982); Ibragimov
(1985); Olver (1986); Bluman & Kumei (1989), while their specific applications
are covered by Berest (1993); Aksenov (1995); Craddock & Platen (2004);
Craddock (2012); Kovalenko et al. (2014), among many others. Reductions of a
given equation to a simpler, solvable form is another powerful method that can
be successfully used in many instances; see, e.g., Chandresekhar (1943); Carr
et al. (2002); Lipton et al. (2014); Lipton (2018), Chapter 9. Finally, Kelvin
waves (or affine ansatz) provide another approach, which is the focus of the
present paper; see also Duffie & Kan (1996); Dai & Singleton (2000); Duffie et
al. (2003); Lipton & Sepp (2008); Filipovic (2009); Lipton (2018), Chapter 12.

###  1.2 Main Results and Paper Structure

This paper builds a coherent mathematical framework using Kelvin waves as a
powerful and versatile tool for studying t.p.d.fs in the context of generic
affine processes. We discover previously uncovered connections among large
classes of apparently unrelated problems from hydrodynamics, molecular
physics, and financial engineering. All these problems require solving affine
differential equations, i.e., equations with coefficients, which linearly
depend on spatial variables. We discuss some classical results and derive new
ones related to 1) small wave-like perturbations of linear flows of ideal and
viscous fluids described by Euler and Navier-Stokes equations, respectively;
2) motions of free and harmonically bound particles described by the Klein-
Kramers equations and the hypoelliptic Kolmogorov equation; and 3) dynamics of
financial assets and market movements driven by Gaussian and non-Gaussian
affine processes, such as the Ornstein-Uhlenbeck and Feller. We augment these
processes by introducing subordinate processes for auxiliary variables, such
as integrals over the original mean-reverting variable. We develop a uniform
mathematical formalism to construct t.p.d.fs for the abovementioned problems.

Our analysis identifies an error in the original solution of the Kolmogorov
equation, which we subsequently rectify. This work, in addition to recovering
and improving several well-known results, also derives new ones. For instance,
we find a considerably more convenient expression for t.p.d.fs in the case of
processes with stochastic volatility. Moreover, our analysis reveals an
unexpected similarity between the propagation of vorticity in two-dimensional
flows of viscous incompressible fluid and the motion of a harmonically bound
particle, which we use to find a new explicit expression for the vorticity of
a two-dimensional flow in terms of the Gaussian density.

Finally, we apply the new methodology to various financial engineering topics,
such as pricing Asian options, volatility and variance swaps, options on
stocks with path-dependent volatility, and bonds and bond options. In contrast
to the classical treatment, we view fixed-income products as path-dependent,
which allows us to gain additional intuition regarding their pricing and risk
management. We also highlight the flexibility of our interdisciplinary
framework to incorporate additional complexities and to have natural
applications to other financial engineering problems, such as jump-diffusion
processes (verbatim) and, more generally, processes driven by affine pseudo
differential equations.

The paper is organized as follows. Section 2 introduces the Euler equations
describing the evolution of the perfect fluid. We show that these equations
have exact solutions with the velocity linearly and pressure quadratically
depending on spatial coordinates. The linearized Euler equations govern small
perturbations of these solutions known as the celebrated Kelvin waves. We
demonstrate how to use Kelvin waves to study the stability of the exact linear
solutions. Subsequently, we utilize Kevin waves as our primary tool in our
methodology and apply it to other areas of study. We also briefly cover the
Navier-Stokes equation for the viscous fluid and show how to modify the ODEs
for Kelvin waves to account for viscosity. In Section 5.5, we return to this
topic, establishing deep connections between the two-dimensional vorticity
equations and the Klein-Kramers equation underpinning molecular physics, which
allows us to derive a new formula expressing the vorticity as a Gaussian
density, and the stream function as a solution of the corresponding Poisson
equation.

Section 3 deals with some representative stochastic problems occurring in
physics. Namely, Section 3.1 introduces the Langevin equation describing the
motion of the underdamped Brownian particle. The Klein-Kramers equation
describing the motion of such particle is derived in Section 3.2. Its
solutions found by Chandrasekhar for a free and a harmonic particle are given
in Section 3.3. The Klein-Kramers equation is naturally degenerate since there
is a source of white noise for the particle velocity but not for its position.
As we shall see later, many path-dependent mathematical finance problems fall
into this category, provided the corresponding financial variables, such as
the geometric price averages naturally occurring as underlyings for Asian
options, are interpreted as path integrals. Section 4 studies the celebrated
degenerate stochastic process introduced by Kolmogoroff (1934), which is
closely linked with the Klein-Kramers equation studied in 3, the associated
Fokker-Planck equation, and its solution proposed by Kolmogorov. Section 4.1
summarizes original Kolmogorov’s results. In Section 4.2, to our surprise, we
show that the Fokker-Planck equation used by Kolmogorov in his original paper
is inconsistent with his assumptions about the underlying process. Moreover,
his proposed solution lacks dimensional consistency and must be rectified to
satisfy the Fokker-Planck equation and the initial condition. Thankfully, not
all is lost, and Kolmogorov’s solution can be corrected, as we demonstrate in
Sections 4.3, 4.4. We start with direct calculation in Section 4.3 and then
arrive at the same corrected Kolmogorov’s solution by using Kelvin waves in
Section 4.4.

Sections 5, 6 constitute the heart of the paper. Section 5 performs a detailed
analysis of t.p.d.fs for Gaussian affine processes, whose drifts are linearly
dependent on the spatial coordinates, while covariance matrices are spatially
(but not necessarily temporary) invariant. Section 6 investigates t.p.d.fs for
non-Gaussian processes with drift and covariance matrix linearly dependent on
spatial coordinates. Our methodology demonstrates that Kelvin waves are
ideally suited for studying such affine processes as the OU and Feller
processes that often occur in finance, asset management, and trading.

Section 5 applies our methodology to Gaussian processes. Section 5.1 derives
the general formula for the Gaussian dynamics, which covers both degenerate
(like in Kolmogorov’s case) and non-degenerate situations. In Section 5.2, we
derive the t.p.d.f. for the Kolmogorov process with time-dependent
coefficients. Section 5.3, investigates the OU process with time-dependent
coefficients and the augmented OU process, which describes the joint dynamics
of the process and its integral. While the results are classical, their
derivation based on Kelvin wave expansions complements them naturally and adds
value to the analysis by providing a different perspective to understand and
derive these well-established results. Section 5.4 investigates the case of
free and harmonically bound particles and compares and contrasts the Kelvin
wave approach with the classical one developed by Chandresekhar (1943).
Finally, Section 5.5 goes back to our departure point described in Section 2.
It shows that the temporal-spatial evolution of the two-dimensional flow’s
vorticity for a viscous fluid is very similar to the dynamics of a
harmonically bound particle. This is a notable and novel result that bridges
seemingly disparate underlying physical mechanisms.

Non-Gaussian processes are covered in Section 6. Section 6.1 derives the
general formula for the non-Gaussian dynamics, which covers both degenerate
and non-degenerate possibilities. Section 6.2 treats the Feller process with
constant coefficients and two augmented Feller processes \- a degenerate
linking the Feller process with its integral and a non-degenerate describing
an arithmetic Brownian motion with stochastic volatility driven by Feller
process. As a notable byproduct of our analysis, we develop an alternative
approach for analyzing finite-time explosions for exponential moments for
t.p.d.fs of augmented Feller processes. Finally, Section 6.3 applies our
methodology to degenerate and non-degenerate arithmetic Brownian motions with
stochastic volatility driven by an OU process.

Section 7 extends our framework to the realm of financial applications.
Section 7.1 briefly outlines the basics of financial engineering based on the
geometric Brownian motion paradigm. Since this paradigm is too narrow to be
practical, we introduce modifications to better capture the complexities and
nuances of financial markets. It turns out that our new methodology is
versatile enough to analyze numerous practically important cases. In Section
7.2, we describe how Asian options with geometric averaging can be priced by
using the corrected version of the Kolmogorov solution. Section 7.3 applies
the results of Sections 5.3 and 6.2 to the pricing of volatility and variance
swaps and swaptions, respectively. In Section 7.4, we introduce a path-
dependent volatility model and show how to use an augmented OU process
discussed in Section 6.3 in the context of moving particles to construct a
novel analytically solvable model. Section 7.5 expands our methodology to
cover the Vasicek-Hull-White and Cox-Ingersoll-Ross (CIR) frameworks for bonds
and bond options. In contrast to the conventional treatment based on solving
the backward Kolmogorov equation for the killed process, we prefer to treat
bonds and bond options as path-dependent instruments. This approach enables us
to utilize formulas derived in Sections 5, 6 to price fixed-income products.
These formulas give us additional insights into pricing and risk management of
bonds and bond options and connect the corresponding pricing formulas with the
ones derived in physics for describing the motion of underdamped Brownian
particles. Section 7.6 briefly discusses the potential applications of Kelvin
waves to other financial engineering problems.

Section 8 contains our conclusions. Appendices A, B, discuss killed affine
processes.

A word about notation is in order. By its very nature, papers like the present
one cannot use the unified notation. The reason is simple - the ground is too
broad and stretches from fluid dynamics to molecular physics, probability
theory, and financial engineering. Each of these fields uses its own set of
notations, which are not coherent across different areas. The best we could do
is to use consistent notation in each Section and, as much as possible, across
sections, but the reader must be alert to keep her train of thought.

##  2 Kelvin Waves

###  2.1 Euler Equations

In what follows, we use Kelvin waves as our primary tool. We dedicate this
section to their brief description.

Recall that the Euler system of partial differential equations (PDEs)
describing the motion of an inviscid, incompressible fluid has the form:

$\begin{array}[]{c}\frac{\partial\mathbf{V}}{\partial
t}+(\mathbf{V\cdot}\nabla)\mathbf{V+}\nabla\left(\frac{P}{\rho}\right)=0,\\\
\\\ \nabla\cdot\mathbf{V=}0\mathbf{,}\end{array}$ (3)

where $t$ is time, $\mathbf{x}$ is the position,
$V$$\left(t,\mathbf{x}\right)$ is the velocity vector,
$P\left(t,\mathbf{x}\right)$ is the pressure, $\rho$ is the constant density,
$\nabla$ is the gradient; see, e.g., Chandrasekhar (1961). In Cartesian
coordinates, Eqs (3) can be written as follows:

$\begin{array}[]{c}\frac{\partial V_{i}}{\partial
t}+\mathop{\displaystyle\sum}\limits_{j=1}^{3}V_{j}\frac{\partial
V_{i}}{\partial x_{j}}+\frac{\partial}{\partial
x_{i}}\left(\frac{P}{\rho}\right)=0,\\\ \\\
\mathop{\displaystyle\sum}\limits_{i=1}^{3}\frac{\partial V_{i}}{\partial
x_{i}}=0.\end{array}$ (4)

It can be shown that Eqs (3) have a family of solutions
$\left(\mathbf{V}\left(t,\mathbf{x}\right),P\left(t,\mathbf{x}\right)\right)$,
linearly depending on spatial coordinates:

$\begin{array}[]{c}\mathbf{V}\left(t,\mathbf{x}\right)=\mathfrak{L}\left(t\right)\mathbf{x,\
\ \
}\frac{P\left(t,\mathbf{x}\right)}{\rho}=\frac{1}{2}\mathfrak{M}\left(t\right)\mathbf{x\cdot
x,}\end{array}$ (5)

where the matrices $\mathfrak{L}\left(t\right)$, $\mathfrak{M}\left(t\right)$,
such that

$\begin{array}[]{c}Tr\left(\mathfrak{L}\left(t\right)\right)=0,\ \ \
\mathfrak{M}\left(t\right)=\mathfrak{M}^{\ast}\left(t\right),\end{array}$ (6)

satisfy the following system of ordinary differential equations (ODEs):

$\begin{array}[]{c}\frac{d\mathfrak{L}\left(t\right)}{dt}+\mathfrak{L}^{2}\left(t\right)+\mathfrak{M}\left(t\right)=0.\end{array}$
(7)

We are interested in the behavior of perturbations of solutions given by Eq.
(5), which we denote by
$\left(\mathbf{v}\left(t,\mathbf{x}\right),p\left(t,\mathbf{x}\right)\right)$.
By neglecting the quadratic term $($$v$$\mathbf{\cdot}\nabla)$$v$, we can
write the system of PDEs for $\left(\mathbf{v,}p\right)$ as follows:

$\begin{array}[]{c}\frac{\partial\mathbf{v}}{\partial
t}+(\mathfrak{L}\left(t\right)\mathbf{x\cdot}\nabla)\mathbf{v+}\mathfrak{L}\left(t\right)\mathbf{v+}\nabla\left(\frac{p}{\rho}\right)=0,\\\
\\\ \nabla\cdot\mathbf{v=}0\mathbf{.}\end{array}$ (8)

It has been known for a long time that linear PDEs (8) have wave-like
solutions of the form:

$\begin{array}[]{c}\left(\mathbf{v}\left(t,\mathbf{x}\right),\frac{p\left(t,\mathbf{x}\right)}{\rho}\right)\\\
\\\
=\left(\mathbf{a}\left(t\right)\exp\left(i\mathbf{\Gamma}\left(t\right)\cdot\left(\mathbf{x-\xi}\left(t\right)\right)\right),b\left(t\right)\exp\left(i\mathbf{\Gamma}\left(t\right)\cdot\left(\mathbf{x-\xi}\left(t\right)\right)\right)\right),\end{array}$
(9)

where $\left(\mathbf{a}\left(t\right),b\left(t\right)\right)$ are time-
dependent amplitudes, and $\mathbf{\Gamma}\left(t\right)$ is the time-
dependent wave vector; see Kelvin (1887); Orr (1907); Craik & Criminale
(1986); Friedlander & Lipton-Lifschitz (2003). We call these solutions the
Kelvin waves. The point we emphasize below is that the so-called affine ansatz
is a special instance of Kelvin wave. This observation allows us to draw
similarities among seemingly unrelated topics, which, in turn, facilitates
their holistic and comprehensive study.

The pair $\mathbf{a}\left(t\right),\mathbf{\Gamma}\left(t\right)$ satisfies
the following system of ODEs:

$\begin{array}[]{c}\frac{d\mathbf{\xi}\left(t\right)}{dt}-\mathfrak{L}\left(t\right)\mathbf{\xi}\left(t\right)=0,\
\ \ \mathbf{\xi}\left(0\right)=\mathbf{\xi}_{0},\\\ \\\
\frac{d\mathbf{\Gamma}\left(t\right)}{dt}+\mathfrak{L}^{T}\left(t\right)\mathbf{\Gamma}\left(t\right)=0,\
\ \ \mathbf{\Gamma}\left(0\right)=\mathbf{\Gamma}_{0},\\\ \\\
\frac{d\mathbf{a}\left(t\right)}{dt}+\mathfrak{L}\left(t\right)\mathbf{a}\left(t\right)-2\frac{\mathfrak{L}\left(t\right)\mathbf{a}\left(t\right)\cdot\mathbf{k}\left(t\right)}{\mathbf{k}\left(t\right)\cdot\mathbf{k}\left(t\right)}\mathbf{k}\left(t\right)=0,\
\mathbf{a}\left(0\right)=\mathbf{a}_{0},\\\ \\\ \
\mathbf{\Gamma}_{0}\cdot\mathbf{a}_{0}=0.\end{array}$ (10)

Here and below, the superscript $T$ stands for transpose, and $\cdot$ denotes
the scalar product. The corresponding $p\left(t\right)$ can be found via the
incompressibility condition. It is easy to show that for $t\geq 0$,

$\begin{array}[]{c}\mathbf{\Gamma}\left(t\right)\cdot\mathbf{\xi}\left(t\right)=\mathbf{\Gamma}_{0}\cdot\mathbf{\xi}_{0},\
\ \ \mathbf{\Gamma}\left(t\right)\cdot\mathbf{a}\left(t\right)=0.\end{array}$
(11)

Typically, Eqs (10) are used to study the stability of the linear flow. Such a
flow is unstable whenever
$\left\|\mathbf{a}\left(t\right)\right\|\rightarrow\theta$ for some choices of
$\mathbf{\Gamma}_{0},\mathbf{a}_{0}$; see Bayly (1986); Bayly et al. (1996);
Lifschitz (1995). Moreover, it can be shown that the same instabilities occur
in general three-dimensional flows; see Friedlander & Vishik (1991); Lifschitz
& Hameiri (1991); Friedlander & Lipton-Lifschitz (2003).

Interestingly, Chandrasekhar (1961) pointed out that the superposition of the
linear flow (5) and the Kelvin wave (9), i.e.,

$\begin{array}[]{c}\mathbf{\tilde{V}}\left(t,\mathbf{x}\right)=\mathfrak{L}\left(t\right)\mathbf{x+v}\left(t,\mathbf{x}\right),\\\
\\\
\frac{\tilde{P}\left(t,\mathbf{x}\right)}{\rho}=\frac{1}{2}\mathfrak{M}\left(t\right)\mathbf{x\cdot
x+}\frac{p\left(t,\mathbf{x}\right)}{\rho},\end{array}$ (12)

satisfies the nonlinear Euler equations (3) since the nonlinear term
$v$$\cdot\nabla$$v$ vanishes identically due to incompressibility.111Thus,
even the greatest minds occasionally can be myopic – it took eighty years for
fluid dynamists to connect the dots and observe that
$\left(\mathbf{\tilde{v}},\tilde{p}/\rho\right)$ solve the nonlinear Euler
equations.

###  2.2 Navier-Stokes Equations

The motion of the incompressible viscous fluid is described by the celebrated
Navier-Stokes equations of the form:

$\begin{array}[]{c}\frac{\partial\mathbf{V}}{\partial
t}+(\mathbf{V\cdot}\nabla)\mathbf{V-}\nu\Delta\mathbf{V+}\nabla\left(\frac{P}{\rho}\right)=0,\\\
\\\ \nabla\cdot\mathbf{V=}0\mathbf{,}\end{array}$ (13)

where $\nu$ is the kinematic viscosity; see, e.g., Chandrasekhar (1961). It is
clear that linear flows, given by Eq. (5), are unaffected by viscosity, hence
they satisfy Eqs (14). Viscosity does affect their small perturbations, which
are governed by the following equations:

$\begin{array}[]{c}\frac{\partial\mathbf{v}}{\partial
t}+(\mathfrak{L}\left(t\right)\mathbf{x\cdot}\nabla)\mathbf{v+}\mathfrak{L}\left(t\right)\mathbf{v\mathbf{-}}\nu\Delta\mathbf{v+}\nabla\left(\frac{p}{\rho}\right)=0,\\\
\\\ \nabla\cdot\mathbf{v=}0\mathbf{.}\end{array}$ (14)

The viscous version of Eqs (10) has the form; see Lifschitz (1991):

$\begin{array}[]{c}\frac{d\mathbf{\xi}\left(t\right)}{dt}-\mathfrak{L}\left(t\right)\mathbf{\xi}\left(t\right)=0,\
\ \ \mathbf{\xi}\left(0\right)=\mathbf{\xi}_{0},\\\ \\\
\frac{d\mathbf{\Gamma}\left(t\right)}{dt}+\mathfrak{L}^{T}\left(t\right)\mathbf{\Gamma}\left(t\right)=0,\
\ \ \mathbf{\Gamma}\left(0\right)=\mathbf{\Gamma}_{0},\\\ \\\
\frac{d\mathbf{a}\left(t\right)}{dt}+\mathfrak{L}\left(t\right)\mathbf{a}\left(t\right)-2\frac{\mathfrak{L}\left(t\right)\mathbf{a}\left(t\right)\cdot\mathbf{k}\left(t\right)}{\mathbf{k}\left(t\right)\cdot\mathbf{k}\left(t\right)}\mathbf{k}\left(t\right)+\nu\left|\mathbf{k}\left(t\right)\right|^{2}\mathbf{a}\left(t\right)=0,\
\mathbf{a}\left(0\right)=\mathbf{a}_{0},\\\ \\\ \
\mathbf{\Gamma}_{0}\cdot\mathbf{a}_{0}=0.\end{array}$ (15)

We shall show in Section 5.5 that in the two-dimensional case, the Navier-
Stokes equations for small perturbations of linear flows are more or less
identical to the Fokker-Planck equations for harmonically bound articles,
which is surprising.

##  3 Klein-Kramers Equation

###  3.1 Langevin Equation

We start with the Langevin equation for particles moving in a potential field
and impacted by random forces; see Langevin (1908). This section uses the same
notation as in Chandresekhar (1943).

Consider an underdamped Brownian particle. In contrast to the standard
Brownian motion, which is overdamped, we assume that the frictions are finite,
so that we must treat the momentum as an independent degree of freedom, and
describe a particle’s state as a pair $\left(\mathbf{r},\mathbf{u}\right)$,
where $\mathbf{r}$ and $\mathbf{u}$ are its position and velocity,
respectively. We consider a $d$-dimensional space, with $d=1$ and $d=3$ of
particular interest. Hence, the corresponding Langevin equations are

$\begin{array}[]{c}\frac{d\mathbf{r}}{dt}=\mathbf{u},\\\ \\\
\frac{d\mathbf{u}}{dt}=-\beta\mathbf{u-}\frac{\nabla
V\left(\mathbf{r}\right)}{m}+\sqrt{\frac{2\beta
k_{B}T}{m}}\mathbf{\eta}\left(t\right),\\\ \\\
\left\langle\mathbf{\eta}^{T}\left(t\right)\mathbf{\eta}\left(t^{\prime}\right)\right\rangle=\mathbf{I}\delta\left(t-t^{\prime}\right).\end{array}$
(16)

Here $m$ is the particle mass, $\beta$ is the friction coefficient, $k_{B}$ is
the Boltzmann constant, $T$ is the temperature, $V\left(\mathbf{r}\right)$ is
the external potential, $\mathbf{\eta}\left(t\right)$ is a $d$-dimensional
Gaussian white noise, and $\mathbf{I}$ is the $d$-dimensional identity matrix.
Below we denote the ratio $\beta k_{B}T/m$ as $q$.

Of course, we can rewrite Eqs (16) as a system of stochastic differential
equations (SDEs):

$\begin{array}[]{c}d\mathbf{r}=\mathbf{u}dt,\\\ \\\
d\mathbf{u}=-\beta\mathbf{u}dt\mathbf{-}\frac{\nabla
V\left(\mathbf{r}\right)}{m}dt+\sqrt{2q}d\mathbf{W}\left(t\right),\end{array}$
(17)

where $\mathbf{W}\left(t\right)$ is a standard $d$-dimensional Wiener process.

###  3.2 Klein-Kramers Equation

Fokker, Planck, and their numerous followers derived and studied the forward
parabolic equation for the transition probability density
$W\left(\tau,\mathbf{\rho,\upsilon},t,\mathbf{r,u}\right)$ associated with a
stochastic process. For the stochastic process governed by SDEs (17), the
corresponding equation, called the Klein-Kramers equation, has the form:

$\begin{array}[]{c}W_{t}-qW_{\mathbf{uu}}+\mathbf{u}W_{\mathbf{r}}-\left(\left(\beta\mathbf{u}+\frac{\nabla
V\left(\mathbf{r}\right)}{m}\right)W\right)_{\mathbf{u}}=0,\\\ \\\
W\left(\tau,\mathbf{\rho,\upsilon},\tau,\mathbf{r,u}\right)=\delta\left(\mathbf{r}-\mathbf{\rho}\right)\delta\left(\mathbf{u}-\mathbf{\upsilon}\right).\end{array}$
(18)

The backward parabolic Kolmogorov equation has the form:

$\begin{array}[]{c}W_{\tau}+qW_{\mathbf{\upsilon\upsilon}}+\mathbf{\upsilon}W_{\mathbf{\rho}}-\left(\beta\mathbf{\upsilon}+\frac{\nabla
V\left(\mathbf{\rho}\right)}{m}\right)W_{\mathbf{\upsilon}}=0,\\\ \\\
W\left(t,\mathbf{\rho,\upsilon},t,\mathbf{r,u}\right)=\delta\left(\mathbf{r}-\mathbf{\rho}\right)\delta\left(\mathbf{u}-\mathbf{\upsilon}\right).\end{array}$
(19)

Details are given in Fokker (1914); Planck (1917); Klein (1921); Chapman
(1928); Kolmogoroff (1931, 1933, 1934); Kramers (1940); Chandresekhar (1943);
Risken (1989); Hänggi et al. (1990), as well as a multitude of subsequent
sources. For fascinating historical details see Ebeling et al. (2008).

###  3.3 Chandrasekhar’s Solutions

In a well-known survey article, Chandresekhar (1943) described elegant
solutions of Eq. (18) for a free particle and a harmonically bound particle,
which he derived by using clever changes of coordinates.

First, we shall consider a free particle in $1D$. For such a particle,
Chandresekhar (1943) writes the corresponding Klein-Kramers equation as
follows:

$\begin{array}[]{c}W_{t}-qW_{uu}+uW_{r}-\beta uW_{u}-W=0,\\\ \\\
W\left(\tau,\rho\mathbf{,}\upsilon,\tau,r\mathbf{,}u\right)=\delta\left(r-\rho\right)\delta\left(u-\upsilon\right),\end{array}$
(20)

and shows that

$\begin{array}[]{c}W=\frac{1}{2\pi\left(FG-H^{2}\right)^{1/2}}\exp\left(-\frac{\left(FS^{2}-2HRS+GR^{2}\right)}{2\left(FG-H^{2}\right)}\right),\end{array}$
(21)

where

$\begin{array}[]{c}S=u-e^{-\beta\left(t-\tau\right)}\upsilon,\ \ \
R=r-\rho-\frac{1-e^{-\beta\left(t-\tau\right)}}{\beta}\upsilon,\\\ \\\
F=\frac{q}{\beta^{3}}\left(-3+4e^{-\beta\left(t-\tau\right)}-e^{-2\beta\left(t-\tau\right)}+2\beta\left(t-\tau\right)\right),\\\
\\\ G=\frac{q}{\beta}\left(1-e^{-2\beta\left(t-\tau\right)}\right),\ \ \
H=\frac{q}{\beta^{2}}\left(1-e^{-\beta\left(t-\tau\right)}\right)^{2}.\end{array}$
(22)

Since we assumed that stochastic drivers are uncorrelated, the t.p.d.f.
$W^{\left(3\right)}$ can be presented as a product of three $1D$ t.p.d.f.
$W^{\left(1\right)}$:

$\begin{array}[]{c}W^{\left(3\right)}=W_{1}^{\left(1\right)}W_{2}^{\left(1\right)}W_{3}^{\left(1\right)}=\frac{1}{8\pi^{3}\left(FG-H^{2}\right)^{3/2}}\exp\left(-\frac{\left(F\left|\mathbf{S}\right|^{2}-2H\mathbf{R}\cdot\mathbf{S}+G\left|\mathbf{R}\right|^{2}\right)}{2\left(FG-H^{2}\right)}\right),\end{array}$
(23)

where $\mathbf{u=}\left(u_{1},u_{2},u_{3}\right)$, etc.

To analyze a harmonically bound particle, we can proceed as follows; see
Chandresekhar (1943). As usual, the harmonic potential $V\left(x\right)$ has
the form:

$\begin{array}[]{c}V\left(x\right)=\frac{m\omega^{2}x^{2}}{2}.\end{array}$
(24)

Let $\mu_{\pm}$ be the roots of

$\begin{array}[]{c}\mu^{2}+\beta\mu+\omega^{2}=0,\\\ \\\
\mu_{\pm}=\frac{-\beta\pm\sqrt{\beta^{2}-4\omega^{2}}}{2}.\end{array}$ (25)

Define new variables

$\begin{array}[]{c}\xi=e^{-\mu_{-}\left(t-\tau\right)}\left(\mu_{+}r-u\right),\
\ \ \eta=e^{-\mu_{+}\left(t-\tau\right)}\left(\mu_{-}r-u\right),\\\ \\\
\hat{\xi}=\left(\mu_{+}\rho-\upsilon\right),\ \ \
\hat{\eta}=\left(\mu_{-}\rho-\upsilon\right),\end{array}$ (26)

and coefficients

$\begin{array}[]{c}a=\frac{q}{\mu_{+}}\left(1-e^{-2\mu_{+}\left(t-\tau\right)}\right),\
\ \ b=\frac{q}{\mu_{-}}\left(1-e^{-2\mu_{-}\left(t-\tau\right)}\right),\\\ \\\
h=-\frac{2q}{\left(\mu_{+}+\mu_{-}\right)}\left(1-e^{-\left(\mu_{+}+\mu_{-}\right)\left(t-\tau\right)}\right).\end{array}$
(27)

Then

$\begin{array}[]{c}W=\frac{e^{\beta\left(t-\tau\right)}}{2\pi\sqrt{ab-h^{2}}}\exp\left(-\frac{a\left(\xi-\hat{\xi}\right)^{2}+2h\left(\xi-\hat{\xi}\right)\left(\eta-\hat{\eta}\right)+b\left(\eta-\hat{\eta}\right)^{2}}{2\left(ab-h^{2}\right)}\right).\end{array}$
(28)

While Chandresekhar (1943) stopped at Eq. (28), for practical applications, it
is more useful to represent the exponent as an explicit quadratic form of $r$
and $u$, which we do below.

##  4 Kolmogorov Equation

###  4.1 Summary of Kolmogorov’s paper

In a remarkable (and remarkably concise) note, Kolmogorov considers a system
of particles in $n$-dimensional space with coordinates $q_{1},...,q_{n}$, and
velocities $\dot{q}_{1},...,\dot{q}_{n}$, assumes the probability density
function

$\begin{array}[]{c}g\left(t,q_{1},...,q_{n},\dot{q}_{1},...,\dot{q}_{n},t^{\prime},q_{1}^{\prime},...,q_{n}^{\prime},\dot{q}_{1}^{\prime},...,\dot{q}_{n}^{\prime}\right)\end{array}$

exist for some time $t^{\prime}>t$, and reveals (without any explanation) an
analytical expression for $g$ in the one-dimensional case; see Kolmogoroff
(1934). This note is the third in a series of papers; the previous two being
Kolmogoroff (1931, 1933).

Kolmogorov makes the following natural assumptions:

$\begin{array}[]{c}\mathbf{E}\left|\Delta q_{i}-\dot{q}_{i}\Delta
t\right|=o\left(\Delta t\right),\end{array}$ (29)
$\begin{array}[]{c}\mathbf{E}\left(\Delta q_{i}\right)^{2}=o\left(\Delta
t\right),\end{array}$ (30)

where $\Delta t=t^{\prime}-t$. Eqs (29), (30) imply

$\begin{array}[]{c}\mathbf{E}\left(\Delta q_{i}\right)=\dot{q}_{i}\Delta
t+o\left(\Delta t\right),\end{array}$ (31)
$\begin{array}[]{c}\mathbf{E}\left(\Delta q_{i}\Delta
q_{j}\right)\leq\sqrt{\mathbf{E}\left(\Delta
q_{i}\right)^{2}\mathbf{E}\left(\Delta q_{j}\right)^{2}}=o\left(\Delta
t\right).\end{array}$ (32)

Furthermore, under very general assumptions, the following relationships hold

$\begin{array}[]{c}\mathbf{E}\left(\Delta\dot{q}_{i}\right)=f_{i}\left(t,q,\dot{q}\right)\Delta
t+o\left(\Delta t\right),\end{array}$ (33)
$\begin{array}[]{c}\mathbf{E}\left(\Delta\dot{q}_{i}\right)^{2}=k_{ii}\left(t,q,\dot{q}\right)\Delta
t+o\left(\Delta t\right),\end{array}$ (34)
$\begin{array}[]{c}\mathbf{E}\left(\Delta\dot{q}_{i}\Delta\dot{q}_{j}\right)=k_{ij}\left(t,q,\dot{q}\right)\Delta
t+o\left(\Delta t\right),\end{array}$ (35)

where $f,k$ are continuous functions. Eqs (30), (34) imply

$\begin{array}[]{c}\mathbf{E}\left(\Delta\dot{q}_{i}\Delta\dot{q}_{j}\right)\leq\sqrt{\mathbf{E}\left(\Delta\dot{q}_{i}\right)^{2}\mathbf{E}\left(\Delta\dot{q}_{j}\right)^{2}}=o\left(\Delta
t\right).\end{array}$ (36)

Under some natural physical assumptions, it follows that $g$ satisfies the
following differential equation of the Fokker-Planck type:

$\begin{array}[]{c}\frac{\partial g}{\partial
t^{\prime}}=-\sum\dot{q}_{i}^{\prime}\frac{\partial}{\partial
q_{i}^{\prime}}g-\sum\frac{\partial}{\partial\dot{q}_{i}^{\prime}}\left\\{f_{i}\left(t,q,\dot{q}\right)g\right\\}+\sum\sum\frac{\partial^{2}}{\partial\dot{q}_{i}^{\prime
2}}\left\\{\left\\{k\left(t,q,\dot{q}\right)g\right\\}\right\\}.\end{array}$
(37)

In the one-dimensional case, we have

$\begin{array}[]{c}\frac{\partial g}{\partial
t^{\prime}}=-\dot{q}^{\prime}\frac{\partial}{\partial
q^{\prime}}g-\frac{\partial}{\partial\dot{q}^{\prime}}\left\\{f\left(t,q,\dot{q}\right)g\right\\}+\frac{\partial^{2}}{\partial\dot{q}^{\prime
2}}\left\\{\left\\{k\left(t,q,\dot{q}\right)g\right\\}\right\\}.\end{array}$
(38)

These equations are known as ultra-parabolic Fokker-Plank- Kolmogorov
equations due to their degeneracy.

When $f$ and $k$ are constants, the fundamental solution of Eq. (38) has the
form:

$\begin{array}[]{c}g=\frac{2\sqrt{3}}{\pi
k^{2}\left(t^{\prime}-t\right)^{2}}\exp\left\\{-\frac{\left[\dot{q}^{\prime}-\dot{q}-f\left(t^{\prime}-t\right)\right]^{2}}{4k\left(t^{\prime}-t\right)}-\frac{3\left[q^{\prime}-q-\frac{\dot{q}^{\prime}+\dot{q}}{2}\left(t^{\prime}-t\right)\right]^{2}}{k^{3}\left(t^{\prime}-t\right)^{3}}\right\\}.\end{array}$
(39)

One can see that $\Delta\dot{q}$ is of the order $\left(\Delta
t\right)^{1/2}$. At the same time

$\begin{array}[]{c}\Delta q=\dot{q}\Delta t+O\left(\Delta
t\right)^{3/2}.\end{array}$ (40)

One can prove that a similar relation holds for the general Eq. (37).

The previous Section shows that physicists derived Eqs (37), (38) at least a
decade earlier than Kolmogorov, which puts the Kolmogorov equation into a
broader context.

Kolmogorov equations fascinated mathematicians for a long time and generated a
great deal of research; see, e.g., Weber (1951); Hörmander (1967); Kuptsov
(1972); Lanconelli et al. (2002); Pascucci (2005); Ivasishen & Medynsky
(2010); Duong & Tran (2018), among others.

###  4.2 Challenge and Response

Despite its undoubting brilliance, the original Kolmogorov’s paper has several
issues.

First, Eqs (37), (38) _are not_ the Fokker-Planck equations associated with
the process described by Eqs (33) - (35), since they miss the prefactor $1/2$
in front of the diffusion terms. The corrected multivariate equation has the
form:

$\begin{array}[]{c}\frac{\partial g}{\partial
t^{\prime}}=-\sum\dot{q}_{i}^{\prime}\frac{\partial}{\partial
q_{i}^{\prime}}g-\sum\frac{\partial}{\partial\dot{q}_{i}^{\prime}}\left\\{f_{i}\left(t,q,\dot{q}\right)g\right\\}+\frac{1}{2}\sum\sum\frac{\partial^{2}}{\partial\dot{q}_{i}^{\prime}\partial\dot{q}_{j}^{\prime}}\left\\{\left\\{k\left(t,q,\dot{q}\right)g\right\\}\right\\},\end{array}$
(41)

while the corresponding one-dimensional equation has the form:

$\begin{array}[]{c}\frac{\partial g}{\partial
t^{\prime}}=-\dot{q}^{\prime}\frac{\partial}{\partial
q^{\prime}}g-\frac{\partial}{\partial\dot{q}^{\prime}}\left\\{f\left(t,q,\dot{q}\right)g\right\\}+\frac{1}{2}\frac{\partial^{2}}{\partial\dot{q}^{\prime
2}}\left\\{\left\\{k\left(t,q,\dot{q}\right)g\right\\}\right\\}.\end{array}$
(42)

Alternatively, Eqs (34), (35) can be altered as follows:

$\begin{array}[]{c}\mathbf{E}\left(\Delta\dot{q}_{i}\right)^{2}=2k_{ii}\left(t,q,\dot{q}\right)\Delta
t+o\left(\Delta t\right),\end{array}$ (43)
$\begin{array}[]{c}\mathbf{E}\left(\Delta\dot{q}_{i}\Delta\dot{q}_{j}\right)=2k_{ij}\left(t,q,\dot{q}\right)\Delta
t+o\left(\Delta t\right).\end{array}$ (44)

We prefer to update the Fokker-Planck equation.

Second, $g$ given by Eq. (39) does not solve Eq. (38). It also does not
satisfy the (implicit) initial condition

$\begin{array}[]{c}g\left(t,q,\dot{q},t,q^{\prime},\dot{q}^{\prime}\right)=\delta\left(q^{\prime}-q\right)\delta\left(\dot{q}^{\prime}-\dot{q}\right),\end{array}$
(45)

where $\delta\left(.\right)$ is the Dirac $\delta$-function. The fact that
expression (39) does not solve Eq. (38) can be verified by substitution.
However, it is easier to verify this statement via dimensional analysis. The
dimensions of the corresponding variables and coefficients are as follows:

$\begin{array}[]{c}\left[t\right]=\left[t^{\prime}\right]=T,\ \ \
\left[q\right]=\left[q^{\prime}\right]=L,\ \ \
\left[\dot{q}\right]=\left[\dot{q}^{\prime}\right]=T^{-1}L,\ \ \
\left[g\right]=T^{1}L^{-2},\\\ \\\ \left[f\right]=T^{-2}L,\ \ \
\left[k\right]=T^{-3}L^{2}.\end{array}$ (46)

It is easy to show that $g$ is scale-invariant, so that

$\begin{array}[]{c}g\left(\lambda^{2}t,\lambda^{3}q,\lambda\dot{q},\lambda^{2}t^{\prime},\lambda^{3}q^{\prime},\lambda\dot{q}^{\prime};\lambda^{-1}f,k\right)=\lambda^{-4}g\left(t,q,\dot{q},t^{\prime},q^{\prime},\dot{q}^{\prime};f,k\right).\end{array}$
(47)

The original Kolmogorov formula contains two typos, making it dimensionally
incorrect since the term

$\frac{3\left[q^{\prime}-q-\frac{\dot{q}^{\prime}+\dot{q}}{2}\left(t^{\prime}-t\right)\right]^{2}}{k^{3}\left(t^{\prime}-t\right)^{3}}$

in the exponent is not nondimensional, as it should be, and has dimension
$T^{6}L^{-1}$, while the prefactor

$\begin{array}[]{c}\frac{2\sqrt{3}}{\pi
k^{2}\left(t^{\prime}-t\right)^{2}}\end{array}$

has dimension $T^{4}L^{-1}$, instead of the right dimension $TL^{-2}$.

Third, due to yet another typo, the solution given by Eq. (39) does not
converge to the initial condition in the limit $t^{\prime}\rightarrow t$.
Indeed, asymptotically, we have

$\begin{array}[]{c}g\sim
P^{H}\left(\frac{k^{3}\left(t^{\prime}-t\right)^{3}}{6},q^{\prime}-q\right)P^{H}\left(2k\left(t^{\prime}-t\right),\dot{q}^{\prime}-\dot{q}\right)\rightarrow
4\delta\left(q^{\prime}-q\right)\delta\left(\dot{q}^{\prime}-\dot{q}\right).\end{array}$
(48)

Here $P^{H}$ is the standard heat kernel:

$\begin{array}[]{c}P^{H}\left(\tau,\eta\right)=\frac{\exp\left(-\frac{\eta^{2}}{2\tau}\right)}{\sqrt{2\pi\tau}}.\end{array}$
(49)

However, not all is lost. Dimensional analysis shows that the correct solution
$g\left(t,q,\dot{q},t^{\prime},q^{\prime},\dot{q}^{\prime};f,k\right)$ of Eq.
(38) has the form:

$\begin{array}[]{c}g=\frac{\sqrt{3}}{2\pi
k\left(t^{\prime}-t\right)^{2}}\exp\left\\{-\frac{\left[\dot{q}^{\prime}-\dot{q}-f\left(t^{\prime}-t\right)\right]^{2}}{4k\left(t^{\prime}-t\right)}-\frac{3\left[q^{\prime}-q-\frac{\dot{q}^{\prime}+\dot{q}}{2}\left(t^{\prime}-t\right)\right]^{2}}{k\left(t^{\prime}-t\right)^{3}}\right\\},\end{array}$
(50)

which is not far from Kolmogorov’s formula. The correct solution of Eq. (42)
has the form:

$\begin{array}[]{c}g=\frac{\sqrt{3}}{\pi
k\left(t^{\prime}-t\right)^{2}}\exp\left\\{-\frac{\left[\dot{q}^{\prime}-\dot{q}-f\left(t^{\prime}-t\right)\right]^{2}}{2k\left(t^{\prime}-t\right)}-\frac{6\left[q^{\prime}-q-\frac{\dot{q}^{\prime}+\dot{q}}{2}\left(t^{\prime}-t\right)\right]^{2}}{k\left(t^{\prime}-t\right)^{3}}\right\\}.\end{array}$
(51)

###  4.3 Direct Verification

In order to avoid confusion, from now on, we shall change the notation to make
the formulas easier to read. Specifically, we assume that $x$ represents the
position of a particle at time $t$ and $\xi$ its position at time $\tau$,
while $y$ represents its velocity at time $t$, and $\theta$ its velocity at
time $\tau$, so that

$\begin{array}[]{c}\left(t,q,\dot{q}\right)\rightarrow\left(\tau,\xi,\theta\right),\
\ \
\left(t^{\prime},q^{\prime},\dot{q}^{\prime}\right)\rightarrow\left(t,x,y\right).\end{array}$
(52)

One of our objectives is deriving the (corrected) Kolmogorov formula from
first principles using Kelvin Waves. Subsequently, we show how to use it in
the financial mathematics context.

The governing SDE can be written as

$\begin{array}[]{c}dx_{t}=y_{t}dt,\ \ \ x_{\tau}=\xi,\\\ \\\
dy_{t}=bdt+\sqrt{a}dW_{t},\ \ \ y_{\tau}=\theta.\end{array}$ (53)

The corresponding Fokker-Planck- Kolmogorov problem has the form:

$\begin{array}[]{c}P_{t}\left(\tau,\xi,\theta,t,x,y\right)-\frac{1}{2}aP_{yy}\left(\tau,\xi,\theta,t,x,y\right)\\\
\\\
+yP_{x}\left(\tau,\xi,\theta,t,x,y\right)+bP_{y}\left(\tau,\xi,\theta,t,x,y\right)=0,\\\
\\\
P\left(\tau,x,y,\tau,\xi,\theta\right)=\delta\left(x-\xi\right)\delta\left(y-\theta\right).\end{array}$
(54)

We claim that the solution of Eqs (54) is as follows:

$P\left(\tau,\xi,\theta,t,x,y\right)=\frac{\sqrt{3}}{\pi
at^{2}}\exp\left(-\Phi\left(\tau,\xi,\theta,t,x,y\right)\right),$ (55)

where

$\begin{array}[]{c}\Phi\left(\tau,\xi,\theta,t,x,y\right)=\frac{\left(y-\theta-b\left(t-\tau\right)\right)^{2}}{2at}+\frac{6\left(x-\xi-\frac{\left(y+\theta\right)\left(t-\tau\right)}{2}\right)^{2}}{at^{3}}=\frac{A^{2}}{2}+6B^{2},\end{array}$
(56)

and

$\begin{array}[]{c}A=\frac{\left(y-\theta-b\left(t-\tau\right)\right)}{\sqrt{a\left(t-\tau\right)}},\
\ \left[A\right]=1,\\\ \\\
B=\frac{\left(x-\xi-\frac{\left(y+\theta\right)\left(t-\tau\right)}{2}\right)}{\sqrt{a\left(t-\tau\right)^{3}}},\
\ \ \left[B\right]=1.\end{array}$ (57)

Let us check that $P$ satisfies the Fokker-Plank equation and the initial
conditions. We have:

$\begin{array}[]{c}\Phi_{t}=-\left(\frac{A^{2}}{2\left(t-\tau\right)}+\frac{bA}{\sqrt{a\left(t-\tau\right)}}+\frac{18B^{2}}{\left(t-\tau\right)}+\frac{6\left(y+\theta\right)B}{\sqrt{a\left(t-\tau\right)^{3}}}\right),\\\
\\\ \ \Phi_{x}=\frac{12B}{\sqrt{a\left(t-\tau\right)^{3}}},\ \ \
\Phi_{y}=\frac{A-6B}{\sqrt{a\left(t-\tau\right)}},\ \ \
\Phi_{yy}=\frac{4}{a\left(t-\tau\right)},\end{array}$ (58)
$\begin{array}[]{c}\frac{P_{t}}{P}=-\frac{2}{\left(t-\tau\right)}-\Phi_{t},\ \
\ \frac{P_{x}}{P}=-\Phi_{x},\ \ \ \frac{P_{y}}{P}=-\Phi_{y},\ \ \
\frac{P_{yy}}{P}=-\Phi_{yy}+\Phi_{x}^{2},\end{array}$ (59)

so that

$\begin{array}[]{c}P_{t}^{K}-\frac{1}{2}aP_{yy}^{K}+yP_{x}^{K}+bP_{y}^{K}\\\
\\\
=P^{K}\left(-\frac{2}{\left(t-\tau\right)}-\Phi_{t}+\frac{1}{2}a\Phi_{yy}-\frac{1}{2}a\Phi_{y}^{2}-y\Phi_{x}-b\Phi_{y}\right)\\\
\\\
=P^{K}\left(-\frac{2}{\left(t-\tau\right)}+\frac{A^{2}}{2\left(t-\tau\right)}+\frac{bA}{\sqrt{a\left(t-\tau\right)}}+\frac{18B^{2}}{\left(t-\tau\right)}+\frac{6\left(y+\theta\right)B}{\sqrt{a\left(t-\tau\right)^{3}}}\right.\\\
\\\
\left.+\frac{2}{\left(t-\tau\right)}-\frac{\left(A-6B\right)^{2}}{2\left(t-\tau\right)}-\frac{12yB}{\sqrt{a\left(t-\tau\right)^{3}}}-\frac{b\left(A-6B\right)}{\sqrt{a\left(t-\tau\right)}}\right)\\\
\\\ =0.\end{array}$ (60)

When $\left(t-\tau\right)\rightarrow 0$ we have the following asymptotic
expression:

$\begin{array}[]{c}P^{K}\left(\tau,\xi,\theta,t,x,y\right)\sim
P^{H}\left(\frac{a\left(t-\tau\right)^{3}}{12},x-\xi\right)P^{H}\left(a\left(t-\tau\right),y-\theta\right)\\\
\\\
\rightarrow\delta\left(x-\xi\right)\delta\left(y-\theta\right).\end{array}$
(61)

###  4.4 Solution via Kelvin Waves

Now, we wish to derive Kolmogorov’s formula by using Kelvin waves (or an
affine ansatz). To this end, we start with the problem of the form:

$\begin{array}[]{c}H_{t}\left(\tau,\xi,\theta,t,x,y,k,l\right)-\frac{1}{2}aH_{yy}\left(\tau,\xi,\theta,t,x,y,k,l\right)\\\
\\\
+yH_{x}\left(\tau,\xi,\theta,t,x,y,k,l\right)+bH_{y}\left(\tau,\xi,\theta,t,x,y,k,l\right)=0,\\\
\\\
H\left(\tau,x,y,\tau,\xi,\theta,k,l\right)=\exp\left(ik\left(x-\xi\right)+il\left(y-\theta\right)\right).\end{array}$
(62)

Here

$\begin{array}[]{c}\ \left[k\right]=L^{-1},\ \ \ \left[l\right]=TL^{-1},\ \
\left[H\right]=1.\end{array}$ (63)

By using the well-known results concerning the inverse Fourier transform of
the $\delta$-function, we get:

$\begin{array}[]{c}P\left(\tau,\xi,\theta,t,x,y\right)=\frac{1}{\left(2\pi\right)^{2}}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}H\left(\tau,\xi,\theta,t,x,y,k,l\right)dkdl.\end{array}$
(64)

To calculate $H$, we use the affine ansatz and represent it in the form:

$\begin{array}[]{c}H\left(\tau,\xi,\theta,t,x,y,k,l\right)=\exp\left(\Psi\left(\tau,\xi,\theta,t,x,y,k,l\right)\right),\end{array}$
(65)

where

$\begin{array}[]{c}\Psi\left(\tau,\xi,\theta,t,x,y,k,l\right)=\alpha\left(t\right)+ik\left(x-\xi\right)+i\beta\left(t\right)y-il\theta.\end{array}$
(66)

and

$\begin{array}[]{c}\frac{H_{t}}{H}=\Psi_{t}=\left(\alpha^{\prime}\left(t\right)+i\beta^{\prime}\left(t\right)y\right),\
\ \ \frac{H_{x}}{H}=\Psi_{x}=ik,\\\ \\\
\frac{H_{y}}{H}=\Psi_{y}=i\beta\left(t\right),\ \ \
\frac{H_{yy}}{H}=\Psi_{y}^{2}=-\beta^{2}\left(t\right).\end{array}$ (67)

Using this ansatz in Eqs (62) yields:

$\begin{array}[]{c}\alpha^{\prime}\left(t\right)+\frac{1}{2}a\beta^{2}\left(t\right)+i\beta^{\prime}\left(t\right)y+iky+ib\beta\left(t\right)=0,\\\
\\\ \alpha\left(\tau\right)=0,\ \ \ \beta\left(\tau\right)=l,\end{array}$ (68)

so that

$\begin{array}[]{c}\alpha^{\prime}\left(t\right)+\frac{1}{2}a\beta^{2}\left(t\right)+ib\beta\left(t\right)=0,\
\ \ \alpha\left(\tau\right)=0,\\\ \\\ \beta^{\prime}\left(t\right)+k=0,\ \ \
\beta\left(\tau\right)=l.\end{array}$ (69)

Straightforward calculation shows:

$\begin{array}[]{c}\beta\left(t\right)=-k\left(t-\tau\right)+l,\\\ \\\
\alpha\left(t\right)=-\frac{1}{2}a\left(\frac{k^{2}\left(t-\tau\right)^{3}}{3}-kl\left(t-\tau\right)^{2}+l^{2}\left(t-\tau\right)\right)-ib\left(-\frac{k\left(t-\tau\right)^{2}}{2}+l\left(t-\tau\right)\right).\end{array}$
(70)

Equations (64, 65, 66, 70) yield:

$\begin{array}[]{c}P\left(\tau,\xi,\theta,t,x,y\right)=\frac{1}{\left(2\pi\right)^{2}}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\\\
\\\
\times\exp\left(-\frac{1}{2}a\left(\frac{k^{2}\left(t-\tau\right)^{3}}{3}-kl\left(t-\tau\right)^{2}+l^{2}\left(t-\tau\right)\right)\right.\\\
\\\
\left.+ik\left(x-\xi-y\left(t-\tau\right)+\frac{b\left(t-\tau\right)^{2}}{2}\right)+il\left(y-\theta-b\left(t-\tau\right)\right)\right)dkdl.\end{array}$
(71)

It is clear that $P\left(\tau,\xi,\theta,t,x,y\right)$ can be viewed as the
characteristic function of the Gaussian density in the $\left(k,l\right)$
space, evaluated at the point $\left(x-\xi-y\left(t-\tau\right)\right.$
$\left.+\frac{b\left(t-\tau\right)^{2}}{2},y-\theta-b\left(t-\tau\right)\right)$:

$\begin{array}[]{c}P\left(\tau,\xi,\theta,t,x,y\right)=\frac{\left(\det\left(\mathfrak{C}\right)\right)^{1/2}}{2\pi}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\\\
\\\ \times
G\left(t,k,l\right)\exp\left(ik\left(x-\xi-y\left(t-\tau\right)+\frac{b\left(t-\tau\right)^{2}}{2}\right)+il\left(y-\theta-b\left(t-\tau\right)\right)\right)dkdl,\end{array}$
(72)

where

$\begin{array}[]{c}G\left(t,k,l\right)=\frac{1}{2\pi\left(\det\left(\mathfrak{C}\right)\right)^{1/2}}\exp\left(-\frac{1}{2}\left(\begin{array}[]{c}k\\\
l\end{array}\right)\cdot\mathfrak{C}^{-1}\left(\begin{array}[]{c}k\\\
l\end{array}\right)\right),\end{array}$ (73)

and

$\begin{array}[]{c}\mathfrak{C}\mathbb{=}\left(\begin{array}[]{cc}\frac{12}{a\left(t-\tau\right)^{3}}&\frac{6}{a\left(t-\tau\right)^{2}}\\\
\frac{6}{a\left(t-\tau\right)^{2}}&\frac{4}{a\left(t-\tau\right)}\end{array}\right),\\\
\\\
\det\left(\mathfrak{C}\right)=\frac{12}{a^{2}\left(t-\tau\right)^{4}}.\end{array}$
(74)

As before, $\cdot$ denotes the scalar product. Accordingly,

$\begin{array}[]{c}P\left(\tau,\xi,\theta,t,x,y\right)=\frac{\sqrt{3}}{\pi
a\left(t-\tau\right)^{2}}\exp\left(-\Omega\left(\tau,\xi,\theta,t,x,y\right)\right),\end{array}$
(75)

where

$\begin{array}[]{c}\Omega\left(\tau,\xi,\theta,t,x,y\right)=\\\ \\\
=\frac{1}{2}\left(\begin{array}[]{c}x-\xi-y\left(t-\tau\right)+\frac{b\left(t-\tau\right)^{2}}{2}\\\
y-\theta-b\left(t-\tau\right)\end{array}\right)\cdot\mathfrak{C}\left(\begin{array}[]{c}x-\xi-y\left(t-\tau\right)+\frac{b\left(t-\tau\right)^{2}}{2}\\\
y-\theta-b\left(t-\tau\right)\end{array}\right)\\\ \\\
=\frac{6\left(x-\xi-y\left(t-\tau\right)+\frac{b\left(t-\tau\right)^{2}}{2}\right)^{2}}{a\left(t-\tau\right)^{3}}+\frac{6\left(x-\xi-y\left(t-\tau\right)+\frac{b\left(t-\tau\right)^{2}}{2}\right)\left(y-\theta-b\left(t-\tau\right)\right)}{a\left(t-\tau\right)^{2}}+\frac{2\left(y-\theta-b\left(t-\tau\right)\right)^{2}}{a\left(t-\tau\right)}\\\
\\\ =6\left(B-\frac{A}{2}\right)^{2}+6\left(B-\frac{A}{2}\right)A+2A^{2}\\\
\\\ =6B^{2}-6AB+\frac{3A^{2}}{2}+6AB-3A^{2}+2A^{2}\\\ \\\
=\frac{A^{2}}{2}+6B^{2}\\\ \\\
=\Phi\left(t,x,y\,,\xi,\theta\right),\end{array}$ (76)

as expected. This calculation completes our derivation of the corrected
Kolmogorov formula.

$P$ is a bivariate Gaussian distribution. Completing the square, we get:

$\begin{array}[]{c}\Phi=\frac{\left(y-\theta-b\left(t-\tau\right)\right)^{2}}{2a\left(t-\tau\right)}+\frac{6\left(x-\xi-\frac{\left(y+\theta\right)\left(t-\tau\right)}{2}\right)^{2}}{a\left(t-\tau\right)^{3}}\\\
\\\
=\frac{6}{a\left(t-\tau\right)^{3}}\left(x-p\right)^{2}-\frac{6}{a\left(t-\tau\right)^{2}}\left(x-p\right)\left(y-q\right)+\frac{2}{a\left(t-\tau\right)}\left(y-q\right)^{2},\end{array}$
(77)

and represent $P$ the form:

$\begin{array}[]{c}P\left(\tau,\xi,\theta,t,x,y\right)=\frac{\exp\left(-\frac{1}{2\left(1-\rho^{2}\right)}\left(\frac{\left(x-p\right)^{2}}{\sigma_{x}^{2}}-\frac{2\rho\left(x-p\right)\left(y-q\right)}{\sigma_{x}\sigma_{y}}+\frac{\left(y-q\right)^{2}}{\sigma_{y}^{2}}\right)\right)}{2\pi\sigma_{x}\sigma_{y}\sqrt{1-\rho^{2}}},\end{array}$
(78)

where

$\begin{array}[]{c}\sigma_{x}=\sqrt{\frac{a\left(t-\tau\right)^{3}}{3}},\ \ \
\sigma_{y}=\sqrt{a\left(t-\tau\right)},\ \ \ \rho=\frac{\sqrt{3}}{2},\\\ \\\
p=\xi+\theta\left(t-\tau\right)+\frac{b\left(t-\tau\right)^{2}}{2},\ \
q=\theta+b\left(t-\tau\right).\end{array}$ (79)

##  5 Gaussian Processes

###  5.1 The General Case

We now consider a more general Kolmogorov-type SDE solvable via the Kelvin (or
affine) ansatz:

$\begin{array}[]{c}d\mathbf{x}_{t}=\left(\mathbf{b}^{\left(x\right)}+\mathfrak{B}^{\left(xx\right)}\mathbf{x}_{t}+\mathfrak{B}^{\left(xy\right)}\mathbf{y}_{t}\right)dt,\
\ \ \mathbf{x}_{\tau}=\mathbf{\xi,}\\\ \\\
d\mathbf{y}_{t}=\left(\mathbf{b}^{\left(y\right)}+\mathfrak{B}^{\left(yx\right)}\mathbf{x}_{t}+\mathfrak{B}^{\left(yy\right)}\mathbf{y}_{t}\right)dt+\mathfrak{\mathbf{\Sigma}}^{\left(yy\right)}d\mathbf{W}_{t}^{\left(y\right)},\
\ \ \mathbf{y}_{\tau}=\mathbf{\theta,}\end{array}$ (80)

where $\mathbf{x}_{{}_{t}}=\left(x_{m}\right)$ and
$\mathbf{b}^{\left(x\right)}=\left(b_{m}^{\left(x\right)}\right)$ are
$\left(M\times 1\right)$ column vectors,
$\mathbf{y}_{t}=\left(y_{n}\right)^{T}$ and
$\mathbf{b}^{\left(y\right)}=\left(b_{n}^{\left(y\right)}\right)$ are
$\left(N\times 1\right)$ column vectors,
$\mathfrak{B}^{\left(xx\right)}=\left(b_{mm^{\prime}}^{\left(xx\right)}\right)$
is $\left(M\times M\right)$ matrix,
$\mathfrak{B}^{\left(xy\right)}=\left(b_{mn}^{\left(xy\right)}\right)$ is
$\left(M\times N\right)$ matrix,
$\mathfrak{B}^{\left(yx\right)}=\left(b_{nm}^{\left(xx\right)}\right)$ is
$\left(N\times M\right)$ matrix,
$\mathfrak{B}^{\left(yy\right)}=\left(b_{nn^{\prime}}^{\left(xx\right)}\right)$
and $\mathbf{\Sigma}^{\left(yy\right)}=\left(\sigma_{nn^{\prime}}\right)$ are
$\left(N\times N\right)$ matrices. Below, we assume that the corresponding
coefficients are time-dependent. As usual, $\mathbf{W}_{t}^{\left(y\right)}$
is a standard $N$-dimensional Brownian motion.

More compactly, we can write the system of SDEs as follows:

$\begin{array}[]{cc}d\mathbf{z}_{t}=\left(\mathbf{b}^{\left(z\right)}+\mathfrak{B}^{\left(zz\right)}\mathbf{z}_{t}\right)dt+\left(\begin{array}[]{c}0\\\
\mathbf{\Sigma}^{\left(yy\right)}d\mathbf{W}_{t}^{\left(y\right)}\end{array}\right),&\mathbf{z}_{\tau}=\left(\begin{array}[]{c}\mathbf{\xi}\\\
\mathbf{\theta}\end{array}\right),\end{array}$ (81)

where

$\begin{array}[]{c}\mathbf{z}=\left(\begin{array}[]{c}\mathbf{x}\\\
\mathbf{y}\end{array}\right),\ \ \
\mathbf{b}^{\left(z\right)}=\left(\begin{array}[]{c}\mathbf{b}^{\left(x\right)}\\\
\mathbf{b}^{\left(y\right)}\end{array}\right),\ \ \
\mathfrak{B}^{\left(zz\right)}=\left(\begin{array}[]{cc}\mathfrak{B}^{\left(xx\right)}&\mathfrak{B}^{\left(xy\right)}\\\
\mathfrak{B}^{\left(yx\right)}&\mathfrak{B}^{\left(yy\right)}\end{array}\right),\end{array}$
(82)

so that $\mathbf{z}$ and $\mathbf{b}^{\left(z\right)}$ are $\left(I\times
1\right)$ column vector, and $\mathfrak{B}^{\left(zz\right)}$ is a
$\left(I\times I\right)$ matrix, where $I=M+N$. We also introduce a scalar
$b^{\left(z\right)}=\mathrm{Tr}\left(\mathfrak{B}^{\left(zz\right)}\right)=\mathrm{Tr}\left(\mathfrak{B}^{\left(xx\right)}\right)+\mathrm{Tr}\left(\mathfrak{B}^{\left(yx\right)}\right)$.

The corresponding Fokker-Plank problem has the form:

$\begin{array}[]{c}P_{t}\left(\tau,\mathbf{\zeta,}t,\mathbf{z}\right)-\frac{1}{2}\mathop{\displaystyle\sum}\mathop{\displaystyle\sum}\mathfrak{A}P_{\mathbf{yy}}\left(\tau,\mathbf{\zeta,}t,\mathbf{z}\right)\\\
\\\
+\left(\mathbf{b}^{\left(z\right)}+\mathfrak{B}^{\left(z\right)}\mathbf{z}\right)\cdot
P_{\mathbf{z}}\left(\tau,\mathbf{\zeta,}t,\mathbf{z}\right)+b^{\left(z\right)}P\left(\tau,\mathbf{\zeta,}t,\mathbf{z}\right)=0,\\\
\\\
P\left(\tau,\mathbf{z},\tau,\mathbf{\zeta}\right)=\delta\left(\mathbf{x}-\mathbf{\xi}\right)\delta\left(\mathbf{y}-\mathbf{\theta}\right),\end{array}$
(83)

where $\mathfrak{A}$ is the covariance matrix,

$\begin{array}[]{c}\mathfrak{A}=\left(a_{nn^{\prime}}\right)=\sum\limits_{k=1}^{N}\sigma_{nk}\sigma_{n^{\prime}k}=\mathbf{\Sigma}^{\left(yy\right)}\mathbf{\Sigma}^{\left(yy\right)T}.\end{array}$
(84)

Explicitly,

$\begin{array}[]{c}P_{t}-\frac{1}{2}\sum\limits_{n=1}^{N}\sum\limits_{n^{\prime}=1}^{N}a_{nn^{\prime}}P_{y_{n}y_{n^{\prime}}}+\sum\limits_{i=1}^{I}\left(b_{i}^{\left(z\right)}+\sum\limits_{i^{\prime}=1}^{I}b_{ii^{\prime}}^{\left(zz\right)}x_{i^{\prime}}\right)P_{z_{i}}+b^{\left(z\right)}P=0,\end{array}$
(85)

We use the same Kelvin ansatz as before and represent $P$ in the form:

$\begin{array}[]{c}P\left(\tau,\mathbf{\zeta,}t,\mathbf{z}\right)=\frac{1}{\left(2\pi\right)^{I}}\int_{-\infty}^{\infty}...\int_{-\infty}^{\infty}H\left(\tau,\mathbf{\zeta,}t,\mathbf{z,m}\right)d\mathbf{m},\\\
\\\
H\left(\tau,\mathbf{\zeta,}t,\mathbf{z,m}\right)=\exp\left(\Psi\left(\tau,\mathbf{\zeta,}t,\mathbf{z,m}\right)\right),\\\
\\\
\Psi\left(\tau,\mathbf{\zeta,}t,\mathbf{z,m}\right)=\alpha\left(t\right)+i\Upsilon\left(t\right)\cdot\mathbf{z}-i\mathbf{m}\cdot\mathbf{\zeta},\end{array}$
(86)

where $\mathbf{m}\mathcal{=}\left(m_{i}\right)$ is an $\left(I\times 1\right)$
column vector, $\mathbf{m}=\left(\mathbf{k,l}\right)^{T}$,
$\mathbf{k}\mathcal{=}\left(k_{m}\right)$ is an $\left(M\times 1\right)$
column vector, $\mathbf{l}=\left(l_{n}\right)~{}$is an $\left(N\times
1\right)$ column vector, $\mathbf{\Upsilon=}\left(\upsilon_{i}\right)$ is an
$\left(I\times 1\right)$ column vector, $\
\mathbf{\Upsilon}=\left(\mathbf{\Gamma},\mathbf{\Delta}\right)^{T}$,
$\mathbf{\Gamma}=\left(\gamma_{m}^{\left(x\right)}\right)~{}$is an
$\left(M\times 1\right)$ column vector,
$\mathbf{\Delta}=\left(\delta_{n}^{\left(y\right)}\right)~{}$is an
$\left(N\times 1\right)$ column vector, and

$\begin{array}[]{c}\alpha\left(\tau\right)=0,\ \ \
\mathbf{\Upsilon}\left(\tau\right)=\left(\mathbf{\Gamma}\left(\tau\right),\mathbf{\Delta}\left(\tau\right)\right)^{T}=\mathbf{m}=\left(\mathbf{k,l}\right)^{T}.\end{array}$
(87)

As before:

$\begin{array}[]{c}\frac{H_{t}}{H}=\Psi_{t}=\left(\alpha^{\prime}\left(t\right)+i\mathbf{\Upsilon}^{\prime}\left(t\right)\cdot\mathbf{z}\right),\
\ \
\frac{H_{\mathbf{x}}}{H}=\Psi_{\mathbf{x}}=i\mathbf{\Gamma}\left(t\right),\\\
\\\
\frac{H_{\mathbf{y}}}{H}=\Psi_{\mathbf{y}}=i\mathbf{\Delta}\left(t\right),\ \
\
\frac{H_{\mathbf{yy}}}{H}=\Psi_{\mathbf{y}}^{2}=-\mathbf{\Delta}\left(t\right)\mathbf{\Delta}^{T}\left(t\right).\end{array}$
(88)

The equations for $\alpha,\mathbf{\Upsilon}$ have the form:

$\begin{array}[]{c}\alpha^{\prime}\left(t\right)+i\mathbf{\Upsilon}^{\prime}\left(t\right)\cdot\mathbf{z}+\frac{1}{2}\mathbf{\Delta}\left(t\right)\cdot\mathfrak{A}\mathbf{\Delta}\left(t\right)+i\mathbf{\Upsilon}\left(t\right)\cdot\left(\mathbf{b}^{\left(z\right)}+\mathfrak{B}^{\left(zz\right)}\mathbf{z}\right)+b^{\left(z\right)}=0.\end{array}$
(89)

Accordingly,

$\begin{array}[]{c}\mathbf{\Upsilon}^{\prime}\left(t\right)+\mathfrak{B}^{\left(zz\right)T}\mathbf{\Upsilon}\left(t\right)=0,\
\ \
\mathbf{\Upsilon}\left(\tau\right)=\mathbf{m}=\left(\mathbf{k,l}\right)^{T},\end{array}$
(90)

and

$\begin{array}[]{c}\alpha^{\prime}\left(t\right)+\frac{1}{2}\mathbf{\Delta}\left(t\right)\cdot\mathfrak{A}\mathbf{\Delta}\left(t\right)+i\mathbf{\Upsilon}\left(t\right)\cdot\mathbf{b}^{\left(z\right)}+b^{\left(z\right)}=0,\
\ \ \alpha\left(\tau\right)=0.\end{array}$ (91)

Let $\mathfrak{L}\left(\tau,t\right)$ is the fundamental solution of the
homogeneous ODE (91), i.e., the matrix such that

$\begin{array}[]{c}\mathbf{\mathfrak{L}}^{\prime}\left(\tau,t\right)+\mathfrak{B}^{\left(zz\right)T}\mathbf{\mathfrak{L}}\left(\tau,t\right)=0,\
\ \ \mathbf{\mathfrak{L}}\left(\tau,\tau\right)=\mathfrak{I},\end{array}$ (92)

where $\mathfrak{I}$ is the identity matrix. The well-known Liouville’s
formula, which we shall use shortly, yields

$\begin{array}[]{c}\det\left(\mathfrak{L}\left(\tau,t\right)\right)=\exp\left(-\int_{\tau}^{t}b^{\left(z\right)}\left(s\right)ds\right).\end{array}$
(93)

The solution of Eq. (91) has the form:

$\begin{array}[]{c}\mathbf{\Upsilon}\left(t\right)=\mathbf{\mathfrak{L}}\left(\tau,t\right)\mathbf{m}.\end{array}$
(94)

It is convenient to write $\mathfrak{L}\left(\tau,t\right)$ in the block form:

$\begin{array}[]{c}\mathfrak{L}\left(\tau,t\right)=\left(\begin{array}[]{cc}\mathfrak{L}^{\left(xx\right)}\left(\tau,t\right)&\mathfrak{L}^{\left(xy\right)}\left(\tau,t\right)\\\
\mathfrak{L}^{\left(yx\right)}\left(\tau,t\right)&\mathfrak{L}^{\left(yy\right)}\left(\tau,t\right)\end{array}\right).\end{array}$
(95)

It follows from Eq. (90) that

$\begin{array}[]{c}\alpha\left(t\right)=-\frac{1}{2}\mathbf{m}\cdot\mathfrak{C}^{-1}\left(\tau,t\right)\mathbf{m}-i\mathbf{m\cdot
d}^{\left(z\right)}\left(\tau,t\right)-\varpi\left(\tau,t\right),\end{array}$
(96)

where $\mathfrak{C}^{-1}$ is an $I\times I$ positive-definite matrix split
into four blocks of the form:

$\begin{array}[]{c}\mathfrak{C}^{-1}\left(\tau,t\right)\\\ \\\
=\left(\begin{array}[]{cc}\int_{\tau}^{t}\mathfrak{L}^{\left(yx\right)T}\left(\tau,s\right)\mathfrak{A}\left(s\right)\mathfrak{L}^{\left(yx\right)}\left(\tau,s\right)ds&\int_{\tau}^{t}\mathfrak{L}^{\left(yx\right)T}\left(\tau,s\right)\mathfrak{A}\left(s\right)\mathfrak{L}^{\left(yy\right)}\left(\tau,s\right)ds\\\
\int_{\tau}^{t}\mathfrak{L}^{\left(yy\right)T}\left(\tau,s\right)\mathfrak{A}\left(s\right)\mathfrak{L}^{\left(yx\right)}\left(\tau,s\right)ds&\int_{\tau}^{t}\mathfrak{L}^{\left(yy\right)T}\left(\tau,s\right)\mathfrak{A}\left(s\right)\mathfrak{L}^{\left(yy\right)}\left(\tau,s\right)ds\end{array}\right),\end{array}$
(97)

while
$\mathbf{d}^{\left(z\right)}=\left(\mathbf{d}^{\left(x\right)},\mathbf{d}^{\left(y\right)}\right)^{T}$,
$\mathbf{d}^{\left(x\right)}$ and $\mathbf{d}^{\left(y\right)}$ are
$\left(M\times 1\right)$ and $\left(N\times 1\right)$ column vectors, and
$\varpi$ is a scalar:

$\begin{array}[]{c}\mathbf{d}^{\left(z\right)}\left(\tau,t\right)=\int_{\tau}^{t}\mathfrak{L}^{T}\left(\tau,s\right)\mathbf{b}^{\left(z\right)}\left(s\right)ds,\end{array}$
(98)
$\begin{array}[]{c}\varpi\left(\tau,t\right)=\int_{\tau}^{t}b^{\left(z\right)}\left(s\right)ds.\end{array}$
(99)

Accordingly,

$\begin{array}[]{c}\Psi\left(\tau,\mathbf{\zeta,}t,\mathbf{z},\mathbf{m}\right)\\\
\\\
=-\frac{1}{2}\mathbf{m}\cdot\mathfrak{C}^{-1}\left(\tau,t\right)\mathbf{m}+i\mathbf{\mathfrak{L}}\left(\tau,t\right)\mathbf{m}\cdot\mathbf{z}-i\mathbf{m}\cdot\left(\mathbf{d}^{\left(z\right)}\left(\tau,t\right)+\mathbf{\zeta}\right)-\varpi\left(\tau,t\right)\\\
\\\
=-\frac{1}{2}\mathbf{m}\cdot\mathfrak{C}^{-1}\left(\tau,t\right)\mathbf{m}+i\mathbf{m}\cdot\left(\mathfrak{L}^{T}\left(\tau,t\right)\mathbf{z-d}^{\left(z\right)}\left(\tau,t\right)\mathbf{-\zeta}\right)-\varpi\left(\tau,t\right).\end{array}$
(100)

Thus,

$\begin{array}[]{c}P\left(\tau,\mathbf{\zeta,}t,\mathbf{z}\right)=\frac{\det\left(\mathfrak{C}\right)^{1/2}\exp\left(-\varpi\left(\tau,t\right)\right)}{\left(2\pi\right)^{I/2}}\int_{-\infty}^{\infty}...\int_{-\infty}^{\infty}G\left(t,\mathbf{m}\right)\\\
\\\
\times\exp\left(i\mathbf{m}\cdot\left(\mathfrak{L}^{T}\left(\tau,t\right)\mathbf{z-d}^{\left(z\right)}\left(\tau,t\right)\mathbf{-\zeta}\right)\right)d\mathbf{m},\end{array}$
(101)

where $G\left(t,\mathbf{m}\right)$ is the density of a multivariate Gaussian
distribution in the $\mathbf{m}$-space. It is clear that
$P\left(\tau,\mathbf{\zeta,}t,\mathbf{z}\right)$ is proportional to the
characteristic function of $G$ evaluated at the point
$\left(\mathfrak{L}^{T}\left(\tau,t\right)\mathbf{z-d}^{\left(z\right)}\left(\tau,t\right)\mathbf{-\zeta}\right)$,
so that

$\begin{array}[]{c}P\left(\tau,\mathbf{\zeta,}t,\mathbf{z}\right)=\frac{\det\left(\mathfrak{C}\right)^{1/2}\exp\left(-\varpi\left(\tau,t\right)\right)}{\left(2\pi\right)^{I/2}}\\\
\\\
\times\exp\left(-\frac{1}{2}\left(\mathfrak{L}^{T}\left(\tau,t\right)\mathbf{z-d}^{\left(z\right)}\left(\tau,t\right)\mathbf{-\zeta}\right)\cdot\mathfrak{C}\left(\mathfrak{L}^{T}\left(\tau,t\right)\mathbf{z-d}^{\left(z\right)}\left(\tau,t\right)\mathbf{-\zeta}\right)\right).\end{array}$
(102)

Simple calculation using Eq. (93) allows one to rewrite (102) it the standard
Gaussian form:

$\begin{array}[]{c}P\left(\tau,\mathbf{\zeta,}t,\mathbf{z}\right)=\mathrm{N}\left(\mathbf{r}\left(\tau,t\right),\mathfrak{H}\left(\tau,t\right)\right),\end{array}$
(103)

where

$\begin{array}[]{c}\mathfrak{H}\left(\tau,t\right)=\left(\mathfrak{L}^{T}\left(\tau,t\right)\right)^{-1}\mathfrak{C}^{-1}\left(\tau,t\right)\mathfrak{L}^{-1}\left(\tau,t\right),\\\
\\\
\mathbf{r}\left(\tau,t\right)\mathbf{=}\left(\mathfrak{L}^{T}\left(\tau,t\right)\right)^{-1}\left(\mathbf{d}^{\left(z\right)}\left(\tau,t\right)+\mathbf{\zeta}\right).\end{array}$
(104)

###  5.2 Example: Kolmogorov Process

We now wish to extend the Kolmogorov formula to the case when $b$ and $a$ are
functions of time. In order to make connections with other cases a bit easier
to draw, we use the specific notation: $b\left(t\right)$ and
$a\left(t\right)$. The corresponding SDE has the form:

$\begin{array}[]{c}dx_{t}=y_{t}dt,\ \ \ x_{\tau}=\xi,\\\ \\\
dy_{t}=b\left(t\right)dt+\sqrt{a\left(t\right)}dW_{t},\ \ \
y_{\tau}=\theta.\end{array}$ (105)

Accordingly, Eq. (92) can be written as follows:

$\begin{array}[]{c}\mathfrak{L}^{\prime}\left(\tau,t\right)+\left(\begin{array}[]{cc}0&0\\\
1&0\end{array}\right)\mathfrak{L}\left(\tau,t\right)=0,\ \
\mathfrak{L}\left(\tau,\tau\right)=\left(\begin{array}[]{cc}1&0\\\
0&1\end{array}\right),\end{array}$ (106)

so that

$\begin{array}[]{c}\mathfrak{L}\left(\tau,t\right)=\left(\begin{array}[]{cc}1&0\\\
-\left(t-\tau\right)&1\end{array}\right),\\\ \\\
\mathfrak{L}^{-1}\left(\tau,t\right)=\left(\begin{array}[]{cc}1&0\\\
\left(t-\tau\right)&1\end{array}\right).\end{array}$ (107)

Once we know $\mathfrak{L}\left(\tau,t\right)$, we can compute
$\mathfrak{C}^{-1}\left(\tau,t\right)$,
$\mathbf{d}^{\left(z\right)}\left(\tau,t\right)=\left(d^{\left(x\right)}\left(\tau,t\right),d^{\left(y\right)}\left(\tau,t\right)\right)^{T}$,
$\varpi\left(\tau,t\right)$:

$\begin{array}[]{c}\mathfrak{C}^{-1}\left(\tau,t\right)=\left(\begin{array}[]{cc}\psi_{2}\left(\tau,t\right)&-\psi_{1}\left(\tau,t\right)\\\
-\psi_{1}\left(\tau,t\right)&\psi_{0}\left(\tau,t\right)\end{array}\right),\\\
\\\
\mathbf{d}^{\left(z\right)}\left(\tau,t\right)=\left(\begin{array}[]{c}d^{\left(x\right)}\left(\tau,t\right)\\\
d^{\left(y\right)}\left(\tau,t\right)\end{array}\right)=\left(\begin{array}[]{c}-\phi_{1}\left(\tau,t\right)\\\
\phi_{0}\left(\tau,t\right)\end{array}\right),\\\ \\\
\varpi\left(\tau,t\right)=0,\end{array}$ (108)

where

$\begin{array}[]{c}\phi_{i}\left(\tau,t\right)=\int_{\tau}^{t}\left(s-\tau\right)^{i}b\left(s\right)ds,\
\ \ \ \
\psi_{i}\left(\tau,t\right)=\int_{\tau}^{t}\left(s-\tau\right)^{i}a\left(s\right)ds.\end{array}$
(109)

Next, we calculate $\mathfrak{H}\left(\tau,t\right)$, and
$\mathbf{r}\left(\tau,t\right)$:

$\begin{array}[]{c}\mathfrak{H}\left(\tau,t\right)=\left(\mathfrak{L}^{T}\left(\tau,t\right)\right)^{-1}\mathfrak{C}^{-1}\left(\tau,t\right)\mathfrak{L}^{-1}\left(\tau,t\right)\\\
\\\ =\left(\begin{array}[]{cc}1&\left(t-\tau\right)\\\
0&1\end{array}\right)\left(\begin{array}[]{cc}\psi_{2}\left(\tau,t\right)&-\psi_{1}\left(\tau,t\right)\\\
-\psi_{1}\left(\tau,t\right)&\psi_{0}\left(\tau,t\right)\end{array}\right)\left(\begin{array}[]{cc}1&0\\\
\left(t-\tau\right)&1\end{array}\right)\\\ \\\
=\left(\begin{array}[]{cc}\psi_{0}\left(\tau,t\right)\left(t-\tau\right)^{2}-2\psi_{1}\left(\tau,t\right)\left(t-\tau\right)+\psi_{2}\left(\tau,t\right)&\psi_{0}\left(\tau,t\right)\left(t-\tau\right)-\psi_{1}\left(\tau,t\right)\\\
\psi_{0}\left(\tau,t\right)\left(t-\tau\right)-\psi_{1}\left(\tau,t\right)&\psi_{0}\left(\tau,t\right)\end{array}\right),\end{array}$
(110)
$\begin{array}[]{c}\mathbf{r}\left(\tau,t\right)\mathbf{=}\left(\mathfrak{L}^{T}\left(\tau,t\right)\right)^{-1}\left(\mathbf{d}^{\left(z\right)}\left(\tau,t\right)+\mathbf{\zeta}\right)\\\
\\\
=\left(\begin{array}[]{c}-\phi_{1}\left(\tau,t\right)+\xi+\left(\phi_{0}\left(\tau,t\right)+\theta\right)\left(t-\tau\right)\\\
\phi_{0}\left(\tau,t\right)+\theta\end{array}\right).\end{array}$ (111)

Accordingly, $P\left(\tau,\xi,\theta,t,x,y\right)$ is a bivariate Gaussian
distribution of the form (103), with

$\begin{array}[]{c}\sigma_{x}^{2}\left(\tau,t\right)=\psi_{0}\left(\tau,t\right)\left(t-\tau\right)^{2}-2\psi_{1}\left(\tau,t\right)\left(t-\tau\right)+\psi_{2}\left(\tau,t\right),\
\ \ \sigma_{y}^{2}=\psi_{0}\left(\tau,t\right),\\\ \\\
\rho\left(\tau,t\right)=\frac{\left(\psi_{0}\left(\tau,t\right)\left(t-\tau\right)-\psi_{1}\left(\tau,t\right)\right)}{\sqrt{\psi_{0}\left(\tau,t\right)\left(\psi_{0}\left(\tau,t\right)\left(t-\tau\right)^{2}-2\psi_{1}\left(\tau,t\right)\left(t-\tau\right)+\psi_{2}\left(\tau,t\right)\right)}},\\\
\\\
p\left(\tau,t\right)=-\phi_{1}\left(\tau,t\right)+\xi+\left(\phi_{0}\left(\tau,t\right)+\theta\right)\left(t-\tau\right),\
\ \ q\left(\tau,t\right)=\phi_{0}\left(\tau,t\right)+\theta.\end{array}$ (112)

It is easy to see that Eq. (112) coincides with Eq. (78) when $a,b$ are
constant.

###  5.3 Example: OU Process

####  5.3.1 OU Process

For benchmarking purposes, we start with deriving the with the well-known
t.p.d.f. for the OU process:

$\begin{array}[]{c}dy=\left(\chi\left(t\right)-\kappa\left(t\right)y\right)dt+\varepsilon\left(t\right)dZ,\
\ \ y\left(\tau\right)=\theta.\end{array}$ (113)

The corresponding Fokker-Planck problem has the form:

$\begin{array}[]{c}P_{t}\left(\tau,\theta,t,y\right)-\frac{1}{2}\varepsilon^{2}P_{yy}\left(\tau,\theta,t,y\right)+\left(\chi-\kappa
y\right)P_{y}\left(\tau,\theta,t,y\right)-\kappa
P\left(\tau,\theta,t,y\right)=0,\\\ \\\
P\left(\tau,\theta,t,y\right)=\delta\left(y-\theta\right).\end{array}$ (114)

The associated function $H\left(\tau,\theta,t,y,l\right)$ has the form:

$\begin{array}[]{c}H=\exp\left(\alpha\left(t\right)+i\Delta\left(t\right)y-il\theta\right),\end{array}$
(115)

so that

$\begin{array}[]{c}\alpha^{\prime}+\frac{1}{2}\varepsilon^{2}\Delta^{2}+i\chi\Delta-\kappa=0,\
\ \ \alpha\left(\tau\right)=0,\\\ \\\ \Delta^{\prime}-\kappa\Delta=0,\ \ \
\Delta\left(\tau\right)=l.\end{array}$ (116)

It is clear that

$\begin{array}[]{c}\Delta\left(\tau,t\right)=e^{\eta\left(\tau,t\right)}l,\\\
\\\
\alpha\left(\tau,t\right)=-\frac{1}{2}\psi_{0}\left(\tau,t\right)l^{2}-\left(\int_{\tau}^{t}e^{\eta\left(\tau,s\right)}\chi\left(s\right)ds\right)il+\eta\left(\tau,t\right).\end{array}$
(117)

where

$\begin{array}[]{c}\eta\left(\tau,t\right)=\mathop{\displaystyle\int}\limits_{\tau}^{t}\kappa\left(s\right)ds,\end{array}$
(118)

and

$\begin{array}[]{c}\psi_{0}\left(\tau,t\right)=\mathop{\displaystyle\int}\limits_{\tau}^{t}e^{2\eta\left(\tau,s\right)}\varepsilon^{2}\left(s\right)ds.\end{array}$
(119)

Thus, we obtain the following well-known expression:

$\begin{array}[]{c}P\left(\tau,\theta,t,y\right)=\frac{1}{2\pi}\int_{-\infty}^{\infty}\exp\left(-\frac{\psi_{0}\left(\tau,t\right)l^{2}}{2}\right.\\\
\\\
\left.+\left(e^{\eta\left(\tau,t\right)}y-\int_{\tau}^{t}e^{\eta\left(\tau,s\right)}\chi\left(s\right)ds-\theta\right)il+\eta\left(\tau,t\right)\right)dl\\\
\\\
=\frac{\exp\left(\eta\left(\tau,t\right)\right)}{\sqrt{2\pi\psi_{0}\left(\tau,t\right)}}\exp\left(-\frac{\left(e^{\eta\left(\tau,t\right)}y-\int_{\tau}^{t}e^{\eta\left(\tau,s\right)}\chi\left(s\right)ds-\theta\right)^{2}}{2\psi_{0}\left(\tau,t\right)}\right)\\\
\\\
=\frac{1}{\sqrt{2\pi\hat{\psi}_{0}\left(\tau,t\right)}}\exp\left(-\frac{\left(y-\int_{\tau}^{t}e^{-\eta\left(s,t\right)}\chi\left(s\right)ds-e^{-\eta\left(\tau,t\right)}\theta\right)^{2}}{2\hat{\psi}_{0}\left(\tau,t\right)}\right),\end{array}$
(120)

where

$\begin{array}[]{c}\hat{\psi}_{0}\left(\tau,t\right)=\frac{\psi_{0}\left(\tau,t\right)}{\exp\left(2\eta\left(\tau,t\right)\right)}=\mathop{\displaystyle\int}\limits_{\tau}^{t}e^{-2\eta\left(s,t\right)}\varepsilon^{2}\left(s\right)ds.\end{array}$
(121)

For further discussion, see the original paper by Uhlenbeck & Ornstein (1930),
as well as Chandresekhar (1943); Risken (1989), and references therein.

For time-independent parameters, Eq. (120) has the form:

$\begin{array}[]{c}P\left(\tau,\theta,t,y\right)=\frac{1}{\sqrt{\frac{\pi\varepsilon^{2}\left(1-e^{-2\kappa\left(t-\tau\right)}\right)}{\kappa}}}\exp\left(-\frac{\left(y-\frac{\chi}{\kappa}-e^{-\kappa\left(t-\tau\right)}\left(\theta-\frac{\chi}{\kappa}\right)\right)^{2}}{\frac{\varepsilon^{2}\left(1-e^{-2\kappa\left(t-\tau\right)}\right)}{\kappa}}\right).\end{array}$
(122)

####  5.3.2 Degenerate Augmented OU Process

In this subsection, we consider an augmented one-dimensional OU process:

$\begin{array}[]{c}dx_{t}=y_{t}dt,\ \ \ x_{\tau}=\xi,\\\ \\\
dy_{t}=\left(\chi\left(t\right)-\kappa\left(t\right)y_{t}\right)dt+\varepsilon\left(t\right)dW_{t},\
\ \ y_{\tau}=\theta.\end{array}$ (123)

To align our analysis with the existing body of work, we switch from the
general notation, used above, to a specific one customary used for the OU
process. Here and below, the word “augmentation” means that we expand the
original process by incorporating its integral or other path-dependent
characteristics, such as running maximum or minimum as part of the process.
The augmentation is very useful for handling large classes of path-dependent
options; details can be found in Lipton (2001), Chapter 13.

In the case in question, Eq. (92) can be written as follows:

$\begin{array}[]{c}\mathfrak{L}^{\prime}\left(\tau,t\right)+\left(\begin{array}[]{cc}0&0\\\
1&-\kappa\left(t\right)\end{array}\right)\mathfrak{L}\left(\tau,t\right)=0,\ \
\mathfrak{L}\left(\tau,\tau\right)=\left(\begin{array}[]{cc}1&0\\\
0&1\end{array}\right),\end{array}$ (124)

so that

$\begin{array}[]{c}\mathfrak{L}\left(\tau,t\right)=\left(\begin{array}[]{cc}1&0\\\
-\int_{\tau}^{t}e^{\eta\left(s,t\right)}ds&e^{\eta\left(\tau,t\right)}\end{array}\right),\\\
\\\ \mathfrak{L}^{-1}\left(\tau,t\right)=\left(\begin{array}[]{cc}1&0\\\
\int_{\tau}^{t}e^{-\eta\left(\tau,s\right)}ds&e^{-\eta\left(\tau,t\right)}\end{array}\right),\end{array}$
(125)

where $\eta$ is given by Eq. (118).

Now, we can compute $\mathfrak{C}^{-1}\left(\tau,t\right)$,
$\mathbf{d}^{\left(z\right)}\left(\tau,t\right)=\left(d^{\left(x\right)}\left(\tau,t\right),d^{\left(y\right)}\left(\tau,t\right)\right)^{T}$,
$\varpi\left(\tau,t\right)$:

$\begin{array}[]{c}\mathfrak{C}^{-1}\left(\tau,t\right)=\left(\begin{array}[]{cc}\psi_{2}\left(\tau,t\right)&-\psi_{1}\left(\tau,t\right)\\\
-\psi_{1}\left(\tau,t\right)&\psi_{0}\left(\tau,t\right)\end{array}\right),\end{array}$
(126)

where

$\begin{array}[]{c}\psi_{0}\left(\tau,t\right)=\mathop{\displaystyle\int}\limits_{\tau}^{t}e^{2\eta\left(\tau,s\right)}\varepsilon^{2}\left(s\right)ds,\\\
\\\
\psi_{1}\left(\tau,t\right)=-\mathop{\displaystyle\int}\limits_{\tau}^{t}\left(\mathop{\displaystyle\int}\limits_{\tau}^{s}e^{\eta\left(s^{\prime},t\right)}ds^{\prime}\right)e^{\eta\left(\tau,s\right)}\varepsilon^{2}\left(s\right)ds,\\\
\\\
\psi_{2}\left(\tau,t\right)=\mathop{\displaystyle\int}\limits_{\tau}^{t}\left(\mathop{\displaystyle\int}\limits_{\tau}^{s}e^{\eta\left(s^{\prime},t\right)}ds^{\prime}\right)^{2}\varepsilon^{2}\left(s\right)ds.\end{array}$
(127)
$\begin{array}[]{c}\mathbf{d}^{\left(z\right)}\left(\tau,t\right)=\left(\begin{array}[]{c}d^{\left(x\right)}\left(\tau,t\right)\\\
d^{\left(y\right)}\left(\tau,t\right)\end{array}\right)=\left(\begin{array}[]{c}-\int_{\tau}^{t}\left(\int_{\tau}^{s}e^{\eta\left(s^{\prime},s\right)}ds^{\prime}\right)\chi\left(s\right)ds\\\
\int_{\tau}^{t}e^{\eta\left(\tau,s\right)}\chi\left(s\right)ds,\end{array}\right),\end{array}$
(128)
$\begin{array}[]{c}\varpi\left(\tau,t\right)=-\eta\left(\tau,t\right).\end{array}$
(129)

Next, we calculate $\mathfrak{H}\left(\tau,t\right)$, and
$\mathbf{r}\left(\tau,t\right)$:

$\begin{array}[]{c}\mathfrak{H}\left(\tau,t\right)=\left(\mathfrak{L}^{T}\left(\tau,t\right)\right)^{-1}\mathfrak{C}^{-1}\left(\tau,t\right)\mathfrak{L}^{-1}\left(\tau,t\right)\\\
\\\
=\left(\begin{array}[]{cc}1&\int_{\tau}^{t}e^{-\eta\left(\tau,s\right)}ds\\\
0&e^{-\eta\left(\tau,t\right)}\end{array}\right)\left(\begin{array}[]{cc}\psi_{2}\left(\tau,t\right)&-\psi_{1}\left(\tau,t\right)\\\
-\psi_{1}\left(\tau,t\right)&\psi_{0}\left(\tau,t\right)\end{array}\right)\left(\begin{array}[]{cc}1&0\\\
\int_{\tau}^{t}e^{-\eta\left(\tau,s\right)}ds&e^{-\eta\left(\tau,t\right)}\end{array}\right)\\\
\\\
=\left(\begin{array}[]{cc}h_{0}\left(\tau,t\right)&h_{1}\left(\tau,t\right)\\\
h_{1}\left(\tau,t\right)&h_{2}\left(\tau,t\right)\end{array}\right),\end{array}$
(130)

where

$\begin{array}[]{c}h_{0}\left(\tau,t\right)=\psi_{0}\left(\int_{\tau}^{t}e^{-\eta\left(\tau,s\right)}ds\right)^{2}-2\psi_{1}\int_{\tau}^{t}e^{-\eta\left(\tau,s\right)}ds+\psi_{2},\\\
\\\
h_{1}\left(\tau,t\right)=\left(\psi_{0}\int_{\tau}^{t}e^{-\eta\left(\tau,s\right)}ds-\psi_{1}\right)e^{-\eta\left(\tau,t\right)},\\\
\\\
h_{2}\left(\tau,t\right)=\psi_{0}e^{-2\eta\left(\tau,t\right)}.\end{array}$
(131)
$\begin{array}[]{c}\mathbf{r}\left(\tau,t\right)=\left(\begin{array}[]{c}p\left(\tau,t\right)\\\
q\left(\tau,t\right)\end{array}\right)=\left(\mathfrak{L}^{T}\left(\tau,t\right)\right)^{-1}\left(\mathbf{d}^{\left(z\right)}\left(\tau,t\right)+\mathbf{\zeta}\right)\\\
\\\
=\left(\begin{array}[]{c}-\int_{\tau}^{t}\left(\int_{\tau}^{s}e^{\eta\left(s^{\prime},s\right)}ds^{\prime}\right)\chi\left(s\right)ds+\xi\\\
\\\
+\int_{\tau}^{t}e^{-\eta\left(\tau,s\right)}ds\left(\int_{\tau}^{t}e^{\eta\left(\tau,s\right)}\chi\left(s\right)ds+\theta\right),\\\
\\\
e^{-\eta\left(\tau,t\right)}\left(\int_{\tau}^{t}e^{\eta\left(\tau,s\right)}\chi\left(s\right)ds+\theta\right),\end{array}\right).\end{array}$
(132)

Thus, $P\left(\tau,\xi,\theta,t,x,y\right)$ is a bivariate Gaussian
distribution of the form (103) with the covariance matrix $\mathfrak{H}$,
given by Eq. (131) centered at the point $\mathbf{r=}\left(p,q\right)^{T}$
given by Eq. (132). Explicitly, we have

$\begin{array}[]{c}\sigma_{x}^{2}\left(\tau,t\right)=h_{0}\left(\tau,t\right),\
\ \ \sigma_{y}^{2}\left(\tau,t\right)=h_{2}\left(\tau,t\right),\ \ \
\rho\left(\tau,t\right)=\frac{h_{1}\left(\tau,t\right)}{\sqrt{h_{0}\left(\tau,t\right)h_{2}\left(\tau,t\right)}}.\end{array}$
(133)

When $\chi,\kappa,\varepsilon$ are constant, the above formulas become
significantly simpler. Namely,

$\begin{array}[]{c}\mathfrak{L}\left(\tau,t\right)=\left(\begin{array}[]{cc}1&0\\\
\frac{1-e^{\kappa\left(t-\tau\right)}}{\kappa}&e^{\kappa\left(t-\tau\right)}\end{array}\right),\\\
\\\ \mathfrak{L}^{-1}\left(\tau,t\right)=\left(\begin{array}[]{cc}1&0\\\
\frac{1-e^{-\kappa\left(t-\tau\right)}}{\kappa}&e^{-\kappa\left(t-\tau\right)}\end{array}\right).\end{array}$
(134)
$\begin{array}[]{c}\mathfrak{C}^{-1}\left(\tau,t\right)=\left(\begin{array}[]{cc}\psi_{2}\left(\tau,t\right)&-\psi_{1}\left(\tau,t\right)\\\
-\psi_{1}\left(\tau,t\right)&\psi_{0}\left(\tau,t\right)\end{array}\right),\\\
\\\
\psi_{0}\left(\tau,t\right)=\frac{\varepsilon^{2}\left(e^{2\kappa\left(t-\tau\right)}-1\right)}{2\kappa},\\\
\\\
\psi_{1}\left(\tau,t\right)=\frac{\varepsilon^{2}\left(e^{\kappa\left(t-\tau\right)}-1\right)^{2}}{2\kappa^{2}},\\\
\\\
\psi_{2}\left(\tau,t\right)=\frac{\varepsilon^{2}\left(e^{2\kappa\left(t-\tau\right)}-4e^{\kappa\left(t-\tau\right)}+3+2\kappa\left(t-\tau\right)\right)}{2\kappa^{3}},\end{array}$
(135)
$\begin{array}[]{c}\mathbf{d}^{\left(z\right)}\left(\tau,t\right)=\left(\begin{array}[]{c}d^{\left(x\right)}\left(\tau,t\right)\\\
d^{\left(y\right)}\left(\tau,t\right)\end{array}\right)=\left(\begin{array}[]{c}\left(\left(t-\tau\right)-\frac{e^{\kappa\left(t-\tau\right)}-1}{\kappa}\right)\frac{\chi}{\kappa}\\\
\left(e^{\kappa\left(t-\tau\right)}-1\right)\frac{\chi}{\kappa}\end{array}\right),\end{array}$
(136)
$\begin{array}[]{c}\varpi\left(\tau,t\right)=-\kappa\left(t-\tau\right),\end{array}$
(137)
$\begin{array}[]{c}\mathfrak{H}\left(\tau,t\right)=\left(\mathfrak{L}^{T}\left(\tau,t\right)\right)^{-1}\mathfrak{C}^{-1}\left(\tau,t\right)\mathfrak{L}^{-1}\left(\tau,t\right)=\left(\begin{array}[]{cc}h_{0}\left(\tau,t\right)&h_{1}\left(\tau,t\right)\\\
h_{1}\left(\tau,t\right)&h_{2}\left(\tau,t\right)\end{array}\right),\\\ \\\
h_{0}\left(\tau,t\right)=\frac{\varepsilon^{2}}{2\kappa^{3}}\left(-3+4e^{-\kappa\left(t-\tau\right)}-e^{-2\kappa\left(t-\tau\right)}+2\kappa\left(t-\tau\right)\right),\\\
\\\
h_{1}\left(\tau,t\right)=\frac{\varepsilon^{2}}{2\kappa^{2}}\left(1-e^{-\kappa\left(t-\tau\right)}\right)^{2},\\\
\\\
h_{2}\left(\tau,t\right)=\frac{\varepsilon^{2}\left(1-e^{-2\kappa\left(t-\tau\right)}\right)}{2\kappa}.\end{array}$
(138)
$\begin{array}[]{c}\mathbf{r}\left(\tau,t\right)=\left(\begin{array}[]{c}\xi+\frac{\left(\kappa\left(t-\tau\right)-\left(1-e^{-\kappa\left(t-\tau\right)}\right)\right)\chi}{\kappa^{2}}+\frac{\left(1-e^{-\kappa\left(t-\tau\right)}\right)\theta}{\kappa}\\\
\\\
\frac{\left(1-e^{-\kappa\left(t-\tau\right)}\right)\chi}{\kappa}+e^{-\kappa\left(t-\tau\right)}\theta\end{array}\right).\end{array}$
(139)

Thus, $P\left(\tau,\xi,\theta,t,x,y\right)$ is a bivariate Gaussian
distribution of the form (103) with the covariance matrix $\mathfrak{H}$,
given by Eq. (138) and the expected value $\mathbf{r=}\left(p,q\right)^{T}$
given by Eq. (139).

Below, we need to find the marginal distribution of $x$, which we denote by
$Q^{\left(x\right)}\left(\tau,\theta,t,x\right)$. It is well known that
marginal distributions of a multivariate Gaussian distribution are also
Gaussian, so that:

$\begin{array}[]{c}Q^{\left(x\right)}\left(\tau,\theta,t,x\right)=\frac{1}{\sqrt{2\pi
h_{0}\left(\tau,t\right)}}\exp\left(\frac{\left(x-p\left(\tau,t\right)\right)^{2}}{2h_{0}\left(\tau,t\right)}\right),\end{array}$
(140)

where $h_{0}$ is given by Eqs (138). At the same time, the density of marginal
distribution for $y$ has the form:

$\begin{array}[]{c}Q^{\left(y\right)}\left(\tau,\theta,t,y\right)=\frac{1}{\sqrt{2\pi
h_{2}\left(\tau,t\right)}}\exp\left(\frac{\left(y-q\left(\tau,t\right)\right)^{2}}{2h_{2}\left(\tau,t\right)}\right),\end{array}$
(141)

where $h_{2}$ is given by Eqs (138), which is the familiar density of the OU
process derived in the previous section.

###  5.4 Example: Evolution of Free and Harmonically Bound Particles

We can use the above results to revisit the motion of free and harmonically
bound particles considered in Section 3.

To describe a free particle, we assume that $\chi=0$. Eq. (138) does not
change, while Eq. (139) becomes a bit simpler:

$\begin{array}[]{c}\left(p\left(\tau,t\right),q\left(\tau,t\right)\right)^{T}=\left(\xi+\frac{\left(1-e^{-\kappa\left(t-\tau\right)}\right)\theta}{\kappa},e^{-\kappa\left(t-\tau\right)}\theta\right)^{T}.\end{array}$
(142)

It is clear that Eqs (21), (22) and (138), (142) are in agreement.

Analysis of a harmonically bound particle requires additional efforts. In the
case in question, Eq. (92) can be written as follows:

$\begin{array}[]{c}\mathfrak{L}^{\prime}\left(\tau,t\right)+\left(\begin{array}[]{cc}0&-\omega^{2}\\\
1&-\kappa\end{array}\right)\mathfrak{L}\left(\tau,t\right)=0,\ \
\mathfrak{L}\left(\tau,\tau\right)=\left(\begin{array}[]{cc}1&0\\\
0&1\end{array}\right).\end{array}$ (143)

The corresponding second-order equation is

$\begin{array}[]{c}l^{\prime\prime}-\kappa
l^{\prime}+\omega^{2}l=0,\end{array}$ (144)

while the characteristic equation and its solutions can be written as follows:

$\begin{array}[]{c}\lambda^{2}-\kappa\lambda+\omega^{2}=0,\\\ \\\
\lambda_{\pm}=\frac{\kappa\pm\sqrt{\kappa^{2}-4\omega^{2}}}{2}=-\mu_{\mp},\end{array}$
(145)

where $\mu_{\pm}$ are given by Eq. (25). It is easy to check that

$\begin{array}[]{c}\mathfrak{L}\left(\tau,t\right)=\frac{1}{\sqrt{\kappa^{2}-4\omega^{2}}}\left(\begin{array}[]{cc}-\left(\lambda_{-}\mathsf{E}_{+}-\lambda_{+}\mathsf{E}_{-}\right)&\omega^{2}\left(\mathsf{E}_{+}-\mathsf{E}_{-}\right)\\\
-\left(\mathsf{E}_{+}-\mathsf{E}_{-}\right)&\left(\lambda_{+}\mathsf{E}_{+}-\lambda_{-}\mathsf{E}_{-}\right)\end{array}\right),\end{array}$
(146)
$\begin{array}[]{c}\mathfrak{L}^{-1}\left(\tau,t\right)=\frac{1}{\sqrt{\kappa^{2}-4\omega^{2}}}\left(\begin{array}[]{cc}\left(\lambda_{+}\mathsf{E}_{-}^{-1}-\lambda_{-}\mathsf{E}_{+}^{-1}\right)&-\omega^{2}\left(\mathsf{E}_{-}^{-1}-\mathsf{E}_{+}^{-1}\right)\\\
\left(\mathsf{E}_{-}^{-1}-\mathsf{E}_{+}^{-1}\right)&-\left(\lambda_{-}\mathsf{E}_{-}^{-1}-\lambda_{+}\mathsf{E}_{+}^{-1}\right)\end{array}\right),\end{array}$
(147)
$\begin{array}[]{c}\det\mathfrak{L=}\left(\det\mathfrak{L}^{-1}\right)^{-1}\mathfrak{=}\mathsf{E}_{+}\mathsf{E}_{-}=\exp\left(\kappa\left(t-\tau\right)\right),\end{array}$
(148)

where

$\begin{array}[]{c}\mathsf{E}_{+}=e^{\lambda_{\pm}\left(t-\tau\right)}.\end{array}$
(149)

Accordingly,

$\begin{array}[]{c}P\left(\tau,\mathbf{\zeta,}t,\mathbf{z}\right)=\mathrm{N}\left(\mathbf{r}\left(\tau,t\right),\mathfrak{H}\left(\tau,t\right)\right),\end{array}$
(150)

with

$\begin{array}[]{c}\mathfrak{H}\left(\tau,t\right)=\left(\mathfrak{L}^{T}\left(\tau,t\right)\right)^{-1}\mathfrak{C}^{-1}\left(\tau,t\right)\mathfrak{L}^{-1}\left(\tau,t\right),\\\
\\\
\mathbf{r}\left(\tau,t\right)\mathbf{=}\left(\mathfrak{L}^{T}\left(\tau,t\right)\right)^{-1}\mathbf{\zeta}.\end{array}$
(151)

Here

$\begin{array}[]{c}\mathbf{r=}\left(p,q\right)^{T}=\left(\frac{\left(\left(\lambda_{+}\mathsf{E}_{-}^{-1}-\lambda_{-}\mathsf{E}_{+}^{-1}\right)\xi+\left(\mathsf{E}_{-}^{-1}-\mathsf{E}_{+}^{-1}\right)\theta\right)}{\sqrt{\kappa^{2}-4\omega^{2}}},-\frac{\left(\omega^{2}\left(\mathsf{E}_{-}^{-1}-\mathsf{E}_{+}^{-1}\right)\xi+\left(\lambda_{-}\mathsf{E}_{-}^{-1}-\lambda_{+}\mathsf{E}_{+}^{-1}\right)\theta\right)}{\sqrt{\kappa^{2}-4\omega^{2}}}\right)^{T}.\end{array}$
(152)

In the limit $\omega^{2}\rightarrow 0$,

$\begin{array}[]{c}\mathbf{r=}\left(p,q\right)^{T}=\left(\xi+\frac{\left(1-e^{-\kappa\left(t-\tau\right)}\right)\theta}{\kappa},e^{-\kappa\left(t-\tau\right)}\theta\right)^{T}.\end{array}$
(153)

Next,

$\begin{array}[]{c}\mathfrak{C}^{-1}\left(\tau,t\right)=\varepsilon^{2}\left(\begin{array}[]{cc}\int_{\tau}^{t}\mathfrak{L}^{\left(yx\right)T}\left(s\right)\mathfrak{L}^{\left(yx\right)}\left(s\right)ds&\int_{\tau}^{t}\mathfrak{L}^{\left(yx\right)T}\left(s\right)\mathfrak{L}^{\left(yy\right)}\left(s\right)ds\\\
\int_{\tau}^{t}\mathfrak{L}^{\left(yy\right)T}\left(s\right)\mathfrak{L}^{\left(yx\right)}\left(s\right)ds&\int_{\tau}^{t}\mathfrak{L}^{\left(yy\right)T}\left(s\right)\mathfrak{L}^{\left(yy\right)}\left(s\right)ds\end{array}\right)\\\
\\\
=\left(\begin{array}[]{cc}\psi_{2}\left(\tau,t\right)&-\psi_{1}\left(\tau,t\right)\\\
-\psi_{1}\left(\tau,t\right)&\psi_{0}\left(\tau,t\right)\end{array}\right)\equiv\frac{\varepsilon^{2}}{\left(\kappa^{2}-4\omega^{2}\right)}\left(\begin{array}[]{cc}\bar{\psi}_{2}\left(\tau,t\right)&-\bar{\psi}_{1}\left(\tau,t\right)\\\
-\bar{\psi}_{1}\left(\tau,t\right)&\bar{\psi}_{0}\left(\tau,t\right)\end{array}\right),\end{array}$
(154)

where

$\begin{array}[]{c}\bar{\psi}_{0}\left(\tau,t\right)=\int_{\tau}^{t}\left(\lambda_{+}\mathsf{E}_{+}\left(s\right)-\lambda_{-}\mathsf{E}_{-}\left(s\right)\right)^{2}ds\\\
\\\
=\frac{1}{2\kappa}\left(\kappa\lambda_{+}\mathsf{E}_{+}^{2}\left(t\right)-4\omega^{2}\mathsf{E}_{+}\left(t\right)\mathsf{E}_{-}\left(t\right)+\kappa\lambda_{-}\mathsf{E}_{-}^{2}\left(t\right)-\left(\kappa^{2}-4\omega^{2}\right)\right),\end{array}$
(155)
$\begin{array}[]{c}\bar{\psi}_{1}\left(\tau,t\right)=\int_{\tau}^{t}\left(\mathsf{E}_{+}\left(s\right)-\mathsf{E}_{-}\left(s\right)\right)\left(\lambda_{+}\mathsf{E}_{+}\left(s\right)-\lambda_{-}\mathsf{E}_{-}\left(s\right)\right)ds\\\
\\\
=\frac{1}{2}\left(\mathsf{E}_{+}^{2}\left(t\right)-2\mathsf{E}_{+}\left(t\right)\mathsf{E}_{-}\left(t\right)+\mathsf{E}_{-}^{2}\left(t\right)\right),\end{array}$
(156)
$\begin{array}[]{c}\bar{\psi}_{2}\left(\tau,t\right)=\int_{\tau}^{t}\left(\mathsf{E}_{+}\left(s\right)-\mathsf{E}_{-}\left(s\right)\right)^{2}ds\\\
\\\
=\frac{1}{2\kappa\omega^{2}}\left(\kappa\lambda_{-}\mathsf{E}_{+}^{2}\left(t\right)-4\omega^{2}\mathsf{E}_{+}\left(t\right)\mathsf{E}_{-}\left(t\right)+\kappa\lambda_{+}\mathsf{E}_{-}^{2}\left(t\right)-\left(\kappa^{2}-4\omega^{2}\right)\right).\end{array}$
(157)

Further,

$\begin{array}[]{c}\mathfrak{H}\left(\tau,t\right)=\left(\mathfrak{L}^{T}\left(\tau,t\right)\right)^{-1}\mathfrak{C}^{-1}\left(\tau,t\right)\mathfrak{L}^{-1}\left(\tau,t\right)\\\
\\\
=\left(\begin{array}[]{cc}h_{0}\left(\tau,t\right)&h_{1}\left(\tau,t\right)\\\
h_{1}\left(\tau,t\right)&h_{2}\left(\tau,t\right)\end{array}\right)\equiv\frac{\varepsilon^{2}}{\left(\kappa^{2}-4\omega^{2}\right)^{2}}\left(\begin{array}[]{cc}\bar{h}_{0}\left(\tau,t\right)&\bar{h}_{1}\left(\tau,t\right)\\\
\bar{h}_{1}\left(\tau,t\right)&\bar{h}_{2}\left(\tau,t\right)\end{array}\right),\end{array}$
(158) $\begin{array}[]{c}\left(\begin{array}[]{cc}\bar{h}_{0}&\bar{h}_{1}\\\
\bar{h}_{1}&\bar{h}_{2}\end{array}\right)=\left(\begin{array}[]{cc}\left(\lambda_{+}\mathsf{E}_{-}^{-1}-\lambda_{-}\mathsf{E}_{+}^{-1}\right)&\left(\mathsf{E}_{-}^{-1}-\mathsf{E}_{+}^{-1}\right)\\\
-\omega^{2}\left(\mathsf{E}_{-}^{-1}-\mathsf{E}_{+}^{-1}\right)&-\left(\lambda_{-}\mathsf{E}_{-}^{-1}-\lambda_{+}\mathsf{E}_{+}^{-1}\right)\end{array}\right)\\\
\\\ \times\left(\begin{array}[]{cc}\bar{\psi}_{2}&-\bar{\psi}_{1}\\\
-\bar{\psi}_{1}&\bar{\psi}_{0}\end{array}\right)\left(\begin{array}[]{cc}\left(\lambda_{+}\mathsf{E}_{-}^{-1}-\lambda_{-}\mathsf{E}_{+}^{-1}\right)&-\omega^{2}\left(\mathsf{E}_{-}^{-1}-\mathsf{E}_{+}^{-1}\right)\\\
\left(\mathsf{E}_{-}^{-1}-\mathsf{E}_{+}^{-1}\right)&-\left(\lambda_{-}\mathsf{E}_{-}^{-1}-\lambda_{+}\mathsf{E}_{+}^{-1}\right)\end{array}\right).\end{array}$
(159)

Straightforward but very tedious calculation yields

$\begin{array}[]{c}\bar{h}_{0}=\frac{\left(\kappa^{2}-4\omega^{2}\right)}{2\kappa\omega^{2}}\left(-\kappa\lambda_{-}\mathsf{E}_{+}^{-2}+4\omega^{2}\mathsf{E}_{-}^{-1}\mathsf{E}_{+}^{-1}-\kappa\lambda_{+}\mathsf{E}_{-}^{-2}+\kappa^{2}-4\omega^{2}\right),\\\
\\\
\bar{h}_{1}=\frac{\left(\kappa^{2}-4\omega^{2}\right)}{2}\left(\mathsf{E}_{-}^{-2}-2\mathsf{E}_{+}^{-1}\mathsf{E}_{-}^{-1}+\mathsf{E}_{+}^{-2}\right),\\\
\\\
\bar{h}_{2}=\frac{\left(\kappa^{2}-4\omega^{2}\right)}{2\kappa}\left(-\kappa\lambda_{+}\mathsf{E}_{+}^{-2}+4\omega^{2}\mathsf{E}_{-}^{-1}\mathsf{E}_{+}^{-1}-\kappa\lambda_{-}\mathsf{E}_{-}^{-2}+\kappa^{2}-4\omega^{2}\right),\end{array}$
(160)

and

$\begin{array}[]{c}h_{0}=\frac{\varepsilon^{2}}{2\kappa\omega^{2}\left(\kappa^{2}-4\omega^{2}\right)}\left(-\kappa\lambda_{-}\mathsf{E}_{+}^{-2}+4\omega^{2}\mathsf{E}_{-}^{-1}\mathsf{E}_{+}^{-1}-\kappa\lambda_{+}\mathsf{E}_{-}^{-2}+\kappa^{2}-4\omega^{2}\right),\\\
\\\
h_{1}=\frac{\varepsilon^{2}}{2\left(\kappa^{2}-4\omega^{2}\right)}\left(\mathsf{E}_{-}^{-2}-2\mathsf{E}_{+}^{-1}\mathsf{E}_{-}^{-1}+\mathsf{E}_{+}^{-2}\right),\\\
\\\
h_{2}=\frac{\varepsilon^{2}}{2\kappa\left(\kappa^{2}-4\omega^{2}\right)}\left(-\kappa\lambda_{+}\mathsf{E}_{+}^{-2}+4\omega^{2}\mathsf{E}_{-}^{-1}\mathsf{E}_{+}^{-1}-\kappa\lambda_{-}\mathsf{E}_{-}^{-2}+\kappa^{2}-4\omega^{2}\right).\end{array}$
(161)

In the limit $\omega^{2}\rightarrow 0$,

$\begin{array}[]{c}h_{0}=\frac{\varepsilon^{2}}{2\kappa^{3}}\left(-3+2\kappa\left(t-\tau\right)+4e^{-\kappa\left(t-\tau\right)}-e^{-2\kappa\left(t-\tau\right)}\right),\\\
\\\
h_{1}=\frac{\varepsilon^{2}}{2\kappa^{2}}\left(\mathsf{1}-2e^{-\kappa\left(t-\tau\right)}+e^{-2\kappa\left(t-\tau\right)}\right),\\\
\\\
h_{2}=\frac{\varepsilon^{2}}{2\kappa}\left(1-e^{-2\kappa\left(t-\tau\right)}\right),\end{array}$
(162)

so that Eqs. (138) and (162) are in agreement.

Moreover, while it is easy to show that Chandrasekhar’s solution given by Eq.
(28) is in agreement with our solution given by Eq. (150), our solution is
more convenient from a practical standpoint, since it is _explicitly_ written
as a Gaussian density in the $\left(x,y\right)$ space.

###  5.5 Example: Evolution of Vorticity of Two-Dimensional Flows

We shall briefly return to our point of departure and consider strictly two-
dimensional flows; see Friedlander & Lipton-Lifschitz (2003). Velocity fields
of such flows have the form:

$\begin{array}[]{c}\mathbf{V}\left(t,x_{1},x_{2}\right)=\left(V_{1}\left(t,x_{1},x_{2}\right),V_{2}\left(t,x_{1},x_{2}\right)\right),\\\
\\\
\mathbf{v}\left(t,x_{1},x_{2}\right)=\left(v_{1}\left(t,x_{1},x_{2}\right),v_{2}\left(t,x_{1},x_{2}\right)\right).\end{array}$
(163)

By virtue of incompressibility, we can introduce the so-called stream
functions such that

$\begin{array}[]{c}V_{1}=-\frac{\partial\Psi}{\partial x_{2}},\ \ \
V_{2}=\frac{\partial\Psi}{\partial x_{1}},\ \ \
v_{1}=-\frac{\partial\psi}{\partial x_{2}},\ \ \
v_{2}=\frac{\partial\psi}{\partial x_{1}},\end{array}$ (164)

and define the (scalar) vorticity as follows:

$\begin{array}[]{c}\Omega=\Delta\Psi,\ \ \ \omega=\Delta\psi.\end{array}$
(165)

Contour lines of $\Psi$ are called streamlines of the flow.

By using the above definitions, we write the two-dimensional the Navier-Stokes
equations as equations for the stream and vorticity:

$\begin{array}[]{c}\frac{\partial\Omega}{\partial
t}-\frac{\partial\Psi}{\partial x_{2}}\frac{\partial\Omega}{\partial
x_{1}}+\frac{\partial\Psi}{\partial x_{1}}\frac{\partial\Omega}{\partial
x_{2}}-\nu\Delta\Omega=0,\\\ \\\ \Delta\Psi-\Omega=0.\end{array}$ (166)

It is clear that time-independent stream functions $\Psi$, quadratically
dependent of $\left(x_{1},x_{2}\right)$ generate exact solutions of Eqs (166).
We are interested in fields consisting of pure strain and pure rotation. The
corresponding $\Psi$ have the form:

$\begin{array}[]{c}\Psi\left(x_{1},x_{2}\right)=\frac{1}{4}\left(w\left(x_{1}^{2}+x_{2}^{2}\right)-2sx_{1}x_{2}\right),\end{array}$

where $\omega>s$, to ensure that streamlines are elliptic rather than
hyperbolic, so that

$\begin{array}[]{c}V_{1}=-\frac{\partial\Psi}{\partial
x_{2}}=\frac{1}{2}\left(sx_{1}-wx_{2}\right),\ \ \
V_{2}=\frac{\partial\Psi}{\partial
x_{1}}=\frac{1}{2}\left(wx_{1}-sx_{2}\right).\end{array}$ (167)

Small perturbations $\psi$ of the time-independent quadratic stream function
$\Psi$ satisfy the following equations:

$\begin{array}[]{c}\frac{\partial\omega}{\partial
t}-\frac{\partial\Psi}{\partial x_{2}}\frac{\partial\omega}{\partial
x_{1}}+\frac{\partial\Psi}{\partial x_{1}}\frac{\partial\omega}{\partial
x_{2}}-\nu\Delta\omega=0,\\\ \\\ \Delta\psi-\omega=0.\end{array}$ (168)

We wish to study the first equation (168) in isolation, which we write
explicitly as follows:

$\begin{array}[]{c}\frac{\partial\omega}{\partial
t}+\frac{1}{2}\left(sx_{1}-wx_{2}\right)\frac{\partial\omega}{\partial
x_{1}}+\frac{1}{2}\left(wx_{1}-sx_{2}\right)\frac{\partial\omega}{\partial
x_{2}}-\nu\Delta\omega=0,\end{array}$ (169)

and supply with the initial condition at time $\tau$:

$\begin{array}[]{c}\omega\left(\tau,x_{1},x_{2}\right)=\delta\left(x_{1}-\xi_{1}\right)\delta\left(x_{2}-\xi_{2}\right).\end{array}$
(170)

Once the solution of Eqs (169), (170) is found, we can find $\psi$ by solving
the corresponding Laplace equation.

Surprisingly, we can view this equation as the Fokker-Planck equation
associated with the following SDEs for
$\mathbf{z}_{t}=\left(x_{1t},x_{2t}\right)$:

$\begin{array}[]{c}d\mathbf{z}_{t}=\mathfrak{B}\mathbf{z}_{t}dt+\mathbf{\Sigma}d\mathbf{W}_{t},\
\ \ \mathbf{z}_{\tau}=\left(\xi_{1},\xi_{2}\right)^{T},\end{array}$ (171)

where

$\begin{array}[]{c}\mathfrak{B=}\frac{1}{2}\left(\begin{array}[]{cc}s&-w\\\
w&-s\end{array}\right),\ \ \
\mathbf{\Sigma=}\sqrt{2\nu}\left(\begin{array}[]{cc}1&0\\\
0&1\end{array}\right).\end{array}$ (172)

Thus, we can use Section 5.1 results. Eq. (92) becomes

$\begin{array}[]{c}\mathfrak{L}^{\prime}\left(\tau,t\right)+\frac{1}{2}\left(\begin{array}[]{cc}s&w\\\
-w&-s\end{array}\right)\mathfrak{L}\left(\tau,t\right)=0,\ \
\mathfrak{L}\left(\tau,\tau\right)=\left(\begin{array}[]{cc}1&0\\\
0&1\end{array}\right).\end{array}$ (173)

Simple but tedious calculations omitted for the sake of brevity show that

$\begin{array}[]{c}\mathfrak{L}\left(\tau,t\right)=\left(\begin{array}[]{cc}\mathsf{C}_{1}\left(\tau,t\right)-\frac{s}{2\zeta}\mathsf{S}_{1}\left(\tau,t\right)&-\frac{w}{2\zeta}\mathsf{S}_{1}\left(\tau,t\right)\\\
\frac{w}{2\zeta}\mathsf{S}_{1}\left(\tau,t\right)&\mathsf{C}_{1}\left(\tau,t\right)+\frac{s}{2\zeta}\mathsf{S}_{1}\left(\tau,t\right)\end{array}\right),\end{array}$
(174)

and

$\begin{array}[]{c}\mathfrak{L}^{-1}\left(\tau,t\right)=\left(\begin{array}[]{cc}\mathsf{C}_{1}\left(\tau,t\right)+\frac{s}{2\zeta}\mathsf{S}_{1}\left(\tau,t\right)&\frac{w}{2\zeta}\mathsf{S}_{1}\left(\tau,t\right)\\\
-\frac{w}{2\zeta}\mathsf{S}_{1}\left(\tau,t\right)&\mathsf{C}_{1}\left(\tau,t\right)-\frac{s}{2\zeta}\mathsf{S}_{1}\left(\tau,t\right)\end{array}\right),\end{array}$
(175)

where

$\begin{array}[]{c}\zeta=\frac{\sqrt{w^{2}-s^{2}}}{2},\ \ \
\mathsf{C}_{1}\left(\tau,t\right)=\cos\left(\zeta\left(t-\tau\right)\right),\
\ \
\mathsf{S}_{1}\left(\tau,t\right)=\sin\left(\zeta\left(t-\tau\right)\right),\\\
\\\
\det\left(\mathfrak{L}\right)=\det\left(\mathfrak{L}^{-1}\right)=1.\end{array}$
(176)

Next, Eq. (97) yields

$\begin{array}[]{c}\mathfrak{C}^{-1}\left(\tau,t\right)=\int_{\tau}^{t}\mathfrak{L}^{T}\left(\tau,s\right)\mathfrak{A}\left(s\right)\mathfrak{L}\left(\tau,s\right)ds\\\
\\\
=2\nu\int_{\tau}^{t}\mathfrak{L}^{T}\left(\tau,s\right)\mathfrak{L}\left(\tau,s\right)ds=\left(\begin{array}[]{cc}\psi_{2}\left(\tau,t\right)&-\psi_{1}\left(\tau,t\right)\\\
-\psi_{1}\left(\tau,t\right)&\psi_{0}\left(\tau,t\right)\end{array}\right),\end{array}$
(177)

where

$\begin{array}[]{c}\psi_{0}\left(\tau,t\right)=2\nu\int_{\tau}^{t}\left(1+\frac{s}{2\zeta}\mathsf{S}_{2}\left(\tau,s\right)+\frac{s^{2}}{4\zeta^{2}}\left(1-\mathsf{C}_{2}\left(\tau,s\right)\right)\right)ds\\\
\\\
=2\nu\left(\left(1+\frac{s^{2}}{4\zeta^{2}}\right)\left(t-\tau\right)-\frac{s}{4\zeta^{2}}\mathsf{C}_{2}\left(\tau,t\right)-\frac{s^{2}}{8\zeta^{3}}\mathsf{S}_{2}\left(\tau,t\right)\right),\end{array}$
(178) $\begin{array}[]{c}\psi_{1}\left(\tau,t\right)=-\frac{\nu
sw}{2\zeta^{2}}\int_{\tau}^{t}\left(1-\mathsf{C}_{2}\left(\tau,s\right)\right)ds\\\
\\\ =-\frac{\nu
sw}{2\zeta^{2}}\left(\left(t-\tau\right)-\frac{1}{2\zeta}\mathsf{S}_{2}\left(\tau,t\right)\right),\end{array}$
(179)
$\begin{array}[]{c}\psi_{2}\left(\tau,t\right)=2\nu\int_{\tau}^{t}\left(1-\frac{s}{2\zeta}\mathsf{S}_{2}\left(\tau,s\right)+\frac{s^{2}}{4\zeta^{2}}\left(1-\mathsf{C}_{2}\left(\tau,s\right)\right)\right)ds\\\
\\\
=2\nu\left(\left(1+\frac{s^{2}}{4\zeta^{2}}\right)\left(t-\tau\right)+\frac{s}{4\zeta^{2}}\mathsf{C}_{2}\left(\tau,t\right)-\frac{s^{2}}{8\zeta^{3}}\mathsf{S}_{2}\left(\tau,t\right)\right),\end{array}$
(180)

where

$\begin{array}[]{c}\mathsf{C}_{2}\left(\tau,t\right)=\cos\left(2\zeta\left(t-\tau\right)\right),\
\ \
\mathsf{S}_{2}\left(\tau,t\right)=\sin\left(2\zeta\left(t-\tau\right)\right).\end{array}$
(181)

Finally, Eqs (103), (104) yield:

$\begin{array}[]{c}\omega\left(\tau,\mathbf{\zeta,}t,\mathbf{z}\right)=\mathrm{N}\left(\mathbf{r}\left(\tau,t\right),\mathfrak{H}\left(\tau,t\right)\right),\end{array}$
(182)
$\begin{array}[]{c}\mathfrak{H}\left(\tau,t\right)=\left(\mathfrak{L}^{T}\left(\tau,t\right)\right)^{-1}\mathfrak{C}^{-1}\left(\tau,t\right)\mathfrak{L}^{-1}\left(\tau,t\right)=\left(\begin{array}[]{cc}h_{0}\left(\tau,t\right)&h_{1}\left(\tau,t\right)\\\
h_{1}\left(\tau,t\right)&h_{2}\left(\tau,t\right)\end{array}\right),\end{array}$
(183)

where

$\begin{array}[]{c}h_{0}\left(\tau,t\right)=\left(\frac{w^{2}}{8\zeta^{2}}+\frac{4\zeta^{2}-s^{2}}{8\zeta^{2}}\mathsf{C}_{2}\left(\tau,t\right)+\frac{s}{2\zeta}\mathsf{S}_{2}\left(\tau,t\right)\right)\psi_{2}\left(\tau,t\right)\\\
\\\
+\frac{w}{2\zeta}\left(\frac{s}{2\zeta}\left(1-\mathsf{C}_{2}\left(\tau,t\right)\right)+\mathsf{S}_{2}\left(\tau,t\right)\right)\psi_{1}\left(\tau,t\right)+\frac{w^{2}}{8\zeta^{2}}\left(1-\mathsf{C}_{2}\left(\tau,t\right)\right)\psi_{0}\left(\tau,t\right),\\\
\\\
h_{1}\left(\tau,t\right)=\frac{w}{4\zeta}\left(\frac{s}{2\zeta}\left(1-\mathsf{C}_{2}\left(\tau,t\right)\right)+\mathsf{S}_{2}\left(\tau,t\right)\right)\psi_{2}\left(\tau,t\right)\\\
\\\
-\left(1+\frac{w^{2}}{4\zeta^{2}}\left(1-\mathsf{C}_{2}\left(\tau,t\right)\right)\right)\psi_{1}\left(\tau,t\right)+\frac{w}{4\zeta}\left(\frac{s}{2\zeta}\left(1-\mathsf{C}_{2}\left(\tau,t\right)\right)-\mathsf{S}_{2}\left(\tau,t\right)\right)\psi_{0}\left(\tau,t\right),\\\
\\\
h_{2}\left(\tau,t\right)=\frac{w^{2}}{8\zeta^{2}}\left(1-\mathsf{C}_{2}\left(\tau,t\right)\right)\psi_{2}\left(\tau,t\right)+\frac{w}{2\zeta}\left(\frac{s}{2\zeta}\left(1-\mathsf{C}_{2}\left(\tau,t\right)\right)-\mathsf{S}_{2}\left(\tau,t\right)\right)\psi_{1}\left(\tau,t\right)\\\
\\\
+\left(\frac{w^{2}}{8\zeta^{2}}+\frac{4\zeta^{2}-s^{2}}{8\zeta^{2}}\mathsf{C}_{2}\left(\tau,t\right)-\frac{s}{2\zeta}\mathsf{S}_{2}\left(\tau,t\right)\right)\psi_{0}\left(\tau,t\right),\end{array}$
(184)

and

$\begin{array}[]{c}\mathbf{r}\left(\tau,t\right)\mathbf{=}\left(\mathfrak{L}^{T}\left(\tau,t\right)\right)^{-1}\mathbf{\zeta=}\left(\begin{array}[]{c}r_{1}\left(\tau,t\right)\\\
r_{2}\left(\tau,t\right)\end{array}\right)\\\ \\\
=\left(\begin{array}[]{c}\left(\mathsf{C}_{1}\left(\tau,t\right)+\frac{s}{2\zeta}\mathsf{S}_{1}\left(\tau,t\right)\right)\xi_{1}-\frac{w}{2\zeta}\mathsf{S}_{1}\left(\tau,t\right)\xi_{2}\\\
\frac{w}{2\zeta}\mathsf{S}_{1}\left(\tau,t\right)\xi_{1}+\left(\mathsf{C}_{1}\left(\tau,t\right)-\frac{s}{2\zeta}\mathsf{S}_{1}\left(\tau,t\right)\right)\xi_{2}\end{array}\right).\end{array}$
(185)

Eqs (184) are symmetric, namely $h_{0}\rightarrow h_{2}$ when
$\left(a,b\right)\rightarrow\left(-a,-b\right)$,
$\left(\psi_{0},\psi_{2}\right)\rightarrow\left(\psi_{2},\psi_{0}\right)$. The
second Eq. (168), which is a static Poisson equation, allows us to find
$\psi$. Its analytical solution is not easy to derive and we shall not pursue
this line of research here. The special case of purely rotational flow, $s=0$,
can be done easily; see Eq. (193) below.

It is interesting to note that

$\begin{array}[]{c}\Psi\left(r_{1},r_{2}\right)=\Psi\left(\xi_{1},\xi_{2}\right),\end{array}$
(186)

so that the location of the Gaussian distribution $\omega$ moves along
streamlines of the flow defined by the stream function $\Psi$.

When the flow is purely rotational, so that $s=0$, the above formulas
considerably simplify. Specifically, we have:

$\begin{array}[]{c}\psi_{0}\left(\tau,t\right)=2\nu\left(t-\tau\right),\ \
\psi_{1}\left(\tau,t\right)=0,\ \ \ \ \
\psi_{2}\left(\tau,t\right)=2\nu\left(t-\tau\right),\\\ \\\
h_{0}\left(\tau,t\right)=2\nu\left(t-\tau\right),\ \
\psi_{1}\left(\tau,t\right)=0,\ \ \ \ \
h_{2}\left(\tau,t\right)=2\nu\left(t-\tau\right),\\\ \\\
r_{1}\left(\tau,t\right)=\mathsf{C}_{1}\left(\tau,t\right)\xi_{1}-\mathsf{S}_{1}\left(\tau,t\right)\xi_{2},\
\ \
r_{2}\left(\tau,t\right)=\mathsf{S}_{1}\left(\tau,t\right)\xi_{1}+\mathsf{C}_{1}\left(\tau,t\right)\xi_{2},\end{array}$
(187)

so that

$\begin{array}[]{c}\omega\left(\tau,\xi_{1},\xi_{2},t,x_{1},x_{2}\right)\\\
\\\
=\frac{1}{4\pi\nu\left(t-\tau\right)}\exp\left(-\frac{\left(x_{1}-\mathsf{C}_{1}\left(\tau,t\right)\xi_{1}+\mathsf{S}_{1}\left(\tau,t\right)\xi_{2}\right)^{2}+\left(x_{2}-\mathsf{S}_{1}\left(\tau,t\right)\xi_{1}-\mathsf{C}_{1}\left(\tau,t\right)\xi_{2}\right)^{2}}{4\nu\left(t-\tau\right)}\right).\end{array}$
(188)

We can directly calculate the stream function $\psi$ by solving the
corresponding Poisson equation.222We are grateful to Andrey Itkin for pointing
this out. To start, we notice that _both_ $\omega$ and $\psi$ are rotationally
symmetric around the point $\left(\xi_{1},\xi_{2}\right)$. Thus, we can
introduce

$\begin{array}[]{c}R^{2}=\frac{\left(x_{1}-\mathsf{C}_{1}\left(\tau,t\right)\xi_{1}+\mathsf{S}_{1}\left(\tau,t\right)\xi_{2}\right)^{2}+\left(x_{2}-\mathsf{S}_{1}\left(\tau,t\right)\xi_{1}-\mathsf{C}_{1}\left(\tau,t\right)\xi_{2}\right)^{2}}{2\nu\left(t-\tau\right)},\end{array}$
(189)

representing $\omega$ and $\psi$ as follows:

$\begin{array}[]{c}\omega=\omega\left(R\right)=\frac{1}{4\pi\nu\left(t-\tau\right)}\exp\left(-\frac{R^{2}}{2}\right),\
\ \ \psi=\psi\left(R\right).\end{array}$ (190)

Then $\psi\left(R\right)$ solves a radially symmetric Poisson equation of the
form:

$\begin{array}[]{c}\frac{1}{R}\left(R\psi_{R}\left(R\right)\right)_{R}=\frac{1}{2\pi}\exp\left(-\frac{R^{2}}{2}\right).\end{array}$
(191)

Thus,

$\begin{array}[]{c}R\psi_{R}\left(R\right)=-\frac{1}{2\pi}\exp\left(-\frac{R^{2}}{2}\right)+C,\end{array}$
(192)

where $C$ is an arbitrary constant. Next,

$\begin{array}[]{c}\psi\left(R\right)=\frac{1}{4\pi}E_{1}\left(\frac{R^{2}}{2}\right)+\frac{1}{2\pi}\ln\left(R\right),\end{array}$
(193)

where our choice of $C$ guarantees that $\psi$ has the right behavior when
$R\rightarrow 0$ and $R\rightarrow\infty$. Here $E_{1}\left(\eta\right)$ is
the exponential integral of the form:

$\begin{array}[]{c}E_{1}\left(\eta\right)=\int\limits_{\eta}^{\infty}\frac{e^{-\eta^{\prime}}}{\eta^{\prime}}d\eta^{\prime}.\end{array}$
(194)

##  6 Non-Gaussian Processes

###  6.1 The General Case

For financial applications, it is useful to consider a more general equation
governing SDE of the form:

$\begin{array}[]{c}d\mathbf{x}_{t}=\left(\mathbf{b}^{\left(x\right)}\left(t\right)+\mathfrak{B}^{\left(xx\right)}\left(t\right)\mathbf{x}_{t}+\mathfrak{B}^{\left(xy\right)}\left(t\right)\mathbf{y}_{t}\right)dt,\\\
\\\
d\mathbf{y}_{t}=\left(\mathbf{b}^{\left(y\right)}\left(t\right)+\mathfrak{B}^{\left(yx\right)}\left(t\right)\mathbf{x}_{t}+\mathfrak{B}^{\left(yy\right)}\left(t\right)\mathbf{y}_{t}\right)dt\\\
\\\
+\mathfrak{\mathbf{\Sigma}}^{\left(yy\right)}\left(t\right)\left(\mathrm{diag}\left(\mathbf{d}^{\left(0\right)}\left(t\right)+\mathfrak{D}\left(t\right)\mathbf{z}_{t}\right)\right)^{1/2}d\mathbf{W}_{t}^{\left(y\right)},\\\
\\\ \mathbf{x}_{\tau}=\mathbf{\xi,\ \ \ y}_{\tau}=\mathbf{\theta,}\end{array}$
(195)

or, more compactly,

$\begin{array}[]{c}d\mathbf{z}_{t}=\left(\mathbf{b}^{\left(z\right)}\left(t\right)+\mathfrak{B}^{\left(zz\right)}\left(t\right)\mathbf{z}_{t}\right)dt\\\
\\\ +\left(\begin{array}[]{c}0\\\
\mathbf{\Sigma}^{\left(yy\right)}\left(t\right)\left(\mathrm{diag}\left(\mathbf{d}^{\left(0\right)}\left(t\right)+\mathfrak{D}\left(t\right)\mathbf{z}_{t}\right)\right)^{1/2}d\mathbf{W}_{t}^{\left(y\right)}\end{array}\right),\\\
\\\ \mathbf{z}_{\tau}=\mathbf{\zeta=}\left(\begin{array}[]{c}\mathbf{\xi}\\\
\mathbf{\theta}\end{array}\right),\end{array}$ (196)

Here, in addition to the quantities introduced in the previous section, we
define an $\left(N\times 1\right)$ column vector
$\mathbf{d}^{\left(0\right)}$, and an $\left(N\times I\right)$ matrix
$\mathfrak{D}$. It is convenient to introduce vectors
$\mathbf{d}^{\left(i\right)}$ equal to the $i$-th column of $\mathfrak{D}$.

Since the corresponding $\left(N\times N\right)$ covariance matrix
$\mathfrak{A}$ has the form:

$\begin{array}[]{c}\mathfrak{A=}\mathbf{\Sigma}\left(\mathrm{diag}\left(\mathbf{d}^{\left(0\right)}\left(t\right)+\mathfrak{D}\left(t\right)\mathbf{z}_{t}\right)\right)^{1/2}\left(\mathbf{\Sigma}\left(\mathrm{diag}\left(\mathbf{d}^{\left(0\right)}\left(t\right)+\mathfrak{D}\left(t\right)\mathbf{z}_{t}\right)\right)^{1/2}\right)^{T},\end{array}$
(197)

it linearly depends on $\mathbf{z}$:

$\begin{array}[]{c}\mathfrak{A=}\mathbf{\Sigma}\left(\mathrm{diag}\left(\mathbf{d}^{\left(0\right)}+\mathfrak{D}\mathbf{z}\right)\right)\mathbf{\Sigma}^{T}=\mathfrak{A}^{\left(0\right)}+\mathop{\displaystyle\sum}\limits_{i=1}^{M+N}\mathfrak{A}^{\left(i\right)}z_{i},\end{array}$
(198)

where

$\begin{array}[]{c}\mathfrak{A}^{\left(0\right)}=\mathbf{\Sigma}\mathrm{diag}\left(\mathbf{d}^{\left(0\right)}\right)\Sigma^{T},\
\ \
\mathfrak{A}^{\left(i\right)}=\mathbf{\Sigma}\mathrm{diag}\left(\mathbf{d}^{\left(i\right)}\right)\Sigma^{T}.\end{array}$
(199)

In contrast to the Gaussian case, Eqs (196) have to be defined in the domain
$\boldsymbol{D}$ such that

$\boldsymbol{D=}\left\\{\left.\mathbf{z}\right|\mathbf{d}^{\left(0\right)}+\mathfrak{D}\mathbf{z\geq
0}\right\\},$ (200)

rather than in the whole space.

In financial engineering, covariance matrices of the form (197) were
introduced by Dai & Singleton (2000), they are discussed by Duffie et al.
(2003); Filipovic (2009) and many others.

The corresponding Fokker-Plank problem has the form:

$\begin{array}[]{c}P_{t}\left(\tau,\mathbf{\xi},\mathbf{\theta,}t,\mathbf{x},\mathbf{y}\right)-\frac{1}{2}\mathop{\displaystyle\sum}\mathop{\displaystyle\sum}\mathfrak{A}P_{\mathbf{zz}}\left(\tau,\mathbf{\xi},\mathbf{\theta,}t,\mathbf{x},\mathbf{y}\right)\\\
\\\ +\left(\mathbf{b}+\mathfrak{B}\mathbf{z}_{t}\right)\cdot
P_{\mathbf{z}}\left(\tau,\mathbf{\xi},\mathbf{\theta,}t,\mathbf{x},\mathbf{y}\right)+b^{\left(z\right)}P\left(\tau,\mathbf{\xi},\mathbf{\theta,}t,\mathbf{x},\mathbf{y}\right)=0,\\\
\\\
P\left(\tau,\mathbf{z},\tau,\mathbf{\zeta}\right)=\delta\left(\mathbf{z}-\mathbf{\zeta}\right),\end{array}$
(201)

Integral (86) expressing $P$ in terms of $H$ of the form:

$\begin{array}[]{c}H\left(\tau,\mathbf{\zeta,}t,\mathbf{z,m}\right)=\exp\left(\alpha\left(t\right)+i\Upsilon\left(t\right)\cdot\mathbf{z}-i\mathbf{m}\cdot\mathbf{\zeta}\right),\end{array}$
(202)

holds. The equations for $\alpha,\mathbf{\Upsilon}$ have the form:

$\begin{array}[]{c}\alpha^{\prime}\left(t\right)+i\mathbf{\Upsilon}^{\prime}\left(t\right)\cdot\mathbf{z}+\frac{1}{2}\mathbf{\Delta}\left(t\right)\cdot\mathfrak{A}\mathbf{\Delta}\left(t\right)+i\mathbf{\Upsilon}\left(t\right)\cdot\left(\mathbf{b}^{\left(z\right)}+\mathfrak{B}^{\left(zz\right)}\mathbf{z}\right)+b^{\left(z\right)}=0,\end{array}$
(203)

or, more explicitly,

$\begin{array}[]{c}\alpha^{\prime}\left(t\right)+i\mathbf{\Upsilon}^{\prime}\left(t\right)\cdot\mathbf{z}+\frac{1}{2}\mathbf{\Delta}\left(t\right)\cdot\mathfrak{A}^{\left(0\right)}\mathbf{\Delta}\left(t\right)+\sum\limits_{i=1}^{I}\frac{1}{2}\mathbf{\Delta}\left(t\right)\cdot\mathfrak{A}^{\left(i\right)}\mathbf{\Delta}\left(t\right)z_{i}\\\
\\\
+i\mathbf{\Upsilon}\left(t\right)\cdot\left(\mathbf{b}^{\left(z\right)}+\mathfrak{B}^{\left(zz\right)}\mathbf{z}\right)+b^{\left(z\right)}=0.\end{array}$
(204)

Thus, the system of ODEs for $\alpha,\mathbf{\Upsilon}$ can be written as
follows:

$\begin{array}[]{c}\alpha^{\prime}\left(t\right)+\frac{1}{2}\mathbf{\Delta}\left(t\right)\cdot\mathfrak{A}^{\left(0\right)}\left(t\right)\mathbf{\Delta}\left(t\right)+i\mathbf{\Upsilon}\left(t\right)\cdot\mathbf{b}^{\left(z\right)}\left(t\right)+b^{\left(z\right)}\left(t\right)=0,\
\ \ \alpha\left(\tau\right)=0,\\\ \\\
i\mathbf{\Upsilon}_{i}^{\prime}\left(t\right)+\frac{1}{2}\mathbf{\Delta}\left(t\right)\cdot\mathfrak{A}^{\left(i\right)}\mathbf{\Delta}\left(t\right)+\sum\limits_{i=1}^{I}i\mathfrak{B}_{ij}^{\left(zz\right)}\mathbf{\Upsilon}_{j}\left(t\right)=0,\
\ \ \mathbf{\Upsilon}_{i}\left(\tau\right)=m_{i}.\end{array}$ (205)

In the case in question, the equation for $\mathbf{\Upsilon}$ is no longer
linear. Instead, $\mathbf{\Upsilon}$ satisfies the so-called matrix Riccati
equation. Such equations are important for several applications, such as
optimal control. Solving a matrix Riccati equation is quite hard, so it is
more an art than a science. We shall report some of the results in this
direction below. However, in the one-dimensional case, the corresponding
Riccati equation can be converted into the second-order ODE, which can be
solved explicitly when the coefficients $\mathfrak{A}$, $\mathbf{b}$, $b$ are
time-independent.

###  6.2 Example: Feller Process

####  6.2.1 Feller Process

For benchmarking purposes, we start with deriving the with the well-known
t.p.d.f. for the Feller process; see Feller (1951, 1952):

$\begin{array}[]{c}dy=\left(\chi-\kappa y\right)dt+\varepsilon\sqrt{y}dZ,\ \ \
y\left(\tau\right)=\theta.\end{array}$ (206)

To guarantee that the process does not hit zero, we assume that

$\begin{array}[]{c}\vartheta=\frac{2\chi}{\varepsilon^{2}}-1>0.\end{array}$
(207)

The corresponding Fokker-Planck problem has the form:

$\begin{array}[]{c}P_{t}-\frac{1}{2}\varepsilon^{2}yP_{yy}+\left(\chi-\varepsilon^{2}-\kappa
y\right)P_{y}-\kappa P=0,\\\ \\\
P\left(\tau,\theta,t,y\right)=\delta\left(y-\theta\right).\end{array}$ (208)

The associated function $H\left(\tau,\theta,t,y,l\right)$ has the form:

$\begin{array}[]{c}H=\exp\left(\alpha\left(t\right)+i\Delta\left(t\right)y-il\theta\right),\end{array}$
(209)

so that

$\begin{array}[]{c}\alpha^{\prime}+i\left(\chi-\varepsilon^{2}\right)\Delta-\kappa=0,\
\ \ \alpha\left(\tau\right)=0,\\\ \\\ \ \
i\Delta^{\prime}+\frac{1}{2}\varepsilon^{2}\Delta^{2}-i\kappa\Delta=0,\ \ \
\Delta\left(\tau\right)=l.\end{array}$ (210)

Thus, $\Delta$ solves a nonlinear Riccati equation, which can be linearized
via a familiar substitution

$\begin{array}[]{c}\Delta=\frac{2i\Omega^{\prime}}{\varepsilon^{2}\Omega}.\end{array}$
(211)

As a result, we get:

$\begin{array}[]{c}\Omega^{\prime\prime}-\kappa\Omega^{\prime}=0,\ \ \
\Omega\left(\tau\right)=1,\ \ \
\Omega^{\prime}\left(\tau\right)=-\frac{i\varepsilon^{2}l}{2},\end{array}$
(212)
$\begin{array}[]{c}\alpha^{\prime}-2\left(\frac{\chi}{\varepsilon^{2}}-1\right)\left(\ln\left(\Omega\right)\right)^{\prime}-\kappa=0,\
\ \ \alpha\left(\tau\right)=0.\end{array}$ (213)

Accordingly,

$\begin{array}[]{c}\Omega=\frac{2\kappa-i\varepsilon^{2}\left(e^{\kappa\left(t-\tau\right)}-1\right)l}{2\kappa},\end{array}$
(214)
$\begin{array}[]{c}\Omega^{\prime}=-\frac{i\varepsilon^{2}e^{\kappa\left(t-\tau\right)}l}{2},\end{array}$
(215)
$\begin{array}[]{c}\Delta=\frac{2i\Omega^{\prime}}{\varepsilon^{2}\Omega}=\frac{2\kappa
le^{-\kappa\left(t-\tau\right)}}{2\kappa
e^{-\kappa\left(t-\tau\right)}-i\varepsilon^{2}\left(1-e^{-\kappa\left(t-\tau\right)}\right)l},\end{array}$
(216)
$\begin{array}[]{c}\alpha=\kappa\left(t-\tau\right)+2\left(\frac{\chi}{\varepsilon^{2}}-1\right)\ln\left(\frac{2\kappa-i\varepsilon^{2}l\left(e^{\kappa\left(t-\tau\right)}-1\right)}{2\kappa}\right),\end{array}$
(217)

and

$\begin{array}[]{c}H=\exp\left(\kappa\left(t-\tau\right)+2\left(\frac{\chi}{\varepsilon^{2}}-1\right)\ln\left(\frac{2\kappa-i\varepsilon^{2}l\left(e^{\kappa\left(t-\tau\right)}-1\right)}{2\kappa}\right)\right.\\\
\\\ \left.+il\left(\frac{2\kappa}{2\kappa
e^{-\kappa\left(t-\tau\right)}-i\varepsilon^{2}l\left(1-e^{-\kappa\left(t-\tau\right)}\right)}y-\theta\right)\right).\end{array}$
(218)

New variable, $l\rightarrow\eta$:

$\begin{array}[]{c}\eta=\frac{\varepsilon^{2}\left(1-e^{-\kappa\left(\tau-\tau\right)}\right)}{4\kappa
e^{-\kappa\left(t-\tau\right)}}l,\ \ \ l=\frac{4\kappa
e^{-\kappa\left(t-\tau\right)}}{\varepsilon^{2}\left(1-e^{-\kappa\left(\tau-\tau\right)}\right)}\eta.\end{array}$
(219)

and rescaled $H$, $Hdl\rightarrow\tilde{H}d\eta$:

$\begin{array}[]{c}\tilde{H}=2M\exp\left(2\left(\frac{\chi}{\varepsilon^{2}}-1\right)\ln\left(1-2i\eta\right)+2Mi\eta\left(\frac{y}{1-2i\eta}-\bar{\theta}\right)\right),\end{array}$
(220)

where $M$ appears due to the change of variables, and $N$ is a prefactor:

$\begin{array}[]{c}M=\frac{2\kappa}{\varepsilon^{2}\left(1-e^{-\kappa\left(\tau-\tau\right)}\right)},\\\
\\\ N=\frac{4\kappa
e^{-\kappa\left(t-\tau\right)}}{\varepsilon^{2}\left(1-e^{-\kappa\left(\tau-\tau\right)}\right)}=2MZ,\\\
\\\ Z=e^{-\kappa\left(t-\tau\right)},\\\ \\\
\bar{\theta}=e^{-\kappa\left(t-\tau\right)}\theta=Z\theta.\end{array}$ (221)

Finally,

$\begin{array}[]{c}Q\left(\tau,\theta,t,y\right)=\frac{2M}{2\pi}\int_{-\infty}^{\infty}\exp\left(2\left(\frac{\chi}{\varepsilon^{2}}-1\right)\ln\left(1-2i\eta\right)+2Mi\eta\left(\frac{y}{1-2i\eta}-\bar{\theta}\right)\right)d\eta,\end{array}$
(222)

A well-known formula yields

$\begin{array}[]{c}Q\left(\tau,\theta,t,y\right)=Me^{-M\left(y+\bar{\theta}\right)}\left(\frac{\bar{\theta}}{y}\right)^{-\frac{\vartheta}{2}}I_{-\vartheta}\left(2M\sqrt{y\bar{\theta}}\right)\\\
\\\
=Me^{-M\left(y+\bar{\theta}\right)}\left(\frac{y}{\bar{\theta}}\right)^{\frac{\vartheta}{2}}I_{\vartheta}\left(2M\sqrt{y\bar{\theta}}\right),\end{array}$
(223)

see, e.g., Lipton (2001) and references therein. The density $Q\left(y\right)$
integrates to one:

$\begin{array}[]{c}\int_{0}^{\infty}Q\left(\tau,\theta,t,y\right)dy=\int_{0}^{\infty}e^{-u-v}\left(\frac{v}{u}\right)^{\vartheta/2}I_{\vartheta}\left(2\sqrt{uv}\right)dv=1,\end{array}$
(224)

where $u=M\bar{\theta},\ \ \ v=My$. We shall use this fact below.

This derivation can be extended to cover augmented Feller processes, as we
shall do next. It is worth noting that studying the Feller process with time-
dependent coefficients is rather tricky and remains to be wholly done; see,
however, Masoliver (2016); Giorno & Nobile (2021) and references therein.

####  6.2.2 Degenerate Augmented Feller Process

Equations

$\begin{array}[]{c}dx=ydt,\ \ \ x_{\tau}=\xi,\\\ \\\ dy=\left(\chi-\kappa
y\right)dt+\varepsilon\sqrt{y}dZ,\ \ \ y_{\tau}=\theta.\end{array}$ (225)

The Fokker-Planck equation:

$\begin{array}[]{c}0=P_{t}-\frac{1}{2}\varepsilon^{2}\left(yP\right)_{yy}+yP_{x}+\left(\left(\chi-\kappa
y\right)P\right)_{y}\\\ \\\
=P_{t}-\frac{1}{2}\varepsilon^{2}yP_{yy}-\varepsilon^{2}P_{y}+yP_{x}+\left(\chi-\kappa
y\right)P_{y}-\kappa P\\\ \\\
=P_{t}-\frac{1}{2}\varepsilon^{2}yP_{yy}+yP_{x}+\left(\chi-\varepsilon^{2}-\kappa
y\right)P_{y}-\kappa P.\end{array}$ (226)

The corresponding function $H$ has the form:

$\begin{array}[]{c}H=\exp\left(\alpha\left(t\right)+ik\left(x-\xi\right)+i\Delta\left(t\right)y-il\theta\right),\end{array}$
(227)

where

$\begin{array}[]{c}\alpha^{\prime}+i\Delta^{\prime}y+\frac{1}{2}\varepsilon^{2}\Delta^{2}y+iky+i\Delta\left(\chi-\varepsilon^{2}-\kappa
y\right)-\kappa=0,\\\ \\\ \alpha\left(\tau\right)=0,\ \ \
\Delta\left(\tau\right)=l.\end{array}$ (228)

Thus

$\begin{array}[]{c}\alpha^{\prime}+i\left(\chi-\varepsilon^{2}\right)\Delta-\kappa=0,\
\ \ \alpha\left(\tau\right)=0,\\\ \\\
i\Delta^{\prime}+\frac{1}{2}\varepsilon^{2}\Delta^{2}-i\kappa\Delta+ik=0,\ \ \
\Delta\left(\tau\right)=l,\end{array}$ (229)

As before, we can linearize the Riccati equation for $\Delta$ via a two-step
linearization

$\begin{array}[]{c}\Delta=\frac{2i\Lambda}{\varepsilon^{2}},\ \ \
\Lambda=\frac{\Omega^{\prime}}{\Omega},\end{array}$ (230)

which yields

$\begin{array}[]{c}\alpha^{\prime}-2\left(\frac{\chi}{\varepsilon^{2}}-1\right)\left(\ln\left(\Omega\right)\right)^{\prime}-\kappa=0,\
\ \ \alpha\left(\tau\right)=0,\\\ \\\
\Omega^{\prime\prime}-\kappa\Omega^{\prime}-\frac{i\varepsilon^{2}k}{2}\Omega=0,\
\ \ \Omega\left(\tau\right)=1,\ \ \
\Omega^{\prime}\left(\tau\right)=-\frac{i\varepsilon^{2}l}{2}.\end{array}$
(231)

We can represent $\Omega$ in the form:

$\begin{array}[]{c}\Omega=c_{+}e^{\lambda_{+}\left(t-\tau\right)}+c_{-}e^{\lambda_{+}\left(t-\tau\right)},\end{array}$
(232)

where $\lambda_{\pm}$ are solutions of the characteristic equation:

$\begin{array}[]{c}\lambda^{2}-\kappa\lambda-\frac{i\varepsilon^{2}k}{2}=0,\\\
\\\ \lambda_{\pm}=\mu\pm\zeta,\\\ \\\ \mu=\frac{\kappa}{2},\ \ \
\zeta=\frac{\sqrt{\kappa^{2}+2i\varepsilon^{2}k}}{2},\end{array}$ (233)

and $c_{\pm}$ satisfy the following system of linear equations:

$\begin{array}[]{c}c_{+}+c_{-}=1,\\\ \\\
\lambda_{+}c_{+}+\lambda_{-}c_{-}=-\frac{i\varepsilon^{2}l}{2}.\end{array}$
(234)

so that

$\begin{array}[]{c}c_{+}=-\frac{\left(2\lambda_{-}+i\varepsilon^{2}l\right)}{4\zeta},\
\ \
c_{-}=\frac{\left(2\lambda_{+}+i\varepsilon^{2}l\right)}{4\zeta}.\end{array}$
(235)

It is useful to note that

$\begin{array}[]{c}\lambda_{+}\lambda_{-}=\mu^{2}-\zeta^{2}=-\frac{i\varepsilon^{2}k}{2}.\end{array}$
(236)

For the sake of clarity, we introduce the following notation:

$\begin{array}[]{c}\mathsf{E}_{0}=e^{\mu\left(t-\tau\right)},\ \ \
\mathsf{E}_{\pm}=e^{\pm\zeta\left(t-\tau\right)}.\end{array}$ (237)

Thus,

$\begin{array}[]{c}\Omega=\mathsf{E}_{0}\left(\frac{-\left(2\lambda_{-}+i\varepsilon^{2}l\right)\mathsf{E}_{+}+\left(2\lambda_{+}+i\varepsilon^{2}l\right)\mathsf{E}_{-}}{4\zeta}\right),\end{array}$
(238)
$\begin{array}[]{c}\Omega^{\prime}=\mathsf{E}_{0}\left(\frac{-\lambda_{+}\left(2\lambda_{-}+i\varepsilon^{2}l\right)\mathsf{E}_{+}+\lambda_{-}\left(2\lambda_{+}+i\varepsilon^{2}l\right)\mathsf{E}_{-}}{4\zeta}\right),\end{array}$
(239)
$\begin{array}[]{c}\Delta=\frac{2i\left(-\lambda_{+}\left(2\lambda_{-}+i\varepsilon^{2}l\right)\mathsf{E}_{+}+\lambda_{-}\left(2\lambda_{+}+i\varepsilon^{2}l\right)\mathsf{E}_{-}\right)}{\varepsilon^{2}\left(-\left(2\lambda_{-}+i\varepsilon^{2}l\right)\mathsf{E}_{+}+\left(2\lambda_{+}+i\varepsilon^{2}l\right)\mathsf{E}_{-}\right)},\end{array}$
(240)
$\begin{array}[]{c}\alpha=2\left(\frac{\chi}{\varepsilon^{2}}-1\right)\ln\left(\frac{-\left(2\lambda_{-}+i\varepsilon^{2}l\right)\mathsf{E}_{+}+\left(2\lambda_{+}+i\varepsilon^{2}l\right)\mathsf{E}_{-}}{4\zeta}\right)+\frac{2\chi\mu\left(t-\tau\right)}{\varepsilon^{2}}.\end{array}$
(241)

Accordingly, $H$ can be written in the form:

$\begin{array}[]{c}H=\exp\left(\frac{2\chi\mu\left(t-\tau\right)}{\varepsilon^{2}}+ik\left(x-\xi\right)\right.\\\
\\\
+2\left(\frac{\chi}{\varepsilon^{2}}-1\right)\ln\left(\frac{2\left(-\lambda_{-}\mathsf{E}_{+}+\lambda_{+}\mathsf{E}_{-}\right)-i\varepsilon^{2}l\left(\mathsf{E}_{+}-\mathsf{E}_{-}\right)}{4\zeta}\right)\\\
\\\
\left.+\frac{2\left(2\lambda_{+}\lambda_{-}\left(\mathsf{E}_{+}-\mathsf{E}_{-}\right)+i\varepsilon^{2}l\left(\lambda_{+}\mathsf{E}_{+}-\lambda_{-}\mathsf{E}_{-}\right)\right)y}{\varepsilon^{2}\left(2\left(-\lambda_{-}\mathsf{E}_{+}+\lambda_{+}\mathsf{E}_{-}\right)-i\varepsilon^{2}l\left(\mathsf{E}_{+}-\mathsf{E}_{-}\right)\right)}-il\theta\right).\end{array}$
(242)

Change of variables:

$\begin{array}[]{c}\eta=\frac{\varepsilon^{2}\left(\mathsf{E}_{+}-\mathsf{E}_{-}\right)l}{4\left(-\lambda_{-}\mathsf{E}_{+}+\lambda_{+}\mathsf{E}_{-}\right)},\\\
\\\
l=\frac{4\left(-\lambda_{-}\mathsf{E}_{+}+\lambda_{+}\mathsf{E}_{-}\right)\eta}{\varepsilon^{2}\left(\mathsf{E}_{+}-\mathsf{E}_{-}\right)}\equiv
N\eta.\end{array}$ (243)

allows one to factorize $H$ as follows:

$\begin{array}[]{c}\tilde{H}=\tilde{H}_{1}\tilde{H}_{2},\end{array}$ (244)

where

$\begin{array}[]{c}\tilde{H}_{1}=\exp\left(\frac{2\chi\mu\left(t-\tau\right)}{\varepsilon^{2}}+ik\left(x-\xi\right)+\frac{2\chi}{\varepsilon^{2}}\ln\left(\frac{\left(-\lambda_{-}\mathsf{E}_{+}+\lambda_{+}\mathsf{E}_{-}\right)}{2\zeta}\right)\right.\\\
\\\
\left.+\frac{2\lambda_{+}\lambda_{-}\left(\mathsf{E}_{+}-\mathsf{E}_{-}\right)y}{\varepsilon^{2}\left(-\lambda_{-}\mathsf{E}_{+}+\lambda_{+}\mathsf{E}_{-}\right)}\right),\\\
\\\
\tilde{H}_{2}=2M\exp\left(2\left(\frac{\chi}{\varepsilon^{2}}-1\right)\ln\left(1-2i\eta\right)+2Mi\eta\left(\frac{y}{1-2i\eta}-\bar{\theta}\right)\right).\end{array}$
(245)

and

$\begin{array}[]{c}M=\frac{8\zeta^{2}}{\varepsilon^{2}\left(\mathsf{E}_{+}-\mathsf{E}_{-}\right)\left(-\lambda_{-}\mathsf{E}_{+}+\lambda_{+}\mathsf{E}_{-}\right)},\\\
\\\
Z=\frac{\left(-\lambda_{-}\mathsf{E}_{+}+\lambda_{+}\mathsf{E}_{-}\right)^{2}}{4\zeta^{2}}=\frac{N}{2M},\\\
\\\ \bar{\theta}=Z\theta.\end{array}$ (246)

Integration with respect to $\eta$ can be done analytically

$\begin{array}[]{c}Q_{2}\left(\tau,\theta,t,y\right)=\frac{2M}{2\pi}\int_{-\infty}^{\infty}\tilde{H}_{2}d\eta=Me^{-M\left(y+\bar{\theta}\right)}\left(\frac{y}{\bar{\theta}}\right)^{\frac{\vartheta}{2}}I_{\vartheta}\left(2M\sqrt{y\bar{\theta}}\right).\end{array}$
(247)

which allows us to calculate $P$ via a single inverse Fourier transform:

$\begin{array}[]{c}P=\frac{1}{2\pi}\int_{-\infty}^{\infty}\exp\left(\frac{2\chi\mu\left(t-\tau\right)}{\varepsilon^{2}}+ik\left(x-\xi\right)\right.\\\
\\\
\left.+\frac{2\chi}{\varepsilon^{2}}\ln\left(\frac{\left(-\lambda_{-}\mathsf{E}_{+}+\lambda_{+}\mathsf{E}_{-}\right)}{2\zeta}\right)-Ry\right)\\\
\\\ \times
Me^{-M\left(y+\bar{\theta}\right)}\left(\frac{y}{\bar{\theta}}\right)^{\frac{\vartheta}{2}}I_{\vartheta}\left(2M\sqrt{y\bar{\theta}}\right)dk\\\
\\\
=\frac{1}{2\pi}\int_{-\infty}^{\infty}\exp\left(\frac{2\chi\mu\left(t-\tau\right)}{\varepsilon^{2}}+ik\left(x-\xi\right)\right.\\\
\\\
\left.+\frac{2\chi}{\varepsilon^{2}}\ln\left(\frac{\left(-\lambda_{-}\mathsf{E}_{+}+\lambda_{+}\mathsf{E}_{-}\right)}{2\zeta}\right)\right)\\\
\\\ \times
Me^{-\left(M+R\right)y-M\bar{\theta}}\left(\frac{y}{\bar{\theta}}\right)^{\frac{\vartheta}{2}}I_{\vartheta}\left(2M\sqrt{y\bar{\theta}}\right)dk\end{array}$
(248)

where

$\begin{array}[]{c}R=-\frac{2\lambda_{+}\lambda_{-}\left(\mathsf{E}_{+}-\mathsf{E}_{-}\right)}{\varepsilon^{2}\left(-\lambda_{-}\mathsf{E}_{+}+\lambda_{+}\mathsf{E}_{-}\right)}.\end{array}$
(249)

This expression is not as convenient to deal with due to the presence of the
term $\exp\left(-Ry\right)$. We wish to rewrite it in a form better suited for
further calculations. Further transformations to use the density and transform
extra dependence on $y$ to dependence on $\theta$. We proceed as follows.
Define

$\begin{array}[]{c}Z_{1}=\left(\frac{M}{M+R}\right)^{2}Z.\end{array}$ (250)

Then

$\begin{array}[]{c}Me^{-\left(M+R\right)y-MZ\theta}\left(\frac{y}{Z\theta}\right)^{\frac{\vartheta}{2}}I_{\vartheta}\left(2M\sqrt{yZ\theta}\right)\\\
\\\
=\frac{M}{\left(M+R\right)}\left(M+R\right)e^{-\left(M+R\right)y-\frac{\left(M+R\right)^{2}}{M}Z_{1}\theta}\left(\frac{M}{M+R}\right)^{\vartheta}\\\
\\\
\times\left(\frac{y}{Z_{1}\theta}\right)^{\frac{\vartheta}{2}}I_{\vartheta}\left(2\left(M+R\right)\sqrt{yZ_{1}\theta}\right)\\\
\\\
=\left(\frac{M}{M+R}\right)^{\vartheta+1}e^{-\left(\frac{\left(M+R\right)^{2}}{M}-\left(M+R\right)\right)Z_{1}\theta}\left(M+R\right)e^{-\left(M+R\right)\left(y+Z_{1}\theta\right)}\\\
\\\
\times\left(\frac{y}{Z_{1}\theta}\right)^{\frac{\vartheta}{2}}I_{\vartheta}\left(2\left(M+R\right)\sqrt{yZ_{1}\theta}\right)\\\
\\\
=\left(\frac{M}{M+R}\right)^{\vartheta+1}e^{-\frac{MRZ}{\left(M+R\right)}\theta}\left(M+R\right)e^{-\left(M+R\right)\left(y+Z_{1}\theta\right)}\\\
\\\
\times\left(\frac{y}{Z_{1}\theta}\right)^{\frac{\vartheta}{2}}I_{\vartheta}\left(2\left(M+R\right)\sqrt{yZ_{1}\theta}\right),\end{array}$
(251)

so that

$\begin{array}[]{c}P\left(t,x,y\right)=\frac{1}{2\pi}\int_{-\infty}^{\infty}\exp\left(\frac{2\chi\mu\left(t-\tau\right)}{\varepsilon^{2}}+ik\left(x-\xi\right)\right.\\\
\\\
\left.+\frac{2\chi}{\varepsilon^{2}}\ln\left(\frac{M\left(-\lambda_{-}\mathsf{E}_{+}+\lambda_{+}\mathsf{E}_{-}\right)}{2\zeta\left(M+R\right)}\right)-\frac{MRZ\theta}{M+R}\right)\\\
\\\
\times\left(M+R\right)e^{-\left(M+R\right)\left(y+Z_{1}\theta\right)}\left(\frac{y}{Z_{1}\theta}\right)^{\frac{\vartheta}{2}}I_{\vartheta}\left(2\left(M+R\right)\sqrt{yZ_{1}\theta}\right)dk.\end{array}$
(252)

Explicitly,

$\begin{array}[]{c}M+R=\frac{2\left(\left(\lambda_{+}^{2}+\lambda_{-}^{2}\right)-\lambda_{+}\lambda_{-}\left(\mathsf{E}_{+}^{2}+\mathsf{E}_{-}^{2}\right)\right)}{\varepsilon^{2}\left(\mathsf{E}_{+}-\mathsf{E}_{-}\right)\left(-\lambda_{-}\mathsf{E}_{+}+\lambda_{+}\mathsf{E}_{-}\right)},\\\
\\\
\frac{M}{M+R}=\frac{\left(\lambda_{+}-\lambda_{-}\right)^{2}}{\left(\lambda_{+}^{2}+\lambda_{-}^{2}\right)-\lambda_{+}\lambda_{-}\left(\mathsf{E}_{+}^{2}+\mathsf{E}_{-}^{2}\right)},\\\
\\\
\frac{MRZ}{M+R}=-\frac{2\lambda_{+}\lambda_{-}\left(\mathsf{E}_{+}-\mathsf{E}_{-}\right)}{\varepsilon^{2}\left(\lambda_{+}\mathsf{E}_{+}-\lambda_{-}\mathsf{E}_{-}\right)},\end{array}$
(253)

so that

$\begin{array}[]{c}P\left(t,x,y\right)=\frac{1}{2\pi}\int_{-\infty}^{\infty}\exp\left(\frac{2\chi\mu\left(t-\tau\right)}{\varepsilon^{2}}+ik\left(x-\xi\right)\right.\\\
\\\
\left.+\frac{2\chi}{\varepsilon^{2}}\ln\left(\frac{2\zeta}{\lambda_{+}\mathsf{E}_{+}-\lambda_{-}\mathsf{E}_{-}}\right)+\frac{2\lambda_{+}\lambda_{-}\left(\mathsf{E}_{+}-\mathsf{E}_{-}\right)\theta}{\varepsilon^{2}\left(\lambda_{+}\mathsf{E}_{+}-\lambda_{-}\mathsf{E}_{-}\right)}\right)\\\
\\\