Skip to main content Accessibility help
×
Hostname: page-component-586b7cd67f-g8jcs Total loading time: 0 Render date: 2024-11-30T20:53:47.934Z Has data issue: false hasContentIssue false

Hydrodynamics of Markets

Hidden Links between Physics and Finance

Published online by Cambridge University Press:  26 November 2024

Alexander Lipton
Affiliation:
ADIA Lab

Summary

An intriguing link between a wide range of problems occurring in physics and financial engineering is presented. These problems include the evolution of small perturbations of linear flows in hydrodynamics, the movements of particles in random fields described by the Kolmogorov and Klein-Kramers equations, the Ornstein-Uhlenbeck and Feller processes, and their generalizations. They are reduced to affine differential and pseudo-differential equations and solved in a unified way by using Kelvin waves and developing a comprehensive math framework for calculating transition probabilities and expectations. Kelvin waves are instrumental for studying the well-known Black-Scholes, Heston, and Stein-Stein models and more complex path-dependent volatility models, as well as the pricing of Asian options, volatility and variance swaps, bonds, and bond options. Kelvin waves help to solve several cutting-edge problems, including hedging the impermanent loss of Automated Market Makers for cryptocurrency trading. This title is also available as Open Access on Cambridge Core.
Type
Element
Information
Online ISBN: 9781009503129
Publisher: Cambridge University Press
Print publication: 31 December 2024
Creative Commons
Creative Common License - CCCreative Common License - BYCreative Common License - NC
This content is Open Access and distributed under the terms of the Creative Commons Attribution licence CC-BY-NC 4.0 https://creativecommons.org/cclicenses/

6accdae13eff 7i3l9n4o4qrr4s8t12ux

Letter from Isaac Newton to Henry Oldenburg, 24 October 1676Footnote 1

Bob Montagnet:

Yeah, good choice, Vlad.
Get back to the security system. How does it work?

Vlad:

The way everything works. Mathematics.
“The Good Thief,” screenplay by Neil Jordan, 2002

1 Introduction

1.1 Background

Newton’s discovery of differential equations and calculus was crucial in developing classical mechanics because it allowed for the mathematical description of the motion of objects. This discovery took a groundbreaking step in unifying mathematics with physics, enabling the prediction of planetary orbits, the motion of objects under various forces, and much more, and marked the beginning of a new era in mathematics and science, laying the cornerstone for over three centuries of advancements.

Newton understood the immediate impact of his discoveries and their potential to transform the understanding of the natural world. To establish and protect his intellectual property rights at the same time, he concealed his discovery in the fundamental anagram of calculus, which he included in his 1676 letter to Oldenburg. This anagram contained a Latin statement describing the method of fluxions (his term for calculus) when decoded. The need for an anagram reflected that Newton was competitive and cautious in equal measure by balancing the desire for recognition with the fear of disclosure. The number of occurrences of each Latin character in Newton’s sentence agrees with his anagram, thus proving that the actual sentence was written in 1676.Footnote 2 The original letter is shown in Figure 1.

Figure 1 Newton’s letter to Oldenburg, 1676.

Reproduced by kind permission of the Syndics of Cambridge University Library.

The fact that differential equations are instrumental in mathematics and physics alike was firmly established in the late seventeenth century. However, methods for solving these equations remained ad hoc for more than a century until the work by Lagrange, Laplace, Fourier, and many other mathematicians and physicists. In particular, the Fourier transform stands out as the most potent tool in an applied mathematician’s toolkit, enabling the solving of linear partial differential equations (PDEs) and partial pseudo-differential equations (PPDEs) with spatially constant coefficients; it is also invaluable for analyzing time series and tackling other critical tasks (Reference FourierFourier (1822); Reference Morse and FeshbachMorse & Feschbach (1953)).

At the heart of the n-dimensional Fourier method are wave functions, expressed as follows:

(t,x,k)=a(t)exp(ikx),(1.1)

where x and k are n-dimensional vectors, denotes the scalar product, at is the amplitude, and kx is the phase. Depending on the particular problem at hand, the amplitude at can be a scalar or a vector, hence the notation. Substituting F into a PDE with spatially constant coefficients, one reduces the problem of interest to a system of ordinary differential equations (ODEs) or a single ODE when at is scalar. Of course, this system parametrically depends on k.

This Element studies PDEs and PPDEs with coefficients linearly dependent on x, which are called affine. Hence, one must use a more general approach and consider wave functions with time-dependent wave vectors: Reference KelvinKelvin (1887) and Reference OrrOrr (1907) were the first to use such waves to analyze the stability of the steady motions of an incompressible fluid.

Kt,x,βt=atexpiβtx.(1.2)

Affine problems are not artificial constructs. They appear organically in several situations, for example, when the linear description of the underlying physical mechanism is either exact or provides an excellent approximation to reality or when the evolution in the phase space is studied; see Section 3.

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. The Ornstein–Uhlenbeck (OU) and Feller processes are the simplest but extremely important examples of affine processes; see Reference Uhlenbeck and OrnsteinUhlenbeck and Ornstein (1930), Reference ChandrasekharChandresekhar (1943), and Reference FellerFeller (1951, Reference Feller1952). For financial applications of affine processes see Reference Duffie and KanDuffie and Kan (1996), Reference Duffie, Pan and SingletonDuffie et al. (2000), Reference Dai and SingletonDai and Singleton (2000), Reference LiptonLipton (2001), Reference Duffie, Filipovic and SchachermayerDuffie et al. (2003), Reference SeppSepp (2007), Reference Lipton and SeppLipton and Sepp (2008), and Reference FilipovicFilipovic (2009), among others.

This Element uses Kelvin waves of the form (1.2) to study transition probability density functions (t.p.d.fs) for affine stochastic processes. These processes can be either degenerate, namely, have more independent components than the sources of uncertainty, or nondegenerate, when every component has its source of uncertainty. Recall that the t.p.d.f. for a stochastic process describes the likelihood of a system transitioning from one state to another over a specified period. Knowing the iterated t.p.d.f. is fundamental for understanding the dynamics and behavior of stochastic processes over time and is tantamount to knowing the process itself.

In this Element, Kelvin waves are also used to solve several essential and intricate problems occurring in financial applications. These include pricing options with stochastic volatility, path-dependent options, and Asian options with geometric averaging, among many others.

The main objective is to link various financial engineering topics with their counterparts in hydrodynamics and molecular physics and showcase the interdisciplinary nature of quantitative finance and economic modeling. Finding such connections allows us to understand better how to model, price, and risk-manage various financial instruments, derive several new results, and provide additional intuition regarding their salient features. This Element continues previous efforts in this direction; see Reference Lipton and SeppLipton and Sepp (2008) and Reference LiptonLipton (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, for example, Reference OvsiannikovOvsiannikov (1982), Reference IbragimovIbragimov (1985), Reference OlverOlver (1986), and Reference Bluman and KumeiBluman and Kumei (1989), while their specific applications to affine equations are covered by Reference BerestBerest (1993), Reference AksenovAksenov (1995), Reference Craddock and PlatenCraddock and Platen (2004), Reference CraddockCraddock (2012), and Reference Kovalenko, Stogniy and TertychnyiKovalenko et al. (2014), among many others. However, Lie symmetry techniques are exceedingly cumbersome and might be challenging to use in practice, especially when complicated affine equations are considered.

Laplace transform of spatial variables can be used in some cases, for instance, for Feller processes; see, for example, Reference FellerFeller (1951, Reference Feller1952). However, they are hard to use for solving generic affine equations.

Reductions of a given equation to a simpler, solvable form is another powerful method that can be successfully used in many instances; see, for example, Reference ChandrasekharChandresekhar (1943), Reference Carr, Lipton and MadanCarr et al. (2002), Reference Lipton, Gal and LasisLipton et al. (2014), and Reference LiptonLipton (2018), chapter 9. Although the reduction method is quite powerful, experience suggests it is often hard to use in practice.

Finally, the affine ansatz based on Kelvin waves provides yet another approach, which is the focus of the present Element; see also Reference Duffie and KanDuffie and Kan (1996), Reference Dai and SingletonDai and Singleton (2000), Reference Duffie, Filipovic and SchachermayerDuffie et al. (2003), Reference Lipton and SeppLipton and Sepp (2008), Reference FilipovicFilipovic (2009), and Reference LiptonLipton (2018), chapter 12. Undoubtedly, the affine framework, also known as the affine ansatz, is the most potent among the abovementioned techniques due to its comprehensive nature, versatility, and (relative) ease of use, even in complex situations. In practice, applications of Kelvin waves consist of three steps:

  • Effectively separating variables for the evolution problems with pseudo-differential generators linearly dependent on spatial coordinates;

  • Solving ODEs parametrized by time-dependent wave vectors; see (1.2);

  • Aggregating their solutions together to get the solution to the original problem.

However, despite being a ruthlessly efficient tool, Kelvin waves have limitations – using them to solve evolution problems supplied with external boundary conditions is challenging. This exciting topic is being actively researched now; it will be discussed elsewhere in due course.

1.2 Main Results

This Element develops a coherent, unified mathematical framework using Kelvin waves as a powerful and versatile tool for studying t.p.d.fs in the context of generic affine processes. It discovers previously hidden connections among large classes of apparently unrelated problems from hydrodynamics, molecular physics, and financial engineering. All these problems require solving affine (pseudo-) differential equations, namely, equations with coefficients, which linearly depend on spatial variables. The Element discusses some classical results and derives several original ones related to:

  • small wave-like perturbations of linear flows of ideal and viscous fluids described by Euler and Navier–Stokes equations, respectively;

  • motions of free and harmonically bound particles under the impact of random external white-noise forces described by the Klein–Kramers equations and the hypoelliptic Kolmogorov equation, which play an essential role in statistical physics;

  • Gaussian and non-Gaussian affine processes, such as the Ornstein–Uhlenbeck and Feller processes, which are the archetypal mean-reverting processes, and their generalizations;

  • dynamics of financial markets, particularly derivative products.

To solve some of the more complicated problems, one must augment primary processes by introducing subordinate processes for auxiliary variables, such as integrals over the original stochastic variable, and develop a uniform mathematical formalism to construct t.p.d.fs for the abovementioned processes.

Quite unexpectedly, the analysis identifies and rectifies an error in the original solution of the Kolmogorov equation. The rectified solution is dimensionally correct, properly scales when the process parameters change, and agrees with numerical results.

Furthermore, this Element derives many original results and extends and reinterprets some well-known ones. For instance, it develops a concise and efficient expression for t.p.d.fs in the case of processes with stochastic volatility. Moreover, the 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 is used to find a new explicit expression for the vorticity of a two-dimensional flow in terms of the Gaussian density.

Finally, the Element applies the new methodology to various financial engineering topics, such as pricing options with stochastic volatility, options with path-dependent volatility, Asian options, volatility and variance swaps, options on stocks with path-dependent volatility, and bonds and bond options. In contrast to the classical approach, the Element treats primary fixed-income products, such as bonds and bond options, as path-dependent, allowing us to gain additional intuition regarding such products’ pricing and risk management. It also highlights the flexibility of the interdisciplinary framework by incorporating additional complexities into the picture, such as jump-diffusion processes and, more generally, processes driven by affine pseudo-differential processes frequently used in financial applications.

1.3 Element Structure

Section 2 introduces Kelvin waves. Section 2.1 introduces the Euler equations, which describe the dynamics of a perfect fluid, alongside the Navier–Stokes equation for viscous incompressible fluids. Section 2.2 discusses the exact equilibria of these equations, focusing on states where velocity varies linearly and pressure quadratically with spatial coordinates, referred to as linear flows. Section 2.3 illustrates that the renowned Kelvin waves provide solutions to the linearized Euler and Navier–Stokes equations for small perturbations of the linear flows. This section also explores the use of Kelvin waves in analyzing the stability of these flows.

The Element uses Kelvin waves as a fundamental tool in the analytical arsenal, demonstrating their applicability across various study areas. For instance, they allow one to discover profound and surprising links between the viscous two-dimensional vorticity equations and the Klein–Kramers equation, a cornerstone of stochastic physics; see Section 6.6. This connection results in a novel formula representing vorticity as a Gaussian density and the stream function as the solution to the associated Poisson equation.

Section 3 investigates the degenerate stochastic process introduced by Kolmogorov in 1934, alongside the associated Fokker–Planck equation and its solution proposed by Kolmogorov. Further connections between the Kolmogorov and Klein–Kramers equations are explored in Section 4. To start with, Section 3 summarizes Kolmogorov’s original findings. Surprisingly, the Fokker–Planck equation, as used by Kolmogorov in his seminal paper, is inconsistent with his initial assumptions regarding the underlying process. Moreover, his proposed solution has dimensional inconsistencies and, as a result, does not satisfy the Fokker–Planck equation and initial conditions. However, there is a silver lining; Kolmogorov’s solution can be corrected via several complementary methods, which the section outlines. It concludes with an example of a representative corrected solution to the Kolmogorov problem.

Section 4 explores a selection of representative affine stochastic processes in statistical physics. First, it introduces the Langevin equation, which describes the dynamics of an underdamped Brownian particle in a potential field. Following this, it derives the Klein–Kramers equation, capturing the probabilistic aspects of the motion of such a particle. It turns out that the Kolmogorov equation derived in Section 3 is a particular case of the Klein–Kramers equation. The section presents Chandrasekhar’s solutions to the Klein–Kramers equations describing free and harmonically bound particles. The Klein–Kramers equation is inherently degenerate, with white noise impacting the particle’s velocity but not its position. It is shown in Section 8 that many path-dependent problems share this characteristic in mathematical finance. For instance, financial variables like the geometric price averages, which serve as the underlying instruments for a particular class of Asian options, can be conceptualized as path integrals, fitting into the category of degenerate stochastic processes.

Section 5 describes backward (Kolmogorov) and forward (Fokker–Planck) equations for t.p.d.fs of multidimensional stochastic jump-diffusion processes. The section explains the significance of studying t.p.d.fs. It sets up the general framework for Kolmogorov and Fokker–Planck equations and identifies the subset of affine stochastic processes amenable to analysis using the Kelvin-wave formalism. Subsequently, the section introduces an augmentation technique, providing a natural approach to tackle degenerate problems. Finally, it illustrates methods for transforming specific nonaffine processes into affine form through coordinate transformations, enhancing the scope of problems accessible by the Kelvin-wave methodology.

Section 6 studies Gaussian stochastic processes. It introduces a general formula for regular Gaussian processes, accommodating both degenerate scenarios and nondegenerate cases, as in Kolmogorov’s example. It expands this formula to address the practically significant scenario of killed Gaussian processes, followed by several illustrative examples. Then, the section presents the derivation of the t.p.d.f. for the Kolmogorov process with time-varying coefficients and explores the OU process with time-dependent coefficients and its extension, the augmented OU process, which models the combined dynamics of the process and its integral. Although the results are classical, their derivation through Kelvin-wave expansions provides a novel and enriching angle, offering an alternative viewpoint for understanding and deriving these established results. Next, the section examines free and harmonically bound particles, contrasting the Kelvin-wave method with Chandrasekhar’s classical approach. Finally, it revisits the basic concepts introduced in Section 2, demonstrating the akin nature of the temporal-spatial evolution of vorticity in the two-dimensional flow of a viscous fluid to the dynamics of a harmonically bound particle. This finding is intriguing and unexpected, forging a connection between seemingly unrelated physical phenomena.

Section 7 considers non-Gaussian processes. It starts with a general formula for non-Gaussian dynamics, accommodating degenerate and nondegenerate processes. Then, it expands this formula to killed processes. Several interesting examples are studied. These examples include a Kolmogorov process driven by anomalous diffusion, Feller processes with constant and time-dependent coefficients, and degenerate and nondegenerate augmented Feller processes. A novel method for investigating finite-time explosions of t.p.d.fs for augmented Feller processes is developed as a helpful by-product of the analysis. In addition, arithmetic Brownian motions with path-dependent volatility and degenerate and nondegenerate arithmetic Brownian motions with stochastic volatility are analyzed in detail.

Section 8 illustrates the application of the methodology to financial engineering. To start with, it lays the foundation of financial engineering, providing a primer for the uninitiated. Then, the section introduces the geometric Brownian motion, a staple in financial modeling, and discusses the modifications necessary to reflect the complexities of financial markets better. Several traditional models, such as Bachelier, Black–Scholes, Heston, and Stein–Stein models, and a novel path-dependent volatility model are explored via the Kelvin-wave formalism. In addition, it is shown how to price Asian options with geometric averaging via the Kolmogorov’s solution described in Section 3. Besides, volatility and variance swaps and swaptions, bonds and bond options are investigated by linking financial formulas to those used in physics for underdamped Brownian motion.

Section 9 succinctly outlines potential future expansions of the work presented in this Element and summarizes the conclusions. Finally, this Element is a revised and expanded version of Reference LiptonLipton (2023).

A note on notation: Given the wide-ranging scope of this Element, from hydrodynamics to molecular physics, probability theory, and financial engineering, adopting a unified notation system is impractical. Each field has its conventions carved in stone, leading to inevitable variations in notation. Notation is designed for consistency within and, where possible, across sections. However, readers are encouraged to remain vigilant to maintain coherence in their understanding.

2 Fluid Flows

2.1 Euler and Navier–Stokes Equations

Hydrodynamics studies how fluids (liquids and gases) move, primarily relying on fluid motion’s fundamental equations: the Euler and Navier–Stokes equations, with the Euler equations applicable to inviscid (frictionless) flow and the Navier–Stokes equations describing viscous fluids. Hydrodynamics has numerous applications across various fields, including engineering, astrophysics, oceanography, and climate change, among many others.

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

Vt+(V)V+Pρ=0,V=0;(2.1)

where t is time, x is the position, Vt,x is the velocity vector, Pt,x is the pressure, ρ is the constant density, is the gradient, and denotes the scalar product; see, for example, Reference ChandrasekharChandrasekhar (1961). In Cartesian coordinates, the equations in (2.1) can be written as follows:

Vit+VjVixj+xiPρ=0,Vixi=0.(2.2)

Here and in what follows, Einstein’s summation convention over repeated indices is used.

The motion of the incompressible viscous fluid is described by the classical Navier–Stokes equations of the form:

Vt+(V)VνΔV+Pρ=0,V=0;(2.3)

where ν is the kinematic viscosity; see, for example, Reference ChandrasekharChandrasekhar (1961). Explicitly,

Vit+VjVixjν2Vixjxj+xiPρ=0,Vixi=0.(2.4)

The diffusive term νΔV in (2.4) describes frictions ignored in (2.3). Due to their greater generality, the Navier–Stokes equations are fundamental to understanding important phenomena, such as the transition from laminar to turbulent flow.

2.2 Linear Flows

This section studies exact solutions of the Euler and Navier–Stokes equations known as linear flows. These solutions are valuable for several reasons: (a) exact solutions provide precise, analytical descriptions of fluid flow patterns under specific conditions; (b) they serve as benchmarks for understanding fundamental hydrodynamics phenomena like wave propagation; (c) they provide a bridge which is crucial for more complex studies by simplifying the inherently complex and nonlinear nature of hydrodynamics, and making it possible to understand the behavior of more general fluid flows. Linear solutions of the Euler and Navier–Stokes equations help to study fluid flow stability. This understanding is crucial in predicting and controlling flow behavior in various engineering applications, from aerospace to hydraulic engineering. By starting with linear solutions, one can incrementally introduce nonlinear effects, allowing for a systematic study of nonlinear phenomena in hydrodynamics. This approach can uncover the mechanisms behind complex flows, including turbulence and chaotic flow behaviors. Exact linear solutions of the Euler equations provide a clear, analytical framework for exploring the behavior of fluids and validating more complicated models.

It is easy to show that the equations in (2.1) have a family of solutions Vt,x,Pt,x, linearly depending on spatial coordinates:

Vt,x=Ltx,   Pt,xρ=P0ρ+12Mtxx,(2.5)

where the 3×3 matrices Lt, Mt, are such that

dLtdt+L2t+Mt=0,TrLt=0,   Mt=M*t.(2.6)

It is clear that linear flows, given by (2.5), are unaffected by viscosity, hence they satisfy (2.14).

Flows (2.5) have stagnation points at the origin. Typical examples are planar flows of the form

V1=12sx1wx2,   V2=12wx1sx2,   V3=0,Pρ=P0ρ+14w2s2x12+x22.(2.7)

These flows are elliptic when s<w, and hyperbolic otherwise; see, for example, Reference Friedlander and Lipton-LifschitzFriedlander and Lipton-Lifschitz (2003).

2.3 Kelvin Waves in an Incompressible Fluid

The study of small perturbations of exact solutions of the Euler and Navier–Stokes equations is the core of the stability analysis in fluid dynamics. Examining their behavior is essential for predicting how fluid flows evolve under slight disturbances. One can determine whether a particular flow is stable or unstable by introducing small perturbations to an exact solution and observing the system’s response. If these perturbations grow over time, the flow is considered unstable; if they decay or remain bounded, the flow is stable. One of this analysis’s most critical applications is understanding the transition from laminar (smooth and orderly) to turbulent (chaotic and unpredictable) flows. Small perturbations can exhibit exponential growth, leading to the onset of turbulence. For more detailed investigations, direct numerical simulations of the perturbed Navier–Stokes equations can be used to study the nonlinear evolution of perturbations. This approach can capture the complete transition from initial instability to fully developed turbulence, offering insights into the complex interactions that drive flow dynamics. The study of perturbations offers theoretical insights into the fundamental nature of fluid dynamics, including the mechanisms of flow instability, transition, and turbulence structure. It helps in developing reduced-order models and theories that explain complex fluid phenomena. Here, Kelvin waves are used as the primary tool for studying small perturbations of linear flows. In the rest of this Element, Kelvin waves are used for other purposes. This section is dedicated to their brief description.

It is necessary to study the behavior of perturbations of solutions given by (2.5), which are denoted by vt,x,pt,x. By neglecting the quadratic term (v)v, one can write the system of PDEs for v,p as follows:

vt+(Ltx)v+Ltv+pρ=0,v=0.(2.8)

It has been known for a long time that linear PDEs (2.8) have wavelike solutions of the form:

vt,x,pt,xρ=at,atexpiβtxrt,(2.9)

where at,at are time-dependent amplitudes, and βt is the time-dependent wave vector; see Reference KelvinKelvin (1887), Reference OrrOrr (1907), Reference Craik and CriminaleCraik and Criminale (1986), and Reference Friedlander and Lipton-LifschitzFriedlander and Lipton-Lifschitz (2003). In this Element, these solutions are called the Kelvin waves. It should be emphasized that the so-called affine ansatz is a special instance of Kelvin wave. This observation allows one to discover similarities among seemingly unrelated topics, which, in turn, facilitates their holistic and comprehensive study. An excerpt from Kelvin’s original paper is shown in Figure 2.

Figure 2 An excerpt from Kelvin’s original paper, where Kelvin waves are introduced for the first time; see Reference KelvinKelvin (1887). Public domain.

As one can see from Figure 2, Kelvin considered the special case of the so-called shear linear flow of the form

Vt,x=V1x2,0,0=l12x2,0,0,(2.10)

between two plates, x2=0 and x2=L, the first one at rest and the second one moving in parallel.

The triplet rt, βt, at satisfies the following system of ODEs:

drtdtLtrt=0,   r0=r0,dβtdt+L*tβt=0,   β0=β0,datdt+Ltat2Ltatβtβtβtβt=0, a0=a0, β0a0=0.(2.11)

Here and in what follows, the superscript * stands for transpose. The corresponding pt can be found via the incompressibility condition. It is easy to show that for t0,

βtrt=β0r0,   βtat=0.(2.12)

Thus, the Kelvin-wave formalism results in ingenious separation of variables and allows us to solve a system of ODEs (2.11), rather than PDEs (2.8).

Typically, the equations in (2.11) are used to study the stability of the linear flow. Such a flow is unstable whenever at for some choices of β0,a0; see Reference BaylyBayly (1986), Reference LifschitzLifschitz (1995), and Reference Bayly, Holm and LifschitzBayly et al. (1996). Moreover, it can be shown that the same instabilities occur in general three-dimensional flows, because locally they are equivalent to linear flows; see Reference Lifschitz and HameiriLifschitz and Hameiri (1991a), Reference Friedlander and VishikFriedlander and Vishik (1991), Reference Lifschitz and HameiriLifschitz and Hameiri (1991b), and Reference Friedlander and Lipton-LifschitzFriedlander and Lipton-Lifschitz (2003).

Interestingly, Reference ChandrasekharChandrasekhar (1961) pointed out that the superposition of the linear flow (2.5) and the Kelvin wave (2.9), namely,

V˜t,x=Ltx+vt,x,P˜t,xρ=12Mtxx+pt,xρ,(2.13)

satisfies the nonlinear Euler equations (2.1) since the nonlinear term (v)v vanishes identically due to incompressibility.Footnote 3 Studying secondary instabilities of flows with elliptic streamlines, that is, instabilities of Kelvin waves is an important and intricate topic; see Reference Fabijonas, Holm and LifschitzFabijonas et al. (1997).

Viscosity does affect small perturbations of linear flows. For viscous incompressible fluids, Kelvin waves are governed by the following equations:

vt+(Ltx)v+LtvνΔv+pρ=0,v=0.(2.14)

The viscous version of (2.11) has the following form; see Reference LifschitzLifschitz (1991):

drtdtLtrt=0,    r0=r0,dβtdt+L*tβt=0,   β0=β0,datdt+Ltat2Ltatβtβtβtβt+νβt2at=0, a0=a0,β0a0=0.(2.15)

It is shown in Section 6.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.

The evolution of a typical Kelvin wave parameters triplet rt, βt, at is illustrated in Figure 3. The impact of viscosity is illustrated in Figure 4. These figures show that depending on the initial orientation of the wave vector βt, the amplitude at can be either bounded or unbounded. For elliptic flows, unbounded amplitudes are always present for specific orientations, so all of them are unstable; see Reference BaylyBayly (1986), Reference Bayly, Holm and LifschitzBayly et al. (1996), Reference Friedlander and Lipton-LifschitzFriedlander and Lipton-Lifschitz (2003), and references therein.

Figure 3 Kelvin waves corresponding to two different orientations of the initial wave vector β0 and a0. (a), (b) β(0)=sinπ/4,0,cosπ/4, a0=0,sinπ/4,0; (c), (d) β(0)=sinπ/3,0,cosπ/3, a0=0,sinπ/3,0. Other parameters are as follows: T=100, ω=1, s=0.5. In the first case, at stays bounded, while at explodes in the second case. This explosion means that the underlying elliptic flow is unstable. Author’s graphics.

Figure 4 Kelvin waves in the viscous fluid with viscosity ν=0.07. Other parameters and initial conditions are the same as in Figure 3. Viscosity dampens the instability but, generally, does not suppress it entirely. Author’s graphics.

3 Kolmogorov Stochastic Process

3.1 Background

The Kolmogorov equation studies the evolution of a particle in the phase space. The particle’s position and velocity evolve in time due to the interplay between the deterministic drift and stochastic force affecting only its velocity. Since only the particle’s velocity is affected by the random force, the PDE describing the evolution of the t.p.d.f. in the phase space is degenerate. The Kolmogorov equation is a particular case of the Klein–Kramers equation studied in Section 4.

The significance of the Kolmogorov equation lies in its ability to model the intricate balance between deterministic behavior and stochastic dynamics, providing a basic framework for studying the evolution of systems in phase space. It has important applications in various fields, including physics for understanding particle dynamics, finance for modeling asset prices, and beyond. It demonstrates the profound interplay between stochastic processes and differential equations.

The Kolmogorov equation is hypoelliptic; as such, it serves as a prototype for a broad class of hypoelliptic PDEs. Although it does not meet the exact criteria for ellipticity (due to the second-order derivatives not being present in all directions of the phase space), the solutions to the equation are still smooth, which is particularly important in the context of stochastic processes, where hypoellipticity ensures that the probability density function remains smooth and well-behaved, facilitating the analysis of the system’s dynamics over time.

3.2 Summary of Kolmogorov’s Paper

In a remarkable (and remarkably concise) note, Kolmogorov considers a system of particles in n-dimensional space with coordinates q1,,qn, and velocities q˙1,,q˙n, assuming the probability density function

gt,q1,,qn,q˙1,,q˙n,t,q1,,qn,q˙1,,q˙n

exists for some time t>t, and reveals (without any explanation) an analytical expression for g in the one-dimensional case; see Reference KolmogoroffKolmogoroff (1934).Footnote 4 This note is the third in a series of papers, the previous two being Reference KolmogoroffKolmogoroff (1931), Reference Kolmogoroff(1933).

Kolmogorov makes the following natural assumptions:

EΔqiq˙iΔt=oΔt,(3.1)
EΔqi2=oΔt,(3.2)

where Δt=tt. Equations (3.1) and (3.2) imply

EΔqi=q˙iΔt+oΔt,(3.3)
EΔqiΔqjEΔqi2EΔqj2=oΔt.(3.4)

Furthermore, under very general assumptions, the following relationships hold:

EΔq˙i=fit,q,q˙Δt+oΔt,(3.5)
EΔq˙i2=kiit,q,q˙Δt+oΔt,(3.6)
EΔq˙iΔq˙j=kijt,q,q˙Δt+oΔt,(3.7)

where f and k are continuous functions. Equations (3.2), and (3.6) imply

EΔq˙iΔq˙jEΔq˙i2EΔq˙j2=oΔt.(3.8)

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

gt=q˙igqiq˙ifit,q,q˙g+2q˙iq˙jkt,q,q˙g.(3.9)

In the one-dimensional case, one has

gt=q˙gqq˙ft,q,q˙g+2q˙2kt,q,q˙g.(3.10)

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

When f and k are constants, (3.10) becomes

gt=q˙gqfgq˙+k2gq˙2.(3.11)

The corresponding fundamental solution of has the following form:

g=23πk2tt2expq˙q˙ftt24ktt3qqq˙+q˙2tt2k3tt3.(3.12)

One can see that Δq˙ is of the order Δt1/2. At the same time

Δq=q˙Δt+OΔt3/2.(3.13)

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

Kolmogorov’s original paper is shown in Figure 5.

Reproduced by kind permission of the Editors of Annals of Mathematics.

Figure 5 Kolmogorov’s original paper, presented here for the inquisitive reader to enjoy; see Reference KolmogoroffKolmogoroff (1934).

Kolmogorov equations fascinated mathematicians for a long time and generated a great deal of research; see, for example, Reference WeberWeber (1951), Reference HörmanderHörmander (1967), Reference KuptsovKuptsov (1972), Reference Lanconelli, Pascucci, Polidoro, Sh, Birman, Solonnikov and UraltsevaLanconelli et al. (2002), Reference Pascucci and Parabolic ProblemsPascucci (2005), Reference Ivasishen and MedynskyIvasishen and Medynsky (2010), and Reference Duong and TranDuong & Tran (2018), among others.

It is worth mentioning that physicists derived Equations (3.9) and (3.10) at least a decade earlier than Kolmogorov; see Section 4.

3.3 Challenge and Response

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

First, Equations (3.9) and (3.10) are not the Fokker–Planck equations associated with the process described by Equations (3.5)–(3.7), since they lack the prefactor 1/2 in front of the diffusion terms. The corrected multivariate equation has the form

gt=q˙iqigq˙ifit,q,q˙g+122q˙iq˙jkt,q,q˙g,(3.14)

while the corresponding one-dimensional equation has the form

gt=q˙qgq˙ft,q,q˙g+122q˙2kt,q,q˙g.(3.15)

Alternatively, Equations (3.6) and (3.7) can be altered as follows:

EΔq˙i2=2kiit,q,q˙Δt+oΔt,(3.16)

EΔq˙iΔq˙j=2kijt,q,q˙Δt+oΔt.(3.17)

In the following discussion, the Fokker–Planck equation is updated.

Second, g given by (3.12) does not solve (3.10). It also does not satisfy the (implicit) initial condition

gt,q,q˙,t,q,q˙=δqqδq˙q˙,(3.18)

where δ. is the Dirac δ-function. The fact that expression (3.12) does not solve (3.10) 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:

t=t=T,  q=q=L,  q˙=q˙=LT,  g=TL2,  f=LT2,k=L2T3.(3.19)

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

gλ2t,λ3q,λq˙,λ2t,λ3q,λq˙;λ1f,k=λ4gt,q,q˙,t,q,q˙;f,k.(3.20)

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

3qqq˙+q˙2tt2k3tt3

in the exponent is not nondimensional, as it should be, and has dimension T6L1, while the prefactor

23πk2tt2

has dimension T4L1, instead of the right dimension, TL2.

Third, due to yet another typo, the solution given by (3.12) does not converge to the initial condition in the limit tt. Indeed, asymptotically, one has

gHk3tt36,qqH2ktt,q˙q˙4δqqδq˙q˙,(3.21)

where Hμ,ν is the standard heat kernel:

Hμ,ν=expν22μ2πμ.(3.22)

However, not all is lost. Dimensional analysis shows that the correct solution gt,q,q˙,t,q,q˙;f,k of (3.10) has the following form:

g=32πktt2expq˙q˙ftt24ktt3qqq˙+q˙2tt2ktt3,(3.23)

which is not far from Kolmogorov’s formula. Similarly, the correct solution of (3.15) has the following form:

g=3πktt2expq˙q˙ftt22ktt6qqq˙+q˙2tt2ktt3.(3.24)

3.4 Direct Verification

In order to avoid confusion, from now on, the notation is changed to make the formulas easier to read. Specifically, it is assumed that xˉ represents the position of a particle at time tˉ and x its position at time t, while yˉ represents its velocity at time tˉ, and y its velocity at time t, so that

t,q,q˙t,x,y,   t,q,q˙tˉ,xˉ,yˉ.(3.25)

One of our objectives is deriving the (corrected) Kolmogorov formula from first principles using Kelvin waves. Subsequently, it is shown how to use it in the financial mathematics context. The governing SDE can be written as

dxˆt=yˆtdt,   xˆt=x,dyˆt=bdt+σdWˆt,   yˆt=y.(3.26)

The corresponding Fokker–Planck–Kolmogorov problem for the t.p.d.f. ϖt,x,y,tˉ,xˉ,yˉ has the form:

ϖtˉt,x,y,tˉ,xˉ,yˉ12σ2ϖyˉyˉt,x,y,tˉ,xˉ,yˉ+yˉϖxˉt,x,y,tˉ,xˉ,yˉ+bϖyˉt,x,y,tˉ,xˉ,yˉ=0,ϖt,x,y,t,xˉ,yˉ=δxˉxδyˉy.(3.27)

The solution of (3.27) is as follows:

ϖt,x,y,tˉ,xˉ,yˉ=3πσ2T2expΦt,x,y,tˉ,xˉ,yˉ,(3.28)

where

Φt,x,y,tˉ,xˉ,yˉ=yˉybT22σ2T+6xˉxyˉ+yT22σ2T3=A22+6B2,(3.29)

and

A=yˉybTσ2T,  A=1,   B=xˉxyˉ+yT2σ2T3,   B=1.(3.30)

Here and in what follows, the following shorthand notation is used:

T=tˉt.(3.31)

Let us check that ϖ satisfies the Fokker–Planck equation and the initial conditions. A simple calculation yields:

Φtˉ=A22T+bAσ2T+18B2T+6yˉ+yBσ2T3,Φxˉ=12Bσ2T3,   Φyˉ=A6Bσ2T,   Φyˉyˉ=4σ2T,(3.32)
ϖtˉϖ=2TΦtˉ,   ϖxˉϖ=Φxˉ,   ϖyˉϖ=Φyˉ,   ϖyˉyˉϖ=Φyˉyˉ+Φxˉ2,(3.33)

so that

ϖtˉK12σ2ϖyˉyˉK+yˉϖxˉK+bϖyˉK=ϖK2TΦtˉ+12σ2ΦyˉyˉΦyˉ2yˉΦxˉbΦyˉ=ϖK2T+A22T+bAσ2T+18B2T+6yˉ+yBσ2T3+2TA6B22T12yˉBσ2T3bA6Bσ2T=0.(3.34)

When T0, one has the following asymptotic expression:

ϖKt,x,y,tˉ,xˉ,yˉHσ2T312,xˉxHσ2T,yˉyδxˉxδyˉy.(3.35)

3.5 Solution via Kelvin Waves

Now, Kolmogorov’s formula is derived by using Kelvin waves (or an affine ansatz), which requires solving the problem of the following form:

Ktˉt,x,y,tˉ,xˉ,yˉ,k,l12σ2Kyˉyˉt,x,y,tˉ,xˉ,yˉ,k,l+yˉKxˉt,x,y,tˉ,xˉ,yˉ,k,l+bKyˉt,x,y,tˉ,xˉ,yˉ,k,l=0,Kt,xˉ,yˉ,t,x,y,k,l=expikxˉx+ilyˉy.(3.36)

Here

k=1L,   l=TL,  K=1.(3.37)

By using the well-known results concerning the inverse Fourier transform of the δ-function, one gets the following expression for the t.p.d.f. ϖt,x,y,tˉ,xˉ,yˉ:

ϖt,x,y,tˉ,xˉ,yˉ=12π2Kt,x,y,tˉ,xˉ,yˉ,k,ldkdl.(3.38)

To calculate K, one can use the affine ansatz and represent it in the following form:

Kt,x,y,tˉ,xˉ,yˉ,k,l=expΨt,x,y,tˉ,xˉ,yˉ,k,l,(3.39)

where

Ψt,x,y,tˉ,xˉ,yˉ,k,l=αt,tˉ+ikxˉx+iγt,tˉyˉily.(3.40)

and

KtˉK=Ψtˉ=αtˉt,tˉ+iγtˉt,tˉyˉ,   KxˉK=Ψxˉ=ik,   KyˉK=Ψyˉ=iγt,tˉ,   KyˉyˉK=Ψyˉ2=γ2t,tˉ.(3.41)

Accordingly,

αtˉt,tˉ+12σ2γ2t,tˉ+iγtˉt,tˉyˉ+ikyˉ+ibγt,tˉ=0,αt,t=0,   γt,t=l,(3.42)

so that

αtˉt,tˉ+12σ2γ2t,tˉ+ibγt,tˉ=0,   αt,t=0,γtˉt,tˉ+k=0,   γt,t=l.(3.43)

Straightforward calculation shows that:

γt,tˉ=kT+l,αt,tˉ=12σ2k2T33klT2+l2TibkT22+lT.(3.44)

Equations (3.38), (3.39), (3.40) and (3.44) yield

ϖt,x,y,tˉ,xˉ,yˉ=12π2exp12σ2k2T33klT2+l2T+ ikxˉxyˉT+bT22+ilyˉybTdkdl.(3.45)

It is clear that ϖt,x,y,tˉ,xˉ,yˉ can be viewed as the characteristic function of the Gaussian density in the k,l space, evaluated at the point xˉxyˉT +bT22,yˉybT:

ϖt,x,y,tˉ,xˉ,yˉ=detC1/22πGT,k,l×expikxˉxyˉT+bT22+ilyˉybTdkdl,(3.46)

where

GT,k,l=12πdetC1/2exp12klC1Tkl,(3.47)

and

CT=12σ2T36σ2T26σ2T24σ2T,detCT=12σ4T4.(3.48)

As before, denotes the scalar product. Accordingly,

ϖt,x,y,tˉ,xˉ,yˉ=3πσ2T2expΩt,x,y,tˉ,xˉ,yˉ,(3.49)

where

Ωt,x,y,tˉ,xˉ,yˉ==12xˉxyˉT+bT22yˉybTCxˉxyˉT+bT22yˉybT=6xˉxyˉT+bT222σ2T3+6xˉxyˉT+bT22yˉybTσ2T2+2yˉybT2σ2T=A22+6B2=Φtˉ,xˉ,yˉ,x,y,(3.50)

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

Note that the t.p.d.f. ϖ is a bivariate Gaussian distribution. Completing the square, one can write

Φ=yˉybT22σ2T+6xˉxyˉ+yT22σ2T3=6σ2T3xˉp26σ2T2xˉpyˉq+2σ2Tyˉq2,(3.51)

and represent ϖ the form:

ϖt,x,y,tˉ,xˉ,yˉ=exp121ρ2xˉp2σx22ρxˉpyˉqσxσy+yˉq2σy22πσxσy1ρ2,(3.52)

where

σx=σ2T33,   σy=σ2T,   ρ=32,p=x+yT+bT22,  q=y+bT.(3.53)

Equation (3.28) can be derived by using the Hankel transform. Since

σ2C1T=T33T22T22T=T3/22T1/22T3/2233T1/22*T3/22T1/22T3/2233T1/22,(3.54)

one can introduce

kˉlˉ=T3/22T1/22T3/2233T1/22kl,kl=3T3/23T3/2T1/23T1/2kˉlˉ,(3.55)

and rewrite (3.46) as follows:

ϖt,x,y,tˉ,xˉ,yˉ=32π2T2exp12σ2kˉ2+lˉ2+ikˉlˉ3T3/2T1/23T3/23T1/2xˉxyˉT+bT22yˉybTdkˉdlˉ=32π2T2exp12σ2kˉ2+lˉ2+ikˉ3xˉxyˉT+bT22+yˉybTTT3/2+ilˉ3xˉxyˉT+bT22+3yˉybTTT3/2dkˉdlˉ.(3.56)

Thus, ϖt,x,y,tˉ,xˉ,yˉ is the Fourier transform of a radially symmetric function of kˉ,lˉ*. Accordingly, it can be calculated via the Hankel transform of the function expσ2rˉ2/2:

ϖt,x,y,tˉ,xˉ,yˉ=3πT2H0eσ2rˉ22sˉ=3πσ2T2esˉ22σ2,(3.57)

where

rˉ2=kˉ2+lˉ2,sˉ2=43xˉxyˉT+bT222+3xˉxyˉT+bT22yˉybTT+yˉybT2T2T3.(3.58)

See, for example, Reference Piessens and PoularikasPiessens (2000). As expected, the corresponding expression coincides with the one given by (3.52).

3.6 Solution via Coordinate Transform

This section briefly considers the method of coordinate transformations, reducing the original Fokker–Planck equation for the Kolmogorov problem to a Fokker–Planck equation with spatially independent coefficients. To this end, the following ansatz is used:

x˜,y˜=xˉtˉtyˉ,yˉ.(3.59)

This choice is explained in more detail in Section 6. Straightforward calculation yields

tˉ=tˉy˜x˜,   xˉ=x˜,   y=tˉtx˜+y˜,(3.60)

so that (3.27) becomes

tˉy˜x˜ϖt,x,y,tˉ,x˜,y˜12σ2tˉtx˜+y˜2ϖt,x,y,tˉ,x˜,y˜+y˜ϖx˜t,x,y,tˉ,x˜,y˜+btˉtx˜+y˜ϖt,x,y,tˉ,x˜,y˜=0,ϖt,x,y,t,x˜,y˜=δx˜xδy˜y.(3.61)

Further calculations show that coefficients of the preceding equation are spatially independent:

tˉϖt,x,y,tˉ,x˜,y˜12σ2tˉtx˜+y˜2ϖt,x,y,tˉ,x˜,y˜+btˉtx˜+y˜ϖt,x,y,tˉ,x˜,y˜=0,ϖt,x,y,t,x˜,y˜=δx˜xδy˜y.(3.62)

Accordingly, one can use the classical Fourier transform and represent the solution of (3.62) in the form

ϖt,x,y,tˉ,x˜,y˜=12π2exp12σ2k2T33klT2+l2T+ikx˜x+bT22+ily˜ybTdkdl,(3.63)

similar to (3.45). Thus, ϖ has the form given by (3.52) with xˉ,yˉ replaced by x˜,y˜. The exact form is recovered once x˜,y˜ are expressed in terms of xˉ,yˉ by virtue of (3.59).

3.7 A Representative Example

A typical solution of the Kolmogorov equation is illustrated in Figure 6. This figure clearly shows that there is a good agreement between a Monte Carlo simulation of the stochastic process xˆt,yˆt given by the equations in (3.26) and the corrected Kolmogorov formula (3.24).

Figure 6 A thousand trajectories of a typical Kolmogorov process. Parameters are as follows: T=5, dt=0.01, f=0.2, σ=0.8. (a) xt, (b) yt, (c) xˉT,yˉT, (d) contour lines of ϖ0,0,0,T,x˜,y˜. Author’s graphics.

4 Klein–Kramers Stochastic Process

4.1 Background

The Klein–Kramers equation plays a vital role in statistical physics by offering a detailed mathematical framework for studying the dynamics of particles in a viscous, random medium. Specifically, it describes the evolution of the t.p.d.f. of a particle’s momentum and position in the phase plane, accounting for deterministic forces arising from potential and stochastic thermal forces arising from random collisions with the medium’s molecules. This equation is particularly important for studying nonequilibrium systems, which cannot be analyzed via traditional equilibrium statistical mechanics tools. By incorporating frictional forces, which tend to dampen the motion of particles, potential forces, which push them deterministically, and random thermal forces, which inject randomness into the system, the Klein–Kramers equation bridges the gap between microscopic laws of motion and the macroscopic observable phenomena, such as diffusion, thermal conductivity, and viscosity. Moreover, the Klein–Kramers equation serves as a foundation for exploring more complex phenomena in nonequilibrium statistical mechanics, including the study of transition state theory in macrokinetics of chemical reactions, the behavior of particles in external fields, and the exploration of noise-induced transitions and stochastic resonance in physical and biological systems. It also arises in financial engineering, for instance, in pricing volatility and variance swaps.

4.2 Langevin Equation

Start with the Langevin equation for particles moving in a potential field and impacted by random forces; see Reference LangevinLangevin (1908). This section uses the standard notation, rather that the original notation used in Reference ChandrasekharChandresekhar (1943). Hopefully, the diligent reader will not be easily confused. The stochastic Langevin equation describes the evolution of systems under the influence of deterministic forces and random fluctuations. Because of its versatility, it is widely used in physics and other disciplines to model the dynamics of particles subjected to systematic forces derived from potential energy and random forces representing thermal fluctuations. This equation describes a particle experiencing frictional resistance proportional to its velocity (a deterministic component) and random kicks from the surrounding molecules (a stochastic component capturing the essence of Brownian motion). The Langevin equation thus provides a robust framework for studying the behavior of systems subject to noise, enabling insights into phenomena such as diffusion, thermal equilibrium, and the statistical properties of microscopic systems.

Consider an underdamped Brownian particle. In contrast to the standard Brownian motion, which is overdamped, it is assumed that the frictions are finite, so that one must treat the particle’s velocity as an independent degree of freedom. Hence, the particle’s state is described by a pair x,y, where x and y are its position and velocity, respectively. Consider a d-dimensional space, with d=1 and d=3 of particular interest, and write the corresponding Langevin equations in the following form:

dxˆtdt=yˆt,   xˆt=x,dyˆtdt=κyˆtVxˆtm+2κkBTmdWˆtdt,  yˆt=y,(4.1)

where Wˆt is a standard d-dimensional Wiener process. Here m is the particle mass, κ is the friction coefficient, kB is the Boltzmann constant, T is the temperature, Vx is the external potential, and dWˆt/dt is a d -dimensional Gaussian white noise. Below, the ratio κkBT/m is denoted as a.

Of course, one can rewrite the equations of (4.1) as a system of stochastic differential equations (SDEs):

dxˆt=yˆtdt,   xˆt=x,dyˆt=κyˆtdtVxˆtmdt+2adWˆt,   yˆt=y.(4.2)

For a 1-dimensional particle, (4.2) becomes:

dxˆt=yˆtdt,   xˆt=x,dyˆt=κyˆtdtVxxˆtmdt+2adWˆt,   yˆt=y.(4.3)

It is clear that the Kolmogorov equation (3.11) is a special case of (4.3) with κ=0, Vx=mfx, k=a.

4.3 Klein–Kramers Equation

Fokker, Planck, and their numerous followers derived and studied the forward parabolic equation for the t.p.d.f. ϖt,x,y,tˉ,xˉ,yˉ associated with a stochastic process. For the stochastic process governed by SDEs (4.2), the corresponding equation, called the Klein–Kramers equation, has the following form:

ϖtˉaϖyˉyˉ+yˉϖxˉ((κyˉ+V(xˉ)m)ϖ)yˉ=0,ϖ(t,x,y,t,xˉ,yˉ)=δ(xˉx)δ(yˉy).(4.4)

The backward parabolic Kolmogorov equation can be written as follows:

ϖt+aϖyy+yϖxκy+Vxmϖy=0,ϖtˉ,x,y,tˉ,xˉ,yˉ=δxˉxδyˉy.(4.5)

Details are given in Reference FokkerFokker (1914), Reference PlanckPlanck (1917), Reference KleinKlein (1921), Reference ChapmanChapman (1928), Reference KolmogoroffKolmogoroff (1931, Reference Kolmogoroff1933, Reference Kolmogoroff1934); Reference KramersKramers (1940), Reference ChandrasekharChandresekhar (1943), Reference RiskenRisken (1989), and Reference Hänggi, Talkner and BorkovecHänggi et al. (1990), as well as a multitude of subsequent sources. For fascinating historical details, see Reference Ebeling, Gudowska-Nowak and SokolovEbeling et al. (2008). The Klein–Kramers equation (occasionally called Klein–Kramers–Chandrasekhar equation) describes the dynamics of a particle’s probability distribution in phase space (position and momentum) for systems subjected to friction and random forces, typically at the mesoscopic scale. The Klein–Kramers equation provides a comprehensive framework for modeling and understanding complex systems far from equilibrium, linking microscopic physics with macroscopic observables. Accordingly, it is used in various fields, such as materials science, chemistry, and astrophysics, to predict the evolution of systems over time, accounting for both deterministic dynamics and the effects of randomness.

4.4 Chandrasekhar’s Solutions

In a well-known survey article, Reference ChandrasekharChandresekhar (1943) described elegant solutions of (4.4) for a free particle and a harmonically bound particle, which he derived by using ingenious changes of coordinates. For a free particle, Reference ChandrasekharChandresekhar (1943) writes the corresponding Klein–Kramers equation as follows:

ϖtˉaϖyˉyˉ+yˉϖxˉκyˉϖyˉκϖ=0,ϖt,x,y,t,xˉ,yˉ=δxˉxδyˉy.(4.6)

By using ingenious coordinate transforms, he shows that

ϖ=12πFGH21/2expFR22HRS+GS22FGH2,(4.7)

where

R=yˉeκTy,S=xˉx1eκTκy,F=aκ33+4eκTe2κT+2κT,G=aκ1e2κT,H=aκ21eκT2.(4.8)

Here the original Chandrasekhar’s notation is slightly changed to make the exposition more internally consistent.

Since it is assumed that stochastic drivers are uncorrelated, the t.p.d.f. ϖ3 can be presented as a product of three 1 -dimensional t.p.d.f. ϖ1:

ϖ3=ϖ11ϖ21ϖ31=18π3FGH23/2expFR22HRS+GS22FGH2,(4.9)

where R2=yˉ1eκTy12+yˉ2eκTy22+yˉ3eκTy32, and so on.

Chandrasekhar generalized (4.7) to the case of harmonically bound particles. We shall revisit Chandrasekhar’s formulas for free and bound particles later. While Reference ChandrasekharChandresekhar (1943) stopped at (4.7), for practical applications, it is more useful to represent the exponent as an explicit quadratic form of xˉ and yˉ, which is done in Section 6.5.

5 Transition Probability Densities for Stochastic Processes

5.1 Motivation

The problems considered in Sections 3 and 4 are used in what follows to develop a general theory. For that, one needs to know some foundational information about stochastic processes discussed in this section. Stochastic processes play a crucial role across various scientific disciplines, which is fundamental for modeling systems influenced by randomness and uncertainty. These processes are pivotal in fields ranging from physics and chemistry to biology, economics, and financial engineering. They help to understand phenomena where outcomes are not deterministic but probabilistic, capturing the dynamics of complex systems over time. The analysis of stochastic processes enables scientists and engineers to predict behavior, assess risks, and make informed decisions based on the likelihood of future events.

The backward Kolmogorov and forward Fokker–Planck equations offer a mathematical description of how systems evolve under the influence of stochastic factors. This capability to model the t.p.d.fs of diverse processes underlines the equations’ fundamental importance in scientific research and practical applications across disciplines.

The Kolmogorov and Fokker–Planck equations are adjoint partial differential equations that describe how the probability density of a system’s state evolves in time. The Kolmogorov equation focuses on calculating the expected value at a given time of random outcomes, which become known sometime in the future. Conversely, the Fokker–Planck equation is concerned with the evolution of the conditional probability density function of a process’s state at a future time, given its current state.

The Kolmogorov and Fokker–Planck equations are applied in physics and chemistry to study the random motion of particles in fluids, the statistical behavior of thermodynamic systems, and the kinetics of chemical reactions. In biology, these equations model population dynamics, genetic variation, and the spread of diseases, among other processes, providing insights into how randomness affects biological phenomena. In financial engineering, they are used to model the evolution of asset prices, interest rates, and other economic indicators, underpinning the valuation of derivatives and the management of financial risks.

5.2 Backward and Forward Equations

Start with a jump-diffusion process driven by the SDE of the form

dzˆt=bt,zˆtdt+σt,zˆtdWˆt+υdΠˆtt,zˆt,   zˆt=z,(5.1)

with smooth coefficients b,σ. This process is driven by the standard Wiener process Wˆt and the Poisson process Πˆtt,z with intensity λt,z such that

EdΠtt,zzˆt=z=λt,zdt,(5.2)

while υ is drawn from a distribution with density ϕυ,t,z, which (in general) is t,z-dependent.

More generally, it is possible to consider the so-called general compound or marked Poisson processes, such that υ=υt,z,q , where υ is monotonic in z, and q is a random mark variable drawn from a distribution with density ϕq,t,z, which (in general) is t,z-dependent. However, since this Element is interested in a particular class of stochastic processes, solvable via Kelvin waves ansatz this generalization is not particularly useful.

It is well-known that for suitable test functions u˜z the expectation

ut,z=Eu˜zˆtˉzˆt=z(5.3)

solves the following integro-differential backward Kolmogorov problem:

utt,z+at,zuzzt,z+bt,zuzt,z+λt,zut,z+υϕυ,t,zdυλt,zut,z=0,utˉ,z=u˜z,(5.4)

where

at,z=12σ2t,z.(5.5)

In particular, the t.p.d.f. ϖt,z,tˉ,zˉ such that

Probzˉ<zˆtˉ<zˉ+dzˉzˆt=z=ϖt,z,tˉ,zˉdzˉ,(5.6)

solves the following backward Kolmogorov problem:

ϖtt,z+at,zϖzzt,z+bt,zϖzt,z+λt,zϖt,z+υϕυ,t,zdυλt,zϖt,z=0,ϖtˉ,z,tˉ,zˉ=δzzˉ.(5.7)

It is possible to derive a forward problem for ϖt,z,tˉ,zˉ, which ϖ satisfies as a function of tˉ,zˉ, which is called Fokker–Planck or forward Kolmogorov problem. This problem has the following form:

ϖtˉtˉ,zˉatˉ,zˉϖtˉ,zˉzˉzˉ+btˉ,zˉϖtˉ,zˉzˉλtˉ,zˉυϖtˉ,zˉυϕtˉ,zˉυ,υdυ+λtˉ,zˉϖtˉ,zˉ=0,ϖt,z,t,zˉ=δzˉz.(5.8)

One can generalize backward Kolmogorov and forward Fokker–Planck equation to the multidimensional case. The underlying nz-dimensional process zˆt=zˆi,t, i=1,,nz, has the form

dzˆt=bt,zˆtdt+Σt,zˆtdWˆt+υdΠˆtt,zt,(5.9)

where Wˆt=Wˆj,t is an nW -dimensional Wiener process, j=1,,nW, and Πˆ=Πˆk,t is an nΠ-dimensional state-dependent Poisson process, k=1,,nΠ with intensity λ. The corresponding state-dependent coefficients are as follows:

bt,z=bit,z,   i=1,,nz,Σt,z=Σijt,z,    i=1,,nz,    j=1,,nW,λt,z=λit,z,   i=1,,nΠ,υ=υik,    i=1,,nz,    k=1,,nΠ,(5.10)

while υk are drawn from distributions with densities ϕkυ,t,z, which (in general) are t,z-dependent. Explicitly, the equations in (5.9) can be written as follows:

dzˆi,t=bit,zˆtdt+Σijt,zˆtdWˆj,t+υikdΠˆkt,zˆt.(5.11)

The backward and forward equations for the t.p.d.f. ϖ can be written as follows:

ϖtt,z+aijt,zϖzizjt,z+bit,zϖzit,z+λkt,zϖt,z+υkϕkυk,t,zdυkΛt,zϖt,z=0,ϖtˉ,z,tˉ,zˉ=δzzˉ,At,z=aiit,z=12Σt,zΣ*t,z=12σijt,zσijt,z,Λt,z=k=1nΠλkt,z.(5.12)

For the generic terminal condition u˜z, the corresponding backward problem has the following form:

utt,z+aijt,zuzizjt,z+bit,zuzit,z+λkt,zut,z+υkϕkυk,t,zdυkΛt,zut,z=0,utˉ,z=u˜z,(5.14)

The forward equations for the t.p.d.f. ϖ can be written as follows:

ϖtˉtˉ,zˉaijtˉ,zˉϖtˉ,zˉzˉizˉj+bitˉ,zˉϖtˉ,zˉzˉiλktˉ,zˉυkϖtˉ,zˉυkϕkυk,tˉ,zˉυkdυk+Λtˉ,zˉϖtˉ,zˉ=0,ϖt,z,t,zˉ=δzˉz.(5.15)

Further details can be found in Reference Bharucha-ReidBharucha-Reid (1960), Reference FellerFeller (1971), Reference Gihman and SkorohodGihman and Skorohod (1972), Reference ArnoldArnold (1974), and Reference HansonHanson (2007), among others.

Although, depending on the actual problem at hand, it might be preferable to work with either the backward or the forward problem, experience suggests that in the context of mathematical finance the backward problem is easier to deal with, not least because they are meaningful for the generic terminal value u˜zˉ.

Since the preceding definitions are very general, it is necessary to be more specific in defining the class of problems which can be solved by using Kelvin waves. Consider processes such that

At,z=A0t+ziAit,   bt,z=b0t+zibit,λt,z=λ0t+ziλit,   ϕυ,t,z=ϕυ,t,(5.16)

so that the corresponding backward Kolmogorov problem has the form

utt,z+aij0t+zlaijltuzizjt,z+bi0t+zlbiltuzit,z+λk0t+zlλkltut,z+υkϕkυk,tdυkΛ0t+zlΛltut,z=0,utˉ,z=u˜z.(5.17)

Symbolically, (5.17) can be written as follows:

utt,z+L0ut,z+l=1nzzlLlut,z=0,utˉ,z=u˜z,(5.18)

where L0, Ll are spatially homogeneous operators, with coefficients depending only on time (at most):

L0ut,z=aij0tuzizjt,z+bi0tuzit,z+λk0tut,z+υkϕkυk,tdυkΛ0tut,z,Llut,z=aijltuzizjt,z+biltuzit,z+λkltut,z+υkϕkυk,tdυkΛltut,z.(5.19)

For the t.p.d.f. ϖ, one has

ϖtt,z+L0ϖt,z+l=1nzzlLluut,z=0,ϖtˉ,z,tˉ,zˉ=δzzˉ.(5.20)

Moreover, to cover interesting and important cases, such as anomalous diffusions and the like, generalize (5.18) and consider pseudo-differential operators Llˉ, lˉ=0,,nz. Recall that a translationally invariant pseudo-differential operator L is defined as follows:

Luz=12πnzLmuzeimzzdzdm,(5.21)

where Lm is called the symbol of a pseudo-differential operator; see, for example, Reference CordesCordes (1995) and Reference WongWong (2014). It is clear that all diffusion operators belong to this category, and so do jump-diffusion operators. The symbol of the operator Llˉt

Llˉt,m=aijlˉtmimj+ibilˉtmi+λklˉtψkt,mΛlˉt,(5.22)

where ψkm is the characteristic function of ϕkυ:

ψkt,m=eimυkϕkt,υkdυk.(5.23)

While frequently studied in the pure and applied mathematical context, in the financial engineering context pseudo-differential operators are seldom discussed; see, however, Reference Jacob, Schilling and Barndorff-NielsenJacob and Schilling (2001).

By definition, Fourier and Kelvin modes are eigenfunctions of the operators L0, Ll. Accordingly, when all Ll=0, one can solve the corresponding backward problem via the standard Fourier modes F given by (1.1):

ut,z=12πnzu˜zeαt,tˉ,m+imzzdzdm,(5.24)

where

αtt,tˉ,m+L0t,m=0,   αtˉ,tˉ,m=0,(5.25)

so that

αt,tˉ,m=ttˉL0s,mds.(5.26)

However, in general, one needs to use Kelvin modes K, given by (1.2):

ut,z=12πnzu˜zeαt,tˉ,m+iδt,tˉ,mzimzdzdm,(5.27)

where

αtt,tˉ,m+L0t,δt,tˉ,m=0,   αtˉ,tˉ,m=0,δl,tt,tˉ,m+Llt,δt,tˉ,m=0,   δtˉ,tˉ,m=m.(5.28)

Of course, finding explicit solutions of ODEs (5.28) is possible only in exceptional cases, some of which are discussed below. However, it is always possible to solve them numerically, which is much easier than trying to solve the corresponding PDEs directly.

As mentioned earlier, three archetypal stochastic processes are arithmetic Wiener processes (or Brownian motions), Ornstein-Uhlenbeck (OU) and Feller processes; see Reference Uhlenbeck and OrnsteinUhlenbeck and Ornstein (1930), Reference ChandrasekharChandresekhar (1943), and Reference FellerFeller (1951), Reference Feller(1952). These processes are described by the following SDEs:

dyˆt=χdt+εdWˆt,   yˆt=y,(5.29)
dyˆt=χκyˆtdt+εdWˆt,   yˆt=y,(5.30)
dyˆt=χκyˆtdt+εyˆtdWˆt,   yˆt=y,(5.31)

respectively. It is clear that the corresponding L0,L1 are:

L0u=12ε2uyy+χuy,   L1u=0,(5.32)
L0u=12ε2uyy+χuy,   L1u=κuy,(5.33)
L0u=χuy,   L1u=12ε2uyyκuy.(5.34)

There are important differences among these processes. For an arithmetic Brownian motion, the operator L0 is a second-order differential operator, while L1 is zero, and the process is defined on the whole axis. For an OU process the operator L0 is a second-order differential operator, while L1 is a first-order operator; accordingly, this process is defined on the entire axis. In contrast, for a Feller process L0 is a first-order differential operator, while L1 is a second-order operator; hence, the process is only defined on a positive semiaxis.Footnote 5

5.3 Augmentation Procedure

While covering a lot of useful applications, OU and Feller processes are not sufficient to study all the practically important problems. Hence, one needs to enrich them via the so-called augmentation procedure; see Reference LiptonLipton (2001) . The underlying idea is straightforward. Given a stochastic process, say, an arithmetic Brownian motion, or an OU process, one can expand it by introducing additional stochastic variables depending on the original process. For example, an augmented Brownian motion (5.29) becomes a one-dimensional Kolmogorov process:

dxˆt=yˆtdt,   xˆt=x,dyˆt=χdt+εdWˆt,   yˆt=y.(5.35)

Similarly, one can augment OU and Feller processes as follows:

dxˆt=yˆtdt,   xt=x,dyˆt=χκyˆtdt+εdWˆt,   yt=y,(5.36)
dxˆt=yˆtdt,   xt=x,dyˆt=χκyˆtdt+εyˆtdWˆt,   yt=y,(5.37)

respectively. Of course, many other possibilities are practically important. In what follows, the Element analyzes several practically relevant and mathematically interesting augmented stochastic processes.

5.4 Reduction Procedure

Stochastic processes, which are not inherently affine, can often be transformed into an affine form through appropriate modifications. While some transformations are readily apparent, others demand significant effort and inspiration to identify, as highlighted by Reference Carr, Lipton and MadanCarr et al. (2002) and referenced works.

Consider the geometric Brownian motion, the cornerstone of mathematical finance and other disciplines. The associated stochastic process is not affine and is described by

dXˆt=μtXˆtdt+νtXˆtdWˆt,    Xˆt=X.(5.38)

Applying a logarithmic transformation,

Xˆtxˆt=lnXˆt,(5.39)

converts it into an arithmetic Brownian motion, which is affine:

dxˆt=μt12ν2tdt+νtdWˆt,    xˆt=x=lnX.(5.40)

This example illustrates that, with some ingenuity, even nonaffine processes like the geometric Brownian motion can be adapted for use with the existing analytical frameworks.

Another helpful example is transforming the Rayleigh process into the Feller process. Recall that the Rayleigh process describes a stochastic process on the positive semiaxis. We write this process as follows:

dσˆt=AσˆtBσˆtdt+CdWˆt,   σˆt=σ,(5.41)

where A,B,C>0. Define vˆt=σˆt2; then, according to Ito’s lemma, the dynamics of the process vˆt have the following form:

dvˆt=2A+C22Bvˆtdt+2CvˆtdWˆt,   vˆt=v=σ2.(5.42)

In financial applications considered in Section 8, the pair σ, v represents the volatility and variance of a price process.

6 Gaussian Stochastic Processes

6.1 Regular Gaussian Processes

Consider the governing system of SDEs, which might or might not be degenerate, and write the governing system of SDEs as follows:

dzˆt=b+Bzˆtdt+ΣdWˆt,   zˆt=z,(6.1)

where zˆt,b are M×1 vectors, and B and Σ are M×M matrices. Below, it is assumed that the corresponding coefficients are time-dependent.

The Fokker–Plank equation has the following form:

ϖtˉt,z,tˉ,zˉAϖzˉzˉt,z,tˉ,zˉ+b+Bzˉϖzˉt,z,tˉ,zˉ+bϖt,z,tˉ,zˉ=0,ϖt,z,t,zˉ=δzˉz,(6.2)

where, in agreement with the general (5.13), A is proportional to the covariance matrix,

A=amm=12ΣΣ*=12σmkσmk,   b=TrB=bmm.(6.3)

Recall that Einstein’s summation rule is used throughout the Element. Explicitly,

tˉϖammzˉmzˉmϖ+bm+bmmzˉmzˉmϖ+bϖ=0,ϖt,z,t,zˉ=δzˉz.(6.4)

The general Kolmogorov-type SDE, solvable via the Kelvin (or affine) ansatz, can be written in the following form:

dxˆt=bx+Bxxxˆt+Bxyyˆtdt,   xˆt=x,dyˆt=by+Byxxˆt+Byyyˆtdt+ΣyydWˆty,   yt=y,(6.5)

where xˆ t and bx are K×1 column vectors, yˆt and by are L×1 column vectors, Bxx, Bxy, Byx, Byy, and Σyy are K×K, K×L, L×K, L×L, and L×L matrices, respectively. In what follows, it is assumed that the corresponding coefficients are time-dependent. As usual, Wˆt is a standard L-dimensional Brownian motion.

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

dzˆt=bz+Bzzzˆtdt+0ΣyydWˆty,zˆt=xy,(6.6)

where

zˆt=xˆtyˆt,   bz=bxby,   Bzz=BxxBxyByxByy,(6.7)

so that zˆt and bz are M×1 column vectors, and Bzz is a M×M matrix, with M=K+L. In addition, define a scalar bz=TrBzz=TrBxx+TrByx.

The corresponding Fokker–Plank problem has the following form:

ϖtˉ(t,z,tˉ,zˉ)Aϖyy(t,z,tˉ,zˉ)+(b(z)+B(z)zˉ)ϖzˉ(t,z,tˉ,zˉ)+b(z)ϖ(t,z,tˉ,zˉ)=0,ϖ(t,z,t,zˉ)=δ(xx)δ(yy),(6.8)

where A has the following form:

A=all=12σllˉσllˉ=12ΣyyΣyy*.(6.9)

Explicitly,

tˉϖallzˉK+lzˉK+lϖ+bmz+bmmzzzˉmzˉmϖ+bzϖ=0,(6.10)
6.1.1 Solution via Kelvin Waves

By using the Kelvin-inspired ansatz, one can represent ϖ in the following form:

ϖt,z,tˉ,zˉ=12πMKt,z,tˉ,zˉ,mdm,Kt,z,tˉ,zˉ,m=expΨt,z,tˉ,zˉ,m,Ψt,z,tˉ,zˉ,m=αt,tˉ+iδt,tˉzˉimz,(6.11)

where m is an M×1 column vector, δ is an M×1 column vector, and

αt,t=0,   δt,t=m.(6.12)

Accordingly:

KtˉK=Ψtˉ=αtˉt,tˉ+iδtˉt,tˉzˉ,KzˉK=Ψzˉ=iδt,tˉ,   KzˉzˉK=Ψzˉ2=δt,tˉδ*t,tˉ.(6.13)

The coupled equations for α, δ have the following form:

αtˉt,tˉ+iδtˉt,tˉzˉ+δt,tˉAδt,tˉ+iδt,tˉb+Bzˉ+b=0,(6.14)

so that

αtˉt,tˉ+δt,tˉAδt,tˉ+iδt,tˉb+b=0,   αt,t=0,(6.15)
δtˉt,tˉ+B*δt,tˉ=0,   δt,t=m.(6.16)

Let Lt,tˉ be the fundamental solution of the homogeneous system of ODEs (6.16), namely, the matrix such that

tˉLt,tˉ+B*tˉLt,tˉ=0,   Lt,t=I.(6.17)

The solution of (6.16) has the following form:

δt,tˉ=Lt,tˉm.(6.18)

Thus,

αt,tˉ=12mC1t,tˉmimdt,tˉςt,tˉ,(6.19)

where C1 is an M×M positive-definite matrix of the following form:

C1t,tˉ=2ttˉL*t,sAsLt,sds,(6.20)

while d is an M×1 column vector,

dt,tˉ=ttˉL*t,sbsds,(6.21)

and ς is a scalar,

ςt,tˉ=ttˉbsds.(6.22)

Accordingly,

Ψt,tˉ,zˉ,m=12mC1t,tˉm+imL*t,tˉzˉdt,tˉzςt,tˉ.(6.23)

Thus,

ϖt,z,tˉ,zˉ=detCt,tˉ1/2expςt,tˉ2πM/2×Gt,tˉ,mexpimL*t,tˉzˉdt,tˉzdm,(6.24)

where Gt,tˉ,m is the density of a multivariate Gaussian distribution in the m-space. It is clear that ϖt,z,tˉ,zˉ is proportional to the characteristic function of G evaluated at the point L*t,tˉzˉdt,tˉz, so that

ϖt,z,tˉ,zˉ=detCt,tˉ1/2expςt,tˉ2πM/2×exp12L*t,tˉzˉdt,tˉzCt,tˉL*t,tˉzˉdt,tˉz.(6.25)

Thus, ϖ can be represented in the form:

ϖt,z,tˉ,zˉ=Nrt,tˉ,Ht,tˉ,(6.26)

where

Ht,tˉ=L*t,tˉ1C1t,tˉL1t,tˉ,rt,tˉ=L*t,tˉ1dt,tˉ+z.(6.27)

These results are applicable to the general Kolmogorov-type SDE solvable via the Kelvin (or affine) ansatz, which have the form (6.5). By using the same Kelvin ansatz as before, one can represent ϖ in the form (6.11):

ϖt,z,tˉ,zˉ=12πMKt,z,tˉ,zˉ,mdm,Kt,z,tˉ,zˉ,m=expΨt,z,tˉ,zˉ,m,Ψt,z,tˉ,zˉ,m=αt,tˉ+iδt,tˉzˉimz,(6.28)

where m is an M×1 column vector, m=k,l*, k is a K×1 column vector, l  is an L×1 column vector, δ is an M×1 column vector,  δ=β,γ*, β  is a K×1 column vector, γ  is an L×1 column vector, and

αt,t=0,   δt,t=βt,t,γt,t*=m=k,l*.(6.29)

As before:

KtK=Ψt=αtˉt,tˉ+iδtˉt,tˉzˉ,   KxK=Ψx=iβt,tˉ,KyK=Ψy=iγt,tˉ,   KyyK=Ψy2=γt,tˉγ*t,tˉ.(6.30)

The equations for α,δ have the following form:

αtˉt,tˉ+iδtˉt,tˉzˉ+γt,tˉAγt,tˉ+iδt,tˉbz+Bzzzˉ+bz=0.(6.31)

Accordingly,

αtˉt,tˉ+γt,tˉAγt,tˉ+iδt,tˉbz+bz=0,   αt,t=0,(6.32)
δtˉt,tˉ+Bzz*δt,tˉ=0,   δt,t=m=k,l*.(6.33)

Let Lt,tˉ be the fundamental solution of the homogeneous system of ODEs (6.33), namely, the matrix such that

tˉLt,tˉ+Bzz*tˉLt,tˉ=0,   Lt,t=I,(6.34)

where I is the identity matrix. The well-known Liouville’s formula yields

detLt,tˉ=expttˉbzsds.(6.35)

The solution of (6.32) is

δt,tˉ=Lt,tˉm.(6.36)

It is convenient to write Lt,tˉ in the block form:

Lt,tˉ=Lxxt,tˉLxyt,tˉLyxt,tˉLyyt,tˉ.(6.37)

It follows from (6.33) that

αt,tˉ=12mC1t,tˉmimdzt,tˉςt,tˉ,(6.38)

where C1 is an M×M positive-definite matrix split into four blocks of the form:

C1t,tˉ=2ttˉLyx*t,sAsLyxt,sdsttˉLyx*t,sAsLyyt,sdsttˉLyy*t,sAsLyxt,sdsttˉLyy*t,sAsLyyt,sds,(6.39)

while dz=dx,dy*, dx and dy are M×1 and N×1 column vectors, and ς is a scalar:

dzt,tˉ=ttˉL*t,sbzsds,(6.40)
ςt,tˉ=ttˉbzsds.(6.41)

Accordingly,

Ψt,z,tˉ,zˉ,m=12mC1t,tˉm+iLt,tˉmzˉimdzt,tˉ+zςt,tˉ=12mC1t,tˉm+imL*t,tˉzˉdzt,tˉzςt,tˉ.(6.42)

Thus,

ϖt,z,tˉ,zˉ=detC1/2expςt,tˉ2πM/2Gt,tˉ,m×expimL*t,tˉzˉdzt,tˉzdm,(6.43)

where Gt,tˉ,m is the density of a multivariate Gaussian distribution in the m-space. It is clear that ϖt,z,tˉ,zˉ is proportional to the characteristic function of G evaluated at the point L*t,tˉzˉdzt,tˉz, so that

ϖt,z,tˉ,zˉ=detC1/2expςt,tˉ2πM/2×exp12L*t,tˉzˉdzt,tˉzCL*t,tˉzˉdzt,tˉz.(6.44)

By using (6.35), one can rewrite (6.44) in the standard Gaussian form:

ϖt,z,tˉ,zˉ=Nrt,tˉ,Ht,tˉ,(6.45)

where the covariance matrix H and the mean r are as follows:

Ht,tˉ=L*t,tˉ1C1t,tˉL1t,tˉ,rt,tˉ=L*t,tˉ1dzt,tˉ+z.(6.46)
6.1.2 Solution via Coordinate Transform

Consider the Fokker–Planck problem (6.4). Introduce new variables:

tˉ,zˉtˉ,z˜=tˉ,Rtˉzˉ,   z˜m=rmmtˉzˉm,   rmm0=δmm.(6.47)

Then

tˉ=tˉ+tˉrmmzˉmz˜m,   zˉm=rmmz˜m.(6.48)

The transformed Fokker–Planck problem becomes

tˉϖ˜ammrnmrnmz˜nz˜nϖ˜+bmmzˉm+bmrnm+tˉrnmzˉmz˜nϖ˜+bϖ˜=0,ϖ˜t,z,t,z˜=δz˜z.(6.49)

To simplify the drift term, it is required that

tˉrmmt,tˉ+bnmt,tˉrmnt,tˉ=0,   rmmt,t=δnm.(6.50)

In matrix notation:

tˉRt,tˉ+Rt,tˉBt=0,    Rt,t=I.(6.51)

Thus, R=L*, rmm=lmm, where L is given by (6.34). It is easy to see that ϖ˜ satisfies the Fokker–Planck problem of the following form:

tˉϖ˜a˜nnt,tˉz˜nz˜nϖ˜+b˜nt,tˉz˜nϖ˜+bt,tˉϖ˜=0,ϖ˜t,z,t,z˜=δz˜z,(6.52)

with

a˜nnt,tˉ=lmnt,tˉammt,tˉlmnt,tˉ,b˜ntˉ=lmnt,tˉbmtˉ.(6.53)

In matrix notation:

A˜=L*t,tˉAt,tˉLt,tˉ,   b˜=L*t,tˉb.(6.54)

Accordingly,

ϖ˜t,z,tˉ,z˜=expttˉbsdsNz˜z+ttˉb˜nsds,ttˉC˜sds.(6.55)

Reverting back to the original variables, tˉ,z˜tˉ,zˉ, one recovers (6.45), as expected.

6.2 Killed Gaussian Processes

Consider a process governed by a system of SDEs (6.1), which is killed with intensity cˉ linearly depending on zˉ, namely,

cˉ=c+czˉ,(6.56)

where c is a scalar, and cz is an M×1 column vector. Thus, cˉ is the intensity at which the process goes into a “killed” state at some random time. The Fokker–Planck equation for a killed process has the following form:

ϖt¯(t,z,t¯,z¯)Aϖz¯z¯(t,z,t%¯,z¯)+(b+Bz¯)ϖ%z¯(t,z,t¯,z¯)+(b+c+cz¯)ϖ(t,z,t¯,%z¯)=0,ϖ(t,z¯,t,z)=δ(%z¯z).(6.57)

Explicitly,

ϖtˉammϖzˉmzˉm+bm+bmmzˉmϖzˉm+b+c+cmzˉmϖ=0,ϖt,zˉ,t,z=δzˉz.(6.58)

This problem can be solved by the same technique as before.

6.2.1 Solution via Kelvin Waves

The familiar Kelvin ansatz yields

αtˉt,tˉ+δt,tˉAδt,tˉ+iδt,tˉb+b+c=0,   αt,t=0,(6.59)
δtˉt,tˉ+B*δt,tˉic=0,   δt,t=m.(6.60)

Let Lt,tˉ be the fundamental solution of the homogeneous system of ODEs (6.60), namely, the matrix such that

tˉLt,tˉ+B*tˉLt,tˉ=0,   Lt,t=I,(6.61)

The solution of (6.60) has the following form:

δt,tˉ=Lt,tˉm+iLt,tˉttˉL1t,scsdsLt,tˉm+iet,tˉ,et,tˉ=ttˉL1t,scsds.(6.62)

Thus,

α=12mC1mimdς,(6.63)

where C1 is an M×M positive-definite matrix of the form:

C1t,tˉ=2ttˉL*t,sAsLt,sds,(6.64)

while d is an M×1 column vector,

dt,tˉ=ttˉL*t,sbs+AsLt,sesds,(6.65)

and ς=ς0+ς1 is a scalar,

ς0t,tˉ=ttˉbsds,ς1t,tˉ=ttˉcs12et,sL*t,sAsLt,ses12et,sL*t,sbsds.(6.66)

Accordingly,

Ψt,tˉ,zˉ,m=12mC1t,tˉm+imL*t,tˉzˉdt,tˉzLt,tˉet,tˉzˉςt,tˉ.(6.67)

Thus,

ϖt,z,tˉ,zˉ=detCt,tˉ1/2expLt,tˉet,tˉzˉς0t,tˉς1t,tˉ2πM/2×Gt,tˉ,mexpimL*t,tˉzˉdt,tˉzdm,(6.68)

where Gt,tˉ,m is the density of a multivariate Gaussian distribution in the m-space. It is clear that ϖt,z,tˉ,zˉ is proportional to the characteristic function of G evaluated at the point L*t,tˉzˉdt,tˉz, so that

ϖt,z,tˉ,zˉ=detCt,tˉ1/2expLt,tˉet,tˉzˉς0t,tˉς1t,tˉ2πM/2×exp12L*t,tˉzˉdt,tˉzCt,tˉL*t,tˉzˉdt,tˉz.(6.69)

It is often convenient to rewrite (6.69) as follows:

ϖt,z,tˉ,zˉ=Qt,tˉ,zˉNqt,tˉ,Ht,tˉ,(6.70)

where

Ht,tˉ=L*t,tˉ1C1t,tˉL1t,tˉ,qt,tˉ=L*t,tˉ1dt,tˉ+z,Qt,tˉ,zˉ=expLt,tˉet,tˉzˉς1t,tˉ.(6.71)

As could be expected, the probability ϖ is no longer conserved due to a prefactor Q, reflecting the fact that the process is killed with intensity cˉ.

It is worth noting that Q depends on zˉ but does not depend on z. Completing the square, one can represent ϖ in the form:

ϖt,z,tˉ,zˉ=Rt,z,tˉNrt,tˉ,Ht,tˉ,(6.72)

where

rt,tˉ=L*t,tˉ1dt,tˉ+zC1t,tˉet,tˉ=qt,tˉHt,tˉLt,tˉet,tˉ,Rt,z,tˉ=expet,tˉdt,tˉ+z+12et,tˉC1t,tˉet,tˉς1t,tˉ.(6.73)

It is clear that R depends on z but does not depend on zˉ. Accordingly, (6.72) is easier to use than (6.70) when future expectations are calculated.

The same formulas can be derived via the method of coordinate transforms. Details are left to the interested reader.

6.3 Example: Kolmogorov Process

Extend the Kolmogorov formula to the case when b and σ are functions of time, bt and σt. The corresponding SDE has the following form:

dxˆt=yˆtdt,   xˆt=x,dyˆt=btdt+σtdWˆt,   yˆt=y.(6.74)

Accordingly, (6.34) can be written as follows:

Lt,tˉ+0010Lt,tˉ=0,  Lt,t=1001,(6.75)

so that

Lt,tˉ=10T1,   L1t,tˉ=10T1.(6.76)

Once Lt,tˉ is known, one can compute C1t,tˉ, dzt,tˉ, ςt,tˉ:

C1t,tˉ=ψ2t,tˉψ1t,tˉψ1t,tˉψ0t,tˉ,dzt,tˉ=dxt,tˉdyt,tˉ=ϕ1t,tˉϕ0t,tˉ,ςt,tˉ=0,(6.77)

where

ϕit,tˉ=ttˉstibsds,     ψit,tˉ=ttˉstiσ2sds.(6.78)

Next, the covariance matrix Ht,tˉ, and the mean rt,tˉ are calculated as follows:

Ht,tˉ=L*t,tˉ1C1t,tˉL1t,tˉ=1T01ψ2t,tˉψ1t,tˉψ1t,tˉψ0t,tˉ10T1(6.79)
=ψ0t,tˉT22ψ1t,tˉT+ψ2t,tˉψ0t,tˉTψ1t,tˉψ0t,tˉTψ1t,tˉψ0t,tˉ,rt,tˉ=ϕ1t,tˉ+x+ϕ0t,tˉ+yTϕ0t,tˉ+y.(6.80)

Accordingly, ϖt,x,y,tˉ,xˉ,yˉ is a bivariate Gaussian distribution of the form (6.26), with

σx2t,tˉ=ψ0t,tˉT22ψ1t,tˉT+ψ2t,tˉ,   σy2=ψ0t,tˉ,ρt,tˉ=ψ0t,tˉTψ1t,tˉψ0t,tˉψ0t,tˉT22ψ1t,tˉT+ψ2t,tˉ,pt,tˉ=ϕ1t,tˉ+x+ϕ0t,tˉ+yT,   qt,tˉ=ϕ0t,tˉ+y.(6.81)

It is left to the interested reader to verify that (6.81) coincides with (3.52) when σ and b are constant. Therefore, the classical Kolmogorov solution can be extended to the case of time-dependent parameters.

6.4 Example: OU Process

6.4.1 OU Process

It is worth deriving the well-known t.p.d.f. for the OU process using Kelvin waves for benchmarking purposes. The following SDE governs the OU process:

dyˆt=χtκtyˆtdt+εtdWˆt,   yˉt=y.(6.82)

Equivalently,

dyˆt=κtθtyˆtdt+εtdWˆt,   yˉt=y,(6.83)

where θt=χt/κt.

The corresponding Fokker–Planck problem has the following form:

ϖtˉt,y,tˉ,yˉ12ε2ϖyˉyˉt,y,tˉ,yˉ+χtˉκtˉyˉϖyˉt,y,tˉ,yˉκtˉϖt,y,tˉ,yˉ=0,ϖt,y,tˉ,yˉ=δyˉy.(6.84)

The associated function Kt,y,tˉ,yˉ,l has the following form:

K=expαt,tˉ+iγt,tˉyˉily,(6.85)

so that

αtˉt,tˉ+12ε2tˉγ2t,tˉ+iχtˉγt,tˉκtˉ=0,   αt,t=0,γtˉt,tˉκtˉγt,tˉ=0,   γt,t=l.(6.86)

Thus,

γt,tˉ=eηt,tˉl,αt,tˉ=12ψ0t,tˉl2ttˉeηt,sχsdsil+ηt,tˉ.(6.87)

where

ηt,tˉ=ttˉκsds,(6.88)

ψ0(t,t¯)=tt¯e2η(t,s)ε2(s)ds.(6.89)

Since the same quantities will appear regularly throughout the Element, it is convenient to introduce the following notation:

Aκt,tˉ=eηt,tˉ,   Bκt,tˉ=ttˉeηt,sds,   Bˉκt,tˉ=ttˉeηs,tˉds,Aκt,tˉ=eηt,tˉ,   Bκt,tˉ=ttˉeηt,sds,   Bˉκt,tˉ=ttˉeηs,tˉds,(6.90)

In particular, for constant κ, one has

Aκt,tˉ=eκT=AκT,   Aκt,tˉ=eκT=AκT,Bκt,tˉ=Bˉκt,tˉ=1eκTκ=BκT=BˉκT,Bκt,tˉ=Bˉκt,tˉ=eκT1κ=BκT=BˉκT,(6.91)

and

A0t,tˉ=1,   B0t,tˉ=Bˉ0t,tˉ=T.(6.92)

In this notation, ψ0 can be written as follows:

ψ0t,tˉ=ttˉA2κt,sε2sds.(6.93)

Thus, the following well-known expression is obtained:

ϖt,y,tˉ,yˉ=12πexpψ0t,tˉl22+eηt,tˉyˉttˉeηt,sχsdsyil+ηt,tˉdl=Aκt,tˉ2πψ0t,tˉexpAκt,tˉyˉttˉAκt,sχsdsy22ψ0t,tˉ=12πψˆ0t,tˉexpyˉttˉAκt,sχsdsAκt,tˉy22ψˆ0t,tˉ,(6.94)

where

ψˆ0t,tˉ=A2κt,tˉψ0t,tˉ=ttˉA2κs,tˉε2sds.(6.95)

For further discussion, see the original paper by Reference Uhlenbeck and OrnsteinUhlenbeck and Ornstein (1930), as well as Reference ChandrasekharChandresekhar (1943), Reference RiskenRisken (1989), and references therein.

For time-independent parameters, (6.94) has the form:

ϖt,y,tˉ,yˉ=12πΣ2t,tˉexpyˉθAκTyθ22Σ2t,tˉ,(6.96)

with

Σ2t,tˉ=ε21e2κT2κ=ε2B2κT.(6.97)
6.4.2 Gaussian Augmented OU Process

This subsection considers an augmented one-dimensional OU process of the form:

dxˆt=yˆtdt,   xˆt=x,dyˆt=χtκtyˆtdt+εtdWˆt,   yˆt=y.(6.98)

To align the analysis with the existing body of work, switch from the general notation, used above, to a specific one customarily used for the OU process. Here and in what follows, the word “augmentation” means that one expands the original process by incorporating its integral or other path-dependent characteristics, such as running maximum or minimum as part of the process; see Section 5. The augmentation is a very useful tool. In particular, in financial engineering it is used for handling large classes of path-dependent options; details can be found in Reference LiptonLipton (2001), chapter 13.

For an OU process, (6.34) can be written as follows:

Ltˉt,tˉ+001κtˉLt,tˉ=0,  Lt,t=1001,(6.99)

so that

Lt,tˉ=10Bˉκt,tˉAκt,tˉ,   L1t,tˉ=2pt10Bκt,tˉAκt,tˉ.(6.100)

Now, one can compute C1t,tˉ, dzt,tˉ, and ςt,tˉ:

C1t,tˉ=ψ2t,tˉψ1t,tˉψ1t,tˉψ0t,tˉ,(6.101)

where

ψ0t,tˉ=ttˉAκ2t,sε2sds,ψ1t,tˉ=ttˉBˉκt,sAκt,sε2sds,(6.102)
ψ2t,tˉ=ttˉBˉκ2t,sε2sds.dzt,tˉ=dxt,tˉdyt,tˉ=ttˉBˉκt,sχsdsttˉAκt,sχsds,,(6.103)
ςt,tˉ=ηt,tˉ.(6.104)

Next, one can calculate the covariance matrix Ht,tˉ, and mean vector rt,tˉ as follows:

Ht,tˉ=L*t,tˉ1C1t,tˉL1t,tˉ=1 Bκt,tˉ0 Aκt,tˉψ2t,tˉψ1t,tˉψ1t,tˉψ0t,tˉ10 Bκt,tˉ Aκt,tˉ=h0t,tˉh1t,tˉh1t,tˉh2t,tˉ,(6.105)
rt,tˉ=L*t,tˉ1dzt,tˉ+xy=pt,tˉqt,tˉ.(6.106)

Here

h0t,tˉ=ψ0 Bκ2t,tˉ2ψ1 Bκt,tˉ+ψ2,h1t,tˉ=ψ0 Bκt,tˉψ1 Aκt,tˉ,h2t,tˉ=ψ0Aκ2t,tˉ,(6.107)
pt,tˉ=ttˉBˉκt,sχsds+x+Bκt,tˉttˉAκt,sχsds+y,qt,tˉ=Aκt,tˉttˉAκt,sχsds+y.(6.108)

Thus, ϖt,x,y,tˉ,xˉ,yˉ is a bivariate Gaussian distribution of the form (6.26) with the covariance matrix H, given by (6.105) centered at the point r=p,q* given by (6.106). Explicitly, one has

σx2t,tˉ=h0t,tˉ,   σy2t,tˉ=h2t,tˉ,   ρt,tˉ=h1t,tˉh0t,tˉh2t,tˉ.(6.109)

When χ,κ,θ,ε are constant, the preceding formulas become significantly simpler. Namely,

LT=10BκTAκT,   L1T=10BκTAκT,(6.110)
C1T=ε2κ2B0T2BκT+B2κTε22BκTε22BκTε2B2κT,(6.111)
dzT=TBκTθBκTχ,(6.112)
ςT=κT,(6.113)
HT=ε2κ2B0T2BκT+B2κTε22BκTε22BκTε2B2κT,(6.114)
rT=x+θTBκTθyθAκTθy.(6.115)

Thus, when coefficients are constant, ϖt,x,y,tˉ,xˉ,yˉ is a bivariate Gaussian distribution of the form (6.26) with the covariance matrix H, given by (6.114) and the mean vector r=p,q* given by (6.115).

Calculate the marginal distribution of xˉ, denoted by ϖxt,y,tˉ,xˉ, which is used on several occasions in what follows. It is well known that marginal distributions of a multivariate Gaussian distribution are also Gaussian, so that

ϖxt,y,tˉ,xˉ=12πh0t,tˉexpxˉpt,tˉ22h0t,tˉ,(6.116)

where h0 is given by the equations in (6.114). At the same time, the density of marginal distribution for yˉ has the form

ϖyt,y,tˉ,yˉ=12πh2t,tˉexpyˉqt,tˉ22h2t,tˉ,(6.117)

where h2 is given by the equations in (6.114), which is the familiar density of the OU process derived in the previous section.

6.5 Example: Diffusion of Free and Harmonically Bound Particles

The preceding results can be used to revisit the motion of free and harmonically bound particles considered in Section 3.

To describe a free particle, it is assumed that χ=0. Equation (6.114) does not change, while (6.115) can be simplified as follows:

pTqT=x+BκTyAκTy.(6.118)

It is clear that Equations (4.7), (4.8) and (6.114), (6.118) are in agreement. A typical free particle behavior is illustrated in Figure 7.

Figure 7 A thousand trajectories of a typical free particle. Parameters are as follows: T=5, dt=0.01, κ=0.8, σ=1.0. (a) xt, (b) yt, (c) xˉT,yˉT, (d) contour lines of ϖ0,0,0,T,x˜,y˜. Author’s graphics.

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

Lt,tˉ+0ω21κLt,tˉ=0,  Lt,t=1001.(6.119)

The corresponding characteristic equation and its solutions are as follows:

λ2κλ+ω2=0,(6.120)
λ±=μ±ζ,μ=κ2,   ζ=κ24ω22.(6.121)

Introduce

E0T=eμT=eκT/2,   E±T=e±ζT.(6.122)

It is left to the reader to check that

L=E0κ24ω2λE+λ+Eω2E+EE+Eλ+E+λE,(6.123)
L1=E01κ24ω2λ+E+λEω2E+EE+EλE+λ+E,(6.124)
detL=detL11=E02=eκT,(6.125)

Accordingly,

ϖt,z,tˉ,zˉ=Nrt,tˉ,Ht,tˉ,(6.126)

with

H=L*1C1L1,r=L*1z.(6.127)

Here,

C1=ψ2ψ1ψ1ψ0,(6.128)

where

ψ0=ε2κ24ω2ttˉλ+eλ+stλeλst2ds=ε22κκ24ω2E02κλ+E+24ω2+κλE2κ24ω2,ψ1=ε2κ24ω2ttˉeλ+steλstλ+eλ+stλeλstds=ε22κ24ω2E02E+E2,ψ2=ε2κ24ω2ttˉeλ+steλst2ds=ε22κω2κ24ω2E02κλE+24ω2+κλ+E2κ24ω2.(6.129)

Further,

H=E02κ24ω2λ+E+λEE+Eω2E+EλE+λ+E×ψ2ψ1ψ1ψ0λ+E+λEω2E+EE+EλE+λ+E.(6.130)

Straightforward but tedious calculation yields

h0=ε22κω21E02ω2E+E2+λ+E+λE2κ24ω2,h1=ε22E02E+E2κ24ω2,h2=ε22κ1E02ω2E+E2+λE+λ+E2κ24ω2.(6.131)

In the limit ω20,

h0=ε2κ2B02Bκ+B2κ,   h1=ε22Bκ2,   h2=ε2B2κ,(6.132)

so that Equations (6.114) and (6.132) are in agreement.

Here

r=pq=E01λ+E+λEx+E+Eyκ24ω2E01ω2E+Ex+λE+λ+Eyκ24ω2.(6.133)

In the limit ω20,

r=pq=x+BκTyAκTy.(6.134)

Moreover, while it is easy to show that Chandrasekhar’s solution given in Reference ChandrasekharChandresekhar (1943) is in agreement with the solution given by (6.126), the solution is more convenient from a practical standpoint, since it is explicitly written as a Gaussian density in the xˉ,yˉ space. A typical bounded particle behavior is shown in Figure 8.

Figure 8 A thousand trajectories of a harmonically bounded particle. Parameters are as follows: T=5, dt=0.01, κ=0.2, ω=0.5, σ=0.5. (a) xt, (b) yt, (c) xˉT,yˉT, (d) contour lines of ϖ0,0,0,T,x˜,y˜. Author’s graphics.

6.6 Example: Vorticity of Two-Dimensional Flows

Briefly return to the starting point and consider strictly two-dimensional flows; see Reference Friedlander and Lipton-LifschitzFriedlander and Lipton-Lifschitz (2003). Velocity fields of such flows have the following form:

Vtˉ,xˉ1,xˉ2=V1tˉ,xˉ1,xˉ2,V2tˉ,xˉ1,xˉ2,vtˉ,xˉ1,xˉ2=v1tˉ,xˉ1,xˉ2,v2tˉ,xˉ1,xˉ2.(6.135)

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

V1=Ψxˉ2,   V2=Ψxˉ1,   v1=ψxˉ2,   v2=ψxˉ1,(6.136)

and define the scalar vorticity as follows:

Ω=ΔΨ,   ω=Δψ.(6.137)

Contour lines of Ψ are called streamlines of the flow.

By using the preceding definitions, the two-dimensional Navier–Stokes equations can be written as equations for the stream and vorticity:

ΩtˉΨxˉ2Ωxˉ1+Ψxˉ1Ωxˉ2νΔΩ=0,ΔΨΩ=0.(6.138)

Time-independent quadratic stream functions Ψxˉ1,xˉ2 generate exact equilibrium solutions of the equations in (6.138). Consider fields consisting of pure strain and pure rotation. The corresponding Ψ have the following form:

Ψxˉ1,xˉ2=14wxˉ12+xˉ222sxˉ1xˉ2,(6.139)

where w>s, to ensure that streamlines are elliptic rather than hyperbolic, so that

V1=Ψxˉ2=12sxˉ1wxˉ2,   V2=Ψxˉ1=12wxˉ1sxˉ2.(6.140)

Recall that these flows were introduced in Section 2, Equation (2.7).

Small perturbations ψ of the time-independent quadratic stream function Ψ satisfy the following equations:

ωtˉΨxˉ2ωxˉ1+Ψxˉ1ωxˉ2νΔω=0,Δψω=0.(6.141)

It is helpful to study the first equation (6.141) in isolation, by writing it explicitly as follows:

ωtˉ+12sxˉ1wxˉ2ωxˉ1+12wxˉ1sxˉ2ωxˉ2νΔω=0,(6.142)

and supplying it with the initial condition at time t:

ωt,xˉ1,xˉ2=δxˉ1x1δxˉ2x2.(6.143)

Once the solution of Equations (6.142) and (6.143) is found, one can find ψ by solving the corresponding Laplace equation.

Surprisingly, this equation is identical to the Fokker–Planck equation associated with the following SDEs for zˆt=xˆ1t,xˆ2t:

dzˆt=Bzˆtdt+ΣdWˆt,   zˆt=x1x2,(6.144)

where

B=12swws,   Σ=2ν1001.(6.145)

Thus, one can use Section 6.1 results. Equation (6.34) becomes

Lt,tˉ+12swwsLt,tˉ=0,  Lt,t=1001.(6.146)

The corresponding characteristic equation has the following form:

λ2+14w2s2=0.(6.147)

Its solutions are

λ±=±ζ,   ζ=iw2s22.(6.148)

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

L=c1s2ζs1w2ζs1w2ζs1c1+s2ζs1,   detL=1L1=c1+s2ζs1w2ζs1w2ζs1c1s2ζs1,   detL1=1,(6.149)

where

c1t,tˉ=cosζT,   s1t,tˉ=sinζT.(6.150)

Next, (6.39) yields

C1=2νttˉL*t,sLt,sds=ψ2ψ1ψ1ψ0,(6.151)

where

ψ0=2νttˉ1+s2ζs2t,s+s24ζ21c2t,sds=2ν1+s24ζ2Ts4ζ2c2s28ζ3s2,ψ1=νsw2ζ2ttˉ1c2t,sds=νsw2ζ2T12ζs2,ψ2=2νttˉ1s2ζs2t,s+s24ζ21c2t,sds=2ν1+s24ζ2T+s4ζ2c2s28ζ3s2,(6.152)

and

c2t,tˉ=cos2ζT,   s2t,tˉ=sin2ζT.(6.153)

Finally, Equations (6.26) and (6.27) yield:

ωt,z,tˉ,zˉ=Nrt,tˉ,Ht,tˉ.(6.154)

The corresponding covariance matrix H and mean r are as follows:

H=h0h1h1h2,(6.155)

where

h0=w28ζ2+4ζ2s28ζ2c2+s2ζs2ψ2+w2ζs2ζ1c2+s2ψ1+w28ζ21c2ψ0,h1=w4ζs2ζ1c2+s2ψ21+w24ζ21c2ψ1+w4ζs2ζ1c2s2ψ0,h2=w28ζ21c2ψ2+w2ζs2ζ1c2s2ψ1+w28ζ2+4ζ2s28ζ2c2s2ζs2ψ0,(6.156)

and

r=r1r2=c1+s2ζs1x1w2ζs1x2w2ζs1x1+c1s2ζs1x2.(6.157)

The equations in (6.156) are symmetric, namely h0h2 when a,ba,b and ψ0,ψ2ψ2,ψ0. The second of the equations in (6.138), which is a static Poisson equation, allows us to find ψ, since ω is known. Its analytical solution is not easy to derive and is not presented here due to lack of space. However, the special case of purely rotational flow, s=0, can be done easily; see (6.165).

It is interesting to note that

Ψr1,r2=Ψx1,x2,(6.158)

so that the location of the Gaussian distribution ω moves along streamlines of the flow defined by the stream function Ψ.

When the flow is purely rotational, so that s=0, the preceding formulas considerably simplify. Specifically, one has the following:

ψ0=2νT,  ψ1=0,     ψ2=2νT,h0=2νT,  ψ1=0,     h2=2νT,r1=c1x1s1x2,   r2=s1x1+c1x2,(6.159)

so that

ωt,x1,x2,tˉ,xˉ1,xˉ2=14πνTexpxˉ1c1x1+s1x22+xˉ2s1x1c1x224νT.(6.160)

The stream function ψ can be calculated directly by solving the corresponding Poisson equation.Footnote 6 To start, notice that both ω and ψ are rotational symmetric around the point x1,x2. Thus, ω and ψ have the following form:

ω=ωR=14πνTexpR22,   ψ=ψR,(6.161)

where

R2=xˉ1c1x1+s1x22+xˉ2s1x1c1x222νT.(6.162)

Then ψR solves a radially symmetric Poisson equation of the following form:

1RRψRRR=12πexpR22.(6.163)

Thus,

RψRR=12πexpR22+C,(6.164)

where C is an arbitrary constant. Next,

ψR=12πlnR+12E1R22,(6.165)

where the choice of C guarantees that ψ has the right behavior when R0 and R. Here E1η is the exponential integral of the following form:

E1η=ηeηηdη.(6.166)

7 Non-Gaussian Stochastic Processes

7.1 Regular Non-Gaussian Processes

In many situations, it is useful to consider processes governed by more general SDEs of the following form:

dzˆt=bt+Btzˆtdt+Σtdiagd0t+Dtzˆt1/2dWˆtz,zˆt=z.(7.1)

Here, in addition to the functions bt, Bt introduced in the previous section, define an M×1 column vector d0, and an M×M matrix D. It is convenient to introduce auxiliary vectors di equal to the i th column of D.

Since the corresponding M×M covariance matrix A has the form:

A=12Σdiagd0t+Dtzt1/2Σdiagd0t+Dtzt1/2*,(7.2)

it linearly depends on z:

A=12Σdiagd0+DzΣ*=12A0+12Amzm,(7.3)

where

A0=12Σdiagd0Σ*,   Ai=12ΣdiagdiΣ*.(7.4)

In contrast to the Gaussian case, the equations in (7.2) have to be defined in the domain D such that

D=zd0+Dz0,(7.5)

rather than in the whole space. In financial engineering, covariance matrices of the form (7.2) were introduced by Reference Dai and SingletonDai and Singleton (2000), and discussed by Reference Duffie, Filipovic and SchachermayerDuffie et al. (2003), Reference FilipovicFilipovic (2009), and many others.

The corresponding Fokker–Plank problem has the following form:

ϖtˉ(t,z,tˉ,zˉ)(A0+zˉmA(m))ϖzˉzˉ(t,z,tˉ,zˉ)+(bˆ+Bzˉ)ϖzˉ(t,z,tˉ,zˉ)+bϖ(t,z,tˉ,zˉ)=0,ϖ(t,z,t,zˉ)=δ(zˉz),(7.6)

where

bˆm=bm2ammm+ammm+ammm,  (no summation over m), b=TrB.(7.7)

Equation (6.11) expressing ϖ in terms of K holds. The equations for α,δ have the following form:

αtˉt,tˉ+iδtˉt,tˉzˉ+δt,tˉAδt,tˉ+iδt,tˉbˆ+Bzˉ+b=0,(7.8)

or, more explicitly,

αtˉt,tˉ+iδtˉt,tˉzˉ+δt,tˉA0δt,tˉ+δt,tˉAkδt,tˉzˉk+iδt,tˉbˆ+Bzˉ+b=0.(7.9)

Thus, the system of ODEs for α,δ can be written as follows:

αtˉt,tˉ+δt,tˉA0t,tˉδt,tˉ+iδt,tˉbˆt,tˉ+bt,tˉ=0,  αt,t=0,iδit,tˉ+δt,tˉAiδt,tˉ+iBijδjt,tˉ=0,   δit,t=mi.(7.10)

In the case in question, the equation for δ is no longer linear. Instead, δ 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; some of the results in this direction are reported here. However, in the one-dimensional case, the corresponding Riccati equation can be converted into the second-order ODE, and then solved explicitly when the coefficients A, b, b are time-independent.

In case of an augmented process, one must consider an SDE of the following form:

dxˆt=bxt+Bxxtxˆt+Bxytyˆtdt,dyˆt=byt+Byxtxˆt+Byytyˆtdt+Σyytdiagd0t+Dtzˆt1/2dWˆty,xˆt=x,   yˆt=y,(7.11)

or, more compactly,

dzˆt=bt+Btzˆtdt+0Σyytdiagd0t+Dtzˆt1/2dWˆty,zˆt=z=xy.(7.12)

Here bx, by, b=bx,by*, d0 are column vectors, and Bxx, Bxy, Byx, Byy, B, and D are matrices of appropriate dimensions.

The equations for α,δ=β,γ have the following form:

αtˉt,tˉ+iδtˉt,tˉzˉ+γt,tˉAγt,tˉ+iδt,tˉbˆ+Bzˉ+b=0,(7.13)

or, more explicitly,

αtˉt,tˉ+iδtˉt,tˉzˉ+γt,tˉA0γt,tˉ+γt,tˉAkγt,tˉzˉk+iδt,tˉbz+Bzzzˉ+b=0.(7.14)

Thus, the system of ODEs for α,δ can be written as follows:

αtˉt,tˉ+γt,tˉA0t,tˉγt,tˉ+iδt,tˉbˆzt,tˉ+bt,tˉ=0,αt,t=0,iδi,tˉt,tˉ+γt,tˉAiγt,tˉ+iBijzzδjt,tˉ=0,   δit,t=mi.(7.15)

7.2 Killed Non-Gaussian Processes

The non-Gaussian governing SDE has the following form:

dzˆt=b+Bzˆtdt+Σdiagd0+Dzˆt1/2dWˆt,(7.16)

where zˆt, b, d0 are M×1 vectors, and Σ, B, D are the M×M matrices defined previously. As before, the correlation matrix Σ can be a full-rank (nondegenerate) matrix. Once again, it is assumed that the process is killed with intensity cˉ linearly depending on z, namely,

cˉ=c+cz,(7.17)

where c is a scalar, and cz is an M×1 column vector.

The corresponding Fokker–Plank problem has the following form:

ϖtˉ(t,z,tˉ,zˉ)(A0+zˉiAi)ϖzˉzˉ(t,z,tˉ,zˉ)+(bˆ+Bzˉ)ϖzˉ(t,z,tˉ,zˉ)+(b+c+czˉ)ϖ(t,z,tˉ,zˉ)=0,ϖ(t,z,t,zˉ)=δ(zˉz).(7.18)

The equations for α,δ generalize the equations in (7.10). They can be written in the following form:

αtˉt,tˉ+δt,tˉA0t,tˉδt,tˉ+iδt,tˉbt,tˉ+bt,tˉ+ct,tˉ=0,αt,t=0,iδit,tˉ+δt,tˉAiδt,tˉ+iBijδjt,tˉ+ci=0,   δit,t=mi.(7.19)

As in the case without killing, finding an analytical solution to a multidimensional Riccati equation is generally impossible. However, in the time-independent one-dimensional case, it can be done. Solution becomes particularly simple in the special case when A0=0. The most important case is the killed one-dimensional Feller process, used, for example, to price bonds in the Cox–IngersolI–Ross (CIR) model; see Section 8.

7.3 Example: Anomalous Kolmogorov Process

Anomalous diffusion is a phenomenon in which the random motion of particles or molecules deviates from the classical Brownian motion and, as a result, exhibits non-Gaussian probability distributions, such as power-law or exponential tails. One can distinguish between subdiffusions (slower spreading) and superdiffusions (faster spreading). Anomalous diffusion often involves long-range correlations in particle motion, meaning that the movement of a particle at a one-time step depends on its previous positions over longer time scales. Anomalous diffusion frequently displays scale-invariant properties, meaning that the statistical properties of motion remain the same across different time or spatial scales. Anomalous diffusion has applications in physics, chemistry, financial engineering, biology, and geophysics.

Fractional Brownian motion (fBm) is used to model anomalous diffusion because it possesses several relevant characteristics. In particular, it exhibits long memory, which means that the process’s future values are influenced by its past values over long time scales. Additionally, fBm can produce non-Gaussian behavior while preserving scale-invariance. By adjusting the Hurst exponent and other parameters, fBm can be tailored to model different anomalous diffusions, including both subdiffusions and superdiffusions.

This section studies a fractional Kolmogorov equation of the following form:

ϖtˉt,x,y,tˉ,xˉ,yˉ+a2yˉ2νϖt,x,y,tˉ,xˉ,yˉ+yˉϖxˉt,x,y,tˉ,xˉ,yˉ+bϖyˉt,x,y,tˉ,xˉ,yˉ=0,ϖt,xˉ,yˉ,t,x,y=δxˉxδyˉy,(7.20)

where 0<ν<1. The pseudo-differential operator 2yˉ2ν is defined as follows:

2yˉ2νϖ=F1l2νFϖ.(7.21)

Here F and F1 denote the direct and inverse Fourier transforms, respectively. Despite its complexity, problem (7.20) can be solved by using Kelvin waves. For particular solutions of the form (3.38), (3.39), and (3.40), the corresponding characteristic equations are

αtˉt,tˉ+aγ2νt,tˉ+iγtˉt,tˉyˉ+ikyˉ+ibγt,tˉ=0,αt,t=0,   γt,t=l,(7.22)

so that

αtˉt,tˉ+aγ2νt,tˉ+ibγt,tˉ=0,   αt,t=0,γtˉt,tˉ+k=0,   γt,t=l,(7.23)
γt,tˉ=kT+l,αt,tˉ=attˉkst+l2νdsibkT22+lT.(7.24)

Thus,

Ψ=α+ikxˉx+iγyˉily=attˉkst+l2νds+ikxˉxyˉT+bT22+ilyˉybT.(7.25)

Now, assume that ν=1/2. The key is to calculate the integral

I=ttˉkst+l2νds=0Tks+lds,(7.26)

for different values of k,l. Depending on k,l, this integral can be calculated as follows:

I1=0l/kks+lds+l/kTkslds=kT22lT+l2k,   0k<,0lkT,I2=kT22+lT,   0k<,   kTl<,I3=kT22+lT,   <k0,   0l<,I4=0l/kkslds+l/kTks+lds=kT22+lTl2k,   <k0,kTl0,I5=kT22lT,  <k0,   <lkT,I6=kT22lT,  0k<,   <l0.(7.27)

Thus,

2π2J1=00kTexpakT22lT+l2k+ikaT2ζ+ilaTηdkdl=T010expp+iqkkdχdk=T01p0expp+iqkdkdχ=T01dχpiq2=1a2T30Tdχχf+χf2,(7.28)

where ζ,η are nondimensional variables:

ζ=xˉxyˉT+bT22aT2,   η=yˉybTaT,(7.29)
l=Tχk,   pχ=aT212χ+χ2>0,   qχ=aT2ζ+χη,(7.30)

and f± are roots of the quadratic equation

χ21+iηχ+12iζ=0.(7.31)

One can check that

f±=1+iη±1+iη22+4iζ2=1+iη±iD2,f+f=12iζ,   f++f=1+iη,   f+f=iD,(7.32)

with

D=1+η24iζ2iη.(7.33)

The roots f± are never equal, since D does not vanish when ζ,η are real.

Thus, one has

2π2J1=1a2T301dχχf+χf2=1a2T3f+f2011χf+1χf2dχ=1a2T3f+f211f++1f++11f+1f+2f+flnf1f+f+1f=1a2T3D4D+2iζ+iη12iζ12iζ2iη2iDln2ζ+ηD2ζ+η+D.(7.34)

By symmetry,

J4ζ,η=J1ζ,η=J1ζ,η.(7.35)

Next,

2π2J2=0kTexpakT22+lT+ikaT2ζ+ilaTηdkdl=1aT1iη0expkaT22+ikaT2ζ+ηdk=1a2T31iη12iζ+η.(7.36)

Similarly, it is easy to show that

2π2J3=00expakT22+lT+ikaT2ζ+ilaTηdkdl=1a2T31iη12+iζ,(7.37)

while, by symmetry, one gets

J5ζ,η=12π2a2T31+iη12+iζ+η=J2ζ,η,J6ζ,η=12π2a2T31+iη12iζ=J3ζ,η,(7.38)

so that

ϖ=1π2a2T311+η212ηζ+η1+4ζ+η2+1+2ηζ1+4ζ2+Re1D2D+2iζ+iηD2ζ+η2iDln2ζ+ηD2ζ+η+D,(7.39)

and

ϖxˉ,yˉdxˉdyˉ=1π2a211+η212ηζ+η1+4ζ+η2+1+2ηζ1+4ζ2+Re1D2D+2iζ+iηD2ζ+η2iDln2ζ+ηD2ζ+η+Ddζdηϖζ,ηdζdη,(7.40)

which shows that, as expected, in the nondimensional variables there is no explicit dependence on T.Footnote 7

A typical anomalous Kolmogorov process is depicted in Figure 9. The difference between the anomalous diffusion shown in Figure 9 and the pure diffusion shown in Figure 6 is clear.

Figure 9 Contour lines of ϖ0,0,0,T,x˜,y˜ for an anomalous Kolmogorov process with T=1.5, a=2.5, b=1.5 . Author’s graphics.

It is worth comparing Equations (7.39) and (3.28). To this end, rewrite Φ given by (3.29) in the following form:

Φ=yˉybT22aT+6xˉxyˉ+yT22aT3=η22+32ζ+η22,(7.41)

where ζ,η are nondimensional variables of the form:

ζ=xˉxyˉT+bT22aT3,   η=yˉybTaT, (7.42)

and a is the diffusion coefficient; its dimension is a=L2/T3. Thus,

ϖ=3πaT2expη2+32ζ+η22,(7.43)

and

ϖxˉ,yˉdxˉdyˉ=3πaexpη2+32ζ+η22dζdηϖζ,ηdζdη.(7.44)

Comparing Equations (7.39) and (7.43), one can see that the scaling of ϖ and its asymptotic behavior at infinity is completely different.

7.4 Example: Feller Process

7.4.1 Feller Process
Feller Process with Constant Parameters

For benchmarking purposes, it is useful to start with deriving the well-known t.p.d.f. for the Feller process with constant coefficients; see Reference FellerFeller (1951), Reference Feller(1952):

dyˆt=χκyˆtdt+εyˆtdZˆt,   yˆt=y.(7.45)

Initially, the process with time-independent parameters is considered; the time-dependent case is analyzed later in this section.

To start with, it is assumed that

2χε21ϑ>0.(7.46)

This condition guarantees that the process yˆt does not hit zero, which is one of the main reasons to use the Feller process in practice; it is relaxed shortly.

The corresponding Fokker–Planck problem has the form:

ϖtˉ12ε2yˉϖyˉyˉ+χκyˉϖyˉ=0,ϖt,y,t,yˉ=δyˉy.(7.47)

This equation can be written as a conservation law:

ϖtˉ+Fyˉ=0,ϖt,y,t,yˉ=δyˉy,(7.48)

where the probability flux F is given by

F=12ε2yˉϖyˉ+χκyˉϖ.(7.49)

However, experience suggests that solving the backward Kolmogorov problem is more expedient. It can be formulated as follows:

ϖt+12ε2yϖyy+χκyϖy=0,ϖtˉ,y,tˉ,yˉ=δyyˉ.(7.50)

The associated Kelvin wave function Kt,y,tˉ,yˉ,l has the following form:

K=expαt,tˉ+iγt,tˉyilyˉ,(7.51)

where α,γ solve the following system of backward ODEs:

αtt,tˉ+χiγt,tˉ=0,   αtˉ,tˉ=0,iγtˉt,tˉ12ε2γ2t,tˉκiγt,tˉ=0,   γtˉ,tˉ=l.(7.52)

Thus, γ solves a nonlinear Riccati equation, which can be linearized via the standard substitution

γt,tˉ=2iΩt,tˉε2Ωt,tˉ.(7.53)

As a result, one gets the following equations:

Ωttt,tˉκΩtt,tˉ=0,   Ωtˉ,tˉ=1,   Ωtˉ,tˉ=iε2l2,  (7.54)
αtt,tˉ+2χε2lnΩt,tˉt=0,   αtˉ,tˉ=0.(7.55)

Accordingly,

Ωt,tˉ=1ε22BκTil,(7.56)
Ωt,tˉ=ε22AκTil,(7.57)
γt,tˉ=AκTl1ε22BκTil,(7.58)
αt,tˉ=ϑ+1ln1ε22BκTil,(7.59)

and

K=expϑ+1ln1ε22BκTil+ε22AκT1ε22BκTilyyˉil.(7.60)

To analyze the problem further, it is helpful to define

M=2ε2BκT,(7.61)

introduce a new variable, llˆ:

lˆ=l2M,   l=2Mlˆ,(7.62)

and rescale K, KdlKˆdlˆ:

Kˆ=2MexpMY+yˉ×expϑ+1ln12ilˆ+MY12ilˆ+yˉ12ilˆ,(7.63)

where M appears due to the change of variables, and

Y=eκTy.(7.64)

Finally,

ϖt,y,tˉ,yˉ=MπeMyˉ+Yeϑ+1ln12ilˆ+MY12ilˆ+yˉ12ilˆdlˆ.(7.65)

Equation (7.65) allows us to understand the true meaning of condition (7.46). When this condition is satisfied, the corresponding integral converges absolutely when lˆ±. A well-known formula yields

ϖϑt,y,tˉ,yˉ=MeMyˉ+YyˉYϑ/2Iϑ2MyˉY.(7.66)

See, for example, Reference LiptonLipton (2001) and references therein. The probability flux F has the form

Fϑt,y,tˉ,yˉ=12ε2yˉϖϑt,y,tˉ,yˉyˉ+χκyˉϖϑt,y,tˉ,yˉ=12ε2MYϖϑ+1t,y,tˉ,yˉ+12ε2κMMyˉϖϑt,y,tˉ,yˉ.(7.67)

It is important to note that the density ϖyˉ integrates to one:

0ϖt,y,tˉ,yˉdyˉ=0euvvuϑ/2Iϑ2uvdv=1,(7.68)

where u=MY,   v=Myˉ. This fact is used in the following discussion.

Using the asymptotic expansion of the modified Bessel function, one can show that ϖϑ and F vanish on the boundary, since

ϖϑt,y,tˉ,yˉ=MeMYΓϑ+1Myˉϑ1+Oyˉ,Fϑt,y,tˉ,yˉ=ε221MYϑ+1κMMeMYMyˉϑ+1Γϑ+11+Oyˉ.(7.69)

Now assume that condition (7.46) is violated, so that 1<ϑ<0. In this case, the integral in (7.65) is no longer absolutely convergent, so one needs to regularize it. There are two ways of regularizing the corresponding integral: (I) integration by parts, (II) change of variables. Not surprisingly, they produce different results.

Start with integration by parts and write

Intϑ1πeϑ+1ln12ilˆ+MY12ilˆdeMyˉ12ilˆ2iMyˉ=1πϑ+1Myˉeϑ+2ln12ilˆ+MY12ilˆ+yˉ12ilˆdlˆ+1πYyˉeϑ+3ln12ilˆ+MY12ilˆ+yˉ12ilˆdlˆ,(7.70)

where the integrals are absolutely convergent. Thus, (7.66) yields

Intϑ=yˉYϑ22ϑ+1ZIϑ+1Z+Iϑ+2Z=yˉYϑ2IϑZ,(7.71)

where Z=2MyˉY, and a well-known recurrent relation for the modified Bessel functions is used; Reference Abramowitz and StegunAbramowitz and Stegun (1964), Eq. 9.6.26. Thus, Equations (7.66) and (7.67) hold for 1<ϑ<0:

ϖϑ,I=ϖϑ,   Fϑ,I=Fϑ.(7.72)

It is important to note that ϖϑ,Iyˉ0 when ϑ<0 (the corresponding singularity is integrable), while ϖϑ is bounded at yˉ=0, when ϑ>0. While the t.p.d.f. itself blows up at the natural boundary yˉ=0, the probability flux Fϑ,I vanishes on the boundary, so that the total probability of staying on the positive semiaxis is conserved.

Now, use change of variables to regularize Intϑ. Specifically, introduce l˜, such that

12ilˆ=112il˜,(7.73)

and formally write Intϑ as follows:

Intϑ=1πeϑ+1ln12il˜+MY12il˜+Myˉ12il˜dl˜=Yyˉϑ2Iϑ2MyˉY=yˉYϑ2Iϑ2MyˉY.(7.74)

Accordingly,

ϖϑ,IIt,y,tˉ,yˉ=MeMyˉ+YyˉYϑ2Iϑ2MyˉY.(7.75)

A straightforward calculation yields

Fϑ,IIt,y,tˉ,yˉ=MeMyˉ+YyˉYϑ/212ε2MyˉYIϑ+12MyˉY+12ε2ϑ+12ε2κMMyˉIϑ2MyˉY.(7.76)

It is easy to see that both ϖϑ,II and Fϑ,II are bounded at yˉ=0:

ϖϑ,IIt,y,tˉ,yˉ=MeMYΓϑ+1MYϑ1+Oyˉ,Fϑ,IIt,y,tˉ,yˉ=ε2ϑMeMY2Γϑ+1MYϑ1+Oyˉ.(7.77)

Since there is a probability flux across the natural boundary yˉ=0, the total probability on the positive semiaxis 0, is less than one.

Representative t.p.d.fs for Feller processes with different values of ϑ are illustrated in Figure 10.

Figure 10 T.p.d.fs for three Feller processes with different parameters and regularity conditions. (a) χ=0.1, κ=1.2, ε=0.2, y0=0.15, tˉmax=3; (b), (c) χ=0.1, κ=1.2, ε=0.6, y0=0.15, tˉmax=3. For the first and second processes, the probability of yˉ0 is equal to one. For the third process, this probability, shown as a function of time in (d), is less than one. Author’s graphics.

Feller Process with Time-Dependent Parameters

Surprisingly, studying the Feller process with time-dependent coefficients is viewed as a difficult problem, which remains an active area of research; see, for example, Reference MasoliverMasoliver (2016), Reference Giorno and NobileGiorno and Nobile (2021), and references therein. However, using Kelvin wave formalism allows one to find an expression for the t.p.d.f. in a very natural way.

For the process with time-dependent parameters, the problem of interest has the form:

ϖt+12ε2tyϖyy+χtκtyϖy=0,ϖtˉ,y=δyyˉ.(7.78)

Here it is assumed that the following regularity condition is satisfied:

ϑt=2χtε2t1>0.(7.79)

This condition guarantees that the corresponding integrals converge at infinity.

As usual, ϖ can be written as a superposition of Kelvin waves of the form

K=expαt,tˉ+iγt,tˉyilyˉ,(7.80)

where α,γ solve the following system of backward ODEs:

αtt,tˉ+χt,tˉiγt,tˉ=0,   αtˉ,tˉ=0,iγtˉt,tˉ12ε2tγ2t,tˉκtiγt,tˉ=0,   γtˉ,tˉ=l.(7.81)

Introducing Ωt,tˉ, such that

γt,tˉ=2iΩt,tˉε2tΩt,tˉ,(7.82)

one gets the following second-order equation for Ωt,tˉ:

Ωttt,tˉκ+2lnεΩtt,tˉ=0,   Ωtˉ,tˉ=1,   Ωtˉ,tˉ=ε2tˉ2il.(7.83)

Solving this equation, one gets

Ωt,tˉ=1ε2t2Bˉκt,tˉil,Ωtt,tˉ=ε2t2Aκt,tˉ2εtεtBˉκt,tˉil.(7.84)

Accordingly,

γt,tˉ=Aκt,tˉ2εtεtBˉκt,tˉl1ε2t2Bˉκt,tˉil,(7.85)
αt,tˉ=2χtε2tln1ε2t2Bˉκt,tˉilttˉ2χsε2sln1ε2s2Bˉκs,tˉilds.(7.86)

Thus, the Kelvin wave becomes

K=exp2χtε2tln1ε2t2Bˉκt,tˉilttˉ2χsε2sln1ε2s2Bˉκs,tˉilds+Aκt,tˉ2εtεtBˉκt,tˉ1ε2t2Bˉκt,tˉilyyˉil.(7.87)

By analogy with (7.61), (7.62), and (7.63), define

Mt,tˉ=2ε2tBˉκt,tˉ,   lˆ=l2Mt,tˉ,   l=2Mt,tˉlˆ(7.88)

and represent Kˆ as follows:

Kˆt,y,tˉ,yˉ,lˆ=2MexpMt,tˉY+yˉexp2χtε2tln12ilˆttˉ2χsε2sln1Mt,tˉMs,tˉ2ilˆds+Mt,tˉY12ilˆ+yˉ12ilˆ,(7.89)

where

Y=Aκt,tˉ4εtε3tMt,tˉy.(7.90)

Finally,

ϖt,y,tˉ,yˉ=12πKˆt,y,tˉ,yˉ,lˆdlˆ.(7.91)

Therefore, finding ϖt,y,tˉ,yˉ is reduced to solving some very simple ODEs and calculating a one-dimensional integral, which is theoretically appealing and numerically efficient.

Feller Process with Jumps

Consider a jump-diffusion process yˆt with constant coefficients governed by the following equation:

dyˆt=χκyˆtdt+εyˆtdZˆt+JdΠˆt,   yˆt=y,(7.92)

where Zˆt is a standard Wiener process, and Πˆt is a Poisson process with intensity λ. To preserve tractability, it is assumed that jumps are positive and exponentially distributed with parameter ϕ; for additional insights, see Reference Lipton and SheltonLipton and Shelton (2012).

The backward Kolmogorov problem can be written as

ϖt+12ε2yϖyy+(χky)ϖy+λ(ϕ0ϖ(t,y+J)eϕJdJϖ(t,y))=0,ϖ(tˉ,y,tˉ,yˉ)=δ(yyˉ).(7.93)

The corresponding Kelvin wave has the familiar form:

Kt,y,tˉ,yˉ,l=expαt,tˉ+iγt,tˉyilyˉ,(7.94)

where α,γ satisfy the following system of ODEs:

αtt,tˉ+χiγt,tˉ+λiγt,tˉϕiγt,tˉ,   αtˉ,tˉ=0,iγtˉt,tˉ12ε2γ2t,tˉκiγt,tˉ=0,   γ(tˉ,tˉ)=l.(7.95)

The expression for γ is given by (7.58), while α can be split as follows:

αt,tˉ=α0t,tˉ+λα1t,tˉ.(7.96)

In this setting, α0 has the familiar form:

α0t,tˉ=(ϑ+1)ln1ε22BˉκTil,(7.97)

while α1 can be represented as follows:

α1t,tˉ=ttˉAκtˉsilϕϕε22Bˉκtˉs+Aκtˉsilds=1κϕε22lnϕϕε22BˉκT+AκTilϕil.(7.98)

Thus, jumps do profoundly affect the dynamics of the underlying stochastic process.

7.4.2 Augmented Feller Process, I

This section studies the joint dynamics of a Feller process yˆt and its integral xˆt. The corresponding combined process is described by the following equations:

dxˆt=yˆtdt,   xˆt=x,dyˆt=χκyˆtdt+εyˆtdZˆt,   yˆt=y.(7.99)

Depending on the interpretation, these equations can describe the joint evolution of a particle’s position and its velocity, the integral of variance and variance, among other possibilities.

The forward Fokker–Planck has the following form:

ϖtˉ12ε2yˉϖyˉyˉ+yˉϖxˉ+χκyˉϖyˉ=0,ϖt,x,y,t,xˉ,yˉ=δxˉxδyˉy,(7.100)

while the backward Kolmogorov problem can be written as follows:

ϖt+12ε2yϖyy+yϖx+χκyϖy=0ϖtˉ,x,y,tˉ,xˉ,yˉ=δxxˉδyyˉ.(7.101)

In the following discussion the backward problem is considered, which allows one to derive the desired formula more efficiently. The corresponding function K has the following form:

K=expαt,tˉ+ikxxˉ+iγt,tˉyilyˉ,(7.102)

where

αtt,tˉ+iχγt,tˉ=0,   αtˉ,tˉ=0,iγtt,tˉ12ε2γ2t,tˉiκγt,tˉ+ik=0,   γtˉ,tˉ=l.(7.103)

As before, one can linearize the Riccati equation for γ by using substitution given by (7.53), with Ωt,tˉ solving the second-order equation of the following form:

Ωttt,tˉκΩtt,tˉ+iε2k2Ωt,tˉ=0,   Ωtˉ,tˉ=1,   Ωtˉ,tˉ=iε2l2.(7.104)

One can represent Ωt,tˉ in the following form:

Ωt,tˉ=ω+eλ+tˉt+ωeλ+tˉt,(7.105)

where λ± are solutions of the characteristic equation:

λ2+κλ+iε2k2=0,(7.106)

and ω± satisfy the following system of linear equations:

ω++ω=1,λ+ω++λω=iε2l2.(7.107)

Thus,

λ±=μ±ζ,μ=κ2,   ζ=κ22iε2k2,(7.108)
ω+=2λ+iε2l4ζ,   ω=2λ++iε2l4ζ.(7.109)

It is useful to note that

λ+λ=μ2ζ2=iε2k2.(7.110)

For the sake of brevity, notation (6.122) is used:

Ωt,tˉ=E02λ+iε2lE++2λ++iε2lE4ζ,(7.111)
Ωtt,tˉ=E0λ+2λ+iε2lE+λ2λ++iε2lE4ζ,(7.112)
γ=2iλ+2λ+iε2lE+λ2λ++iε2lEε22λ+iε2lE+2λ++iε2lE,(7.113)
α=χκTε2ϑ+1ln2λ+iε2lE++2λ++iε2lE4ζ.(7.114)

Accordingly, K can be written in the following form:

K=expχκTε2+ikxxˉϑ+1ln2λE++λ+Eiε2lE+E4ζ+22λ+λE+E+iε2lλ+E+λEε22λE++λ+Eiε2lE+Eyilyˉ.(7.115)

Define a new variable lˆ, such that

lˆ=l2M,   l=2Mlˆ,(7.116)

where

M=2λE++λ+Eε2E+E.(7.117)

Rescaled Kˆ can be factorized as follows:

Kˆ=Kˆ1Kˆ2,(7.118)

where

Kˆ1=expχκTε2+ikxxˉϑ+1lnλE++λ+E2ζ+2λ+λE+Eε2λE++λ+Ey,Kˆ2=2MexpMY+yˉ×expϑ+1ln12ilˆ+MY12ilˆ+yˉ12ilˆ,(7.119)

with

Y=4ζ2λE++λ+E2y.(7.120)

Integration with respect to lˆ can be done analytically:

MπKˆ2dlˆ=MeMyˉ+YyˉYϑ2Iϑ2MyˉY,(7.121)

which allows one to calculate ϖ via a single inverse Fourier transform:

ϖ=12πexpχκTε2+ikxxˉϑ+1lnλE++λ+E2ζ+2λ+λE+Eε2λE++λ+EyMeMyˉ+YyˉYϑ2Iϑ2MyˉYdk.(7.122)

A typical t.p.d.f. for a degenerate augmented Feller process is illustrated in Figure 11.

Figure 11 A thousand trajectories of a representative t.p.d.f. for the degenerate augmented Feller process. Parameters are T=3, dt=0.01, χ=0.1, κ=1.2, ε=0.2, x=0 ; y0=0.15. (a) xt, (b) yt, (c) xˉT,yˉT, (d) contour lines of ϖ0,0.15,0,T,x˜,y˜. Author’s graphics.

Since the integral over yˉ is equal to one, one can represent the marginal distribution of xˉ in the following form:

ϖxt,x,y,tˉ,xˉ=12πϝt,y,tˉ,keikxxˉdk,(7.123)

where

ϝt,y,tˉ,k=expχκTε2ϑ+1lnλE++λ+E2ζ+2λ+λE+Eε2λE++λ+Ey,(7.124)

with μ,ζ given by the equations in (7.106). It is easy to check that ϖxt,x,y,tˉ,xˉ integrates to one:

ϖxt,x,y,tˉ,xˉdxˉ=12πϝt,y,tˉ,kexpikxxˉdkdxˉ=ϝt,y,tˉ,kδkdk=ϝt,y,tˉ,0=expχκTε22χε2κT2=1.(7.125)

The expected value of xˉ has the following form:

X=ϖxt,x,y,tˉ,xˉxˉdxˉ=x+ϖxt,x,y,tˉ,xˉxˉxdxˉ=x+12πϝt,y,tˉ,kexpikxˉxxˉxdkdxˉ=x+limϵ0ddϵ12πϝt,y,tˉ,kexpik+ϵxˉxdkdxˉ=x+limϵ0ddϵϝt,y,tˉ,kδkiϵdk=x+ddϵϝt,y,tˉ,iϵϵ=0.(7.126)

A calculation left to the reader yields

X=x+χκTBˉκTχκy,(7.127)

which agrees with (6.115).

It is worth noting that ϖxtˉ,xˉ has fat tails, since some of the exponential moments of xˉ have finite-time explosions; see Reference Andersen and PiterbargAndersen and Piterbarg (2007), Reference Friz and Keller-ResselFriz and Keller-Ressel (2010), and references therein.Footnote 8 Specifically, one needs to analyze if Ipt,tˉ of the following form:

Ipt,tˉ=ϖxt,x,y,tˉ,xˉepxˉdxˉ=epx2πϝt,y,tˉ,kexpikxxˉpxxˉdkdxˉ=epxϝt,y,tˉ,kδikpdk=epxϝt,y,tˉ,ip,(7.128)

blows up for some finite tˉ>t. Indeed,

ϝt,y,tˉ,ip=expχκTε2ϑ+1lnλE++λ+E2ζ+2λ+λE+Eε2λE++λ+Ey,(7.129)

where

λ±=μ±ζ,μ=κ2,   ζ=κ22ε2p2,λ+λ=ε2p2.(7.130)

Thus, when ζ>0 is real:

p=(2ζμsinh(|ζ|T)+ζcosh(|ζ|T))(ϑ+1)exp(χκTε2+px+psinh(ζT)(μsinh(|ζ|T)+ζcosh(|ζ|T))y),(7.131)

and, when ζ=iζ is imaginary:

Ip=ζμsinζT+ζcosζTϑ+1expχκTε2+px+psinζTμsinζT+ζcosζTy.(7.132)

For p,pˆ, ζ is real, and for ppˆ,, it is imaginary. Here

pˆ=κ22ε2>0.(7.133)

There is no blowup when ζ is real. When ζ is imaginary, the blowup time t* is the smallest positive root of the equation

κsin2ε2pκ2t*t+2ε2pκ2cos2ε2pκ2t*t=0,(7.134)
t*=t+πarctan2ε2pκ2κ2ε2pκ2.(7.135)

It is clear that I1 does not blow up. This fact in used in the next section.

The marginal distribution of yˉ, ϖytˉ,yˉ is the standard Feller distribution given by (7.66).

7.4.3 Augmented Feller Process, II

This section studies the joint dynamics of an arithmetic Brownian xˆt whose stochastic variance is driven by a Feller process yˆt, and considers the following system of affine SDEs:

dxˆt=yˆtdWˆt,   xˆt=x,dyˆt=χκyˆtdt+εyˆtdZˆt,   yˆt=y.(7.136)

Studying such a process is very helpful for finding option prices and solving other important problems in the financial engineering context.

The associated forward Fokker–Planck problem can be written as follows:

ϖtˉ12yˉϖxˉxˉρεyˉϖxˉyˉ12ε2yˉϖyˉyˉ+χκyˉϖyˉ=0,ϖt,x,y,t,xˉ,yˉ=δxˉxδyˉy,(7.137)

while the backward Kolmogorov problem has the following form:

ϖt+12yϖxx+ρεyϖxy+12ε2yϖyy+χκyϖy=0,ϖtˉ,x,y,tˉ,xˉ,yˉ=δxxˉδyyˉ.(7.138)

As before, concentrate on problem (7.138).

The Kelvin function K has the form (7.102). The governing ODEs for α,γ are as follows:

αtt,tˉ+iχγt,tˉ=0,   αtˉ,tˉ=0,iγtt,tˉ12ε2γ2t,tˉκiρεkiγt,tˉ12k2=0,   γtˉ,tˉ=l.(7.139)

Formulas (7.111)–(7.114) hold; however, the corresponding characteristic equation is

λ2+κiρεkλε24k2=0,(7.140)

so that

λ±=μ±ζ,μ=12κiρεk,   ζ=12ε2ρˉ2k22iρεκk+κ2,λ+λ=μ2ζ2=ε24k2,(7.141)

where ρˉ2=1ρ2. Subsequent calculations are very similar to the ones performed in the previous subsection, so they are omitted for brevity. The final expressions for ϖ and ϖx are given by Equations (7.122), (7.123), and (7.124), with μ,ζ given by the equations in (7.141). These expressions are similar to the formulas originally derived by Lipton as part of his analysis of the Heston stochastic volatility model; see Reference LiptonLipton (2001).Footnote 9

A typical t.p.d.f. for a nondegenerate augmented Feller process is shown in Figure 12.

Figure 12 A thousand trajectories of a representative nondegenerate augmented Feller process. Parameters are T=3, dt=0.01, χ=0.2, κ=2.0, ε=0.2, ρ=0.5, x=0 , y0=0.15. (a) xt, (b) yt, (c) xˉT,yˉT, (d) contour lines of ϖ0,0.15,0,T,x˜,y˜. Author’s graphics.

As before, ϖxt,x,y,tˉ,xˉ has fat tails. Consider Ipt,tˉ given by (7.128). The corresponding λ± have the following form:

λ±=μ±ζ,μ=12κρεp,   ζ=12ε2ρˉ2p22ρεκp+κ2,λ+λ=μ2ζ2=ε24p2.(7.142)

Thus, when ζ>0 is real,

Ip=ζμsinhζT+ζcoshζTϑ+1exp2χμTε2+px+pp1sinhζT2μsinhζT+ζcoshζTy,(7.143)

and when ζ=iζ is imaginary,

Ip=ζμsinζT+ζcosζTϑ+1exp2χμTε2+px+pp1sinζT2μsinζT+ζcosζTy.(7.144)

One needs to determine when ζ becomes imaginary. The corresponding quadratic equation has the form:

ρˉ2ε2p2+2ρεκpκ2=0,(7.145)

its roots are as follows:

p±=ρεκ±ρ2ε2κ2+ρˉ2ε2κ2ρˉ2ε2=ρ±1κρˉ2ε,(7.146)

so that

p+>1,   p<1.(7.147)

For pp,p+, ζ is real, for pp,p+, it is imaginary. There is no blowup when ζ is real. When ζ is imaginary, the blowup time t* is the smallest positive root of the equation

μsinζt*t+ζcosζt*t=0,(7.148)
t*=t+arctanζμζ,μ>0,t+πarctanζμζ,μ<0.(7.149)

7.5 Example: Path-Dependent Process

Let yˆt be a stochastic process and xˆt be its moving average. Then

xˆt=κteκtsyˆsds.(7.150)

A simple calculation yields

dxˆt=κyˆtxˆtdt.(7.151)

The process yˆt is path-dependent, because its volatility σˆt depends on its moving average xˆt:

σˆt=a0+a1yˆtxˆt,(7.152)

where a0>0, a1<0, in order to capture the effect of leverage. Thus, one can write the governing degenerate system of SDEs as follows:

dxˆt=κyˆtxˆtdt,   xˆt=x,dyˆt=a0+a1yˆtxˆtdWˆt,   yˆt=y.(7.153)

The Fokker–Planck and Kolmogorov problems are

ϖtˉ12a0+a1yˉxˉϖyˉyˉ+κyˉxˉϖxˉ=0,ϖt,x,y,t,xˉ,yˉ=δxˉxδyˉy,(7.154)
ϖt+12a0+a1yxϖyy+κyxϖx=0,ϖtˉ,x,y,tˉ,xˉ,yˉ=δxxˉδyyˉ,(7.155)

respectively.

A representative Kelvin mode has the following form:

K=expαt,tˉ+iβt,tˉxikxˉ+iγt,tˉyilyˉ.(7.156)

The system of backward ODEs for α,γ,β is as follows:

αtt,tˉa02γ2t,tˉ=0,   αtˉ,tˉ=0,iβtt,tˉ+a12γ2t,tˉiκβt,tˉ=0,   βtˉ,tˉ=k,iγtt,tˉa12γ2t,tˉ+iκβt,tˉ=0,   γtˉ,tˉ=l.(7.157)

The equations in (7.157) are matrix Riccati equations, as opposed to the scalar Riccati equations considered earlier. In general, such equations are very difficult to solve. However, the case under consideration is one of the relatively rare instances when a matrix Riccati equation can be solved explicitly. Start with an observation:

γtt,tˉ+βtt,tˉ=0,(7.158)

so that

βt,tˉ=γt,tˉ+k+l.(7.159)

Accordingly,

iγtt,tˉa12γ2t,tˉiκγt,tˉ+iκk+l=0,   γtˉ,tˉ=l.(7.160)

One can use Equations (7.111)–(7.113) with β,k replaced by γ,l, and

λ2+κλ+ia1κk+l2=0,(7.161)

so that

λ±=μ±ζ,μ=κ2,   ζ=κ22ia1κk+l2.(7.162)

Equation (7.159) yields

βt,tˉ=2iΩt,tˉa1Ωt,tˉ+k+l,(7.163)

and

γ2t,tˉ=2ia1γtt,tˉ+κβt,tˉ=2ia1γtt,tˉκγt,tˉ+κκ+k.(7.164)

Thus,

αtt,tˉ=ia0a1γtt,tˉκγt,tˉ+κk+l,   αtˉ,tˉ=0.(7.165)

Accordingly,

αt,tˉ=ia0a1γt,tˉl2a0κa12lnΩt,tˉ+a0κTa1ik+l.(7.166)

Finally,

K=expα+iβt,tˉxikxˉ+iγt,tˉyilyˉ=exp2a0κa12lnΩt,tˉ+iγt,tˉyx+a0a1+ ilxyˉ+a0a1κT1+ikxxˉ+a0κTa1.(7.167)

To make sure that σˆt given by (7.152) and the integrand (7.167) are well defined, it is assumed that

a0+a1yx>0,  a0+a1 yˉxˉ>0.(7.168)

7.6 Example: OU-Like Process

This section considers several instances when an OU-inspired process becomes non-Gaussian. This can happen for a variety of reasons, such as effects of anomalous diffusion, the presence of jumps, effects of augmentation, and the likes.

7.6.1 Anomalous OU Process

This section considers a mean-reverting process driven by a non-Gaussian anomalous diffusion. For brevity, it is assumed that coefficients are time-independent. The fractional forward Fokker–Planck and backward Kolmogorov problems can be written as follows:

ϖtˉ+a2yˉ21/2ϖ+χκyˉϖyˉ=0,ϖt,y,t,yˉ=δyˉy,(7.169)
ϖta2y21/2ϖ+χκyϖy=0,ϖtˉ,y,tˉ,yˉ=δyyˉ,(7.170)

respectively. Here a>0 is the anomalous diffusion coefficient.

As before, one can use Kelvin waves to solve (7.170) by choosing a particular solution of the form (7.51). The corresponding α,γ satisfy the following ODEs:

αtt,tˉaγt,tˉ+iχγt,tˉ=0,   αtˉ,tˉ=0,γtt,tˉκγt,tˉ=0,   γtˉ,tˉ=l,(7.171)

so that

αt,tˉ=BˉκTaliχl,γt,tˉ=eκTl.(7.172)

Accordingly,

ϖt,y,tˉ,yˉ=12πexpBˉκTal+BˉκTχ+eκTyyˉildl=1πBˉκTaBˉκTa2+eκTyχκyˉχκ2.(7.173)

Thus, in sharp contrast to the classical OU process, which is described by a Gaussian distribution, the fractional OU process is described by a Cauchy distribution. This distribution has fat tails and no first and second moments.

7.6.2 Non-Gaussian Augmented OU Process, I

On occasion, problems seemingly not of the type given by (7.11) can be cast in the proper form via a suitable trick. Consider, for example, the following system of SDEs:

dxˆt=yˆt2dt,   xˆt=x,dyˆt=χκyˆtdt+εdZˆt,   yˆt=y.(7.174)

Superficially, it does not belong to the class of processes studied earlier. However, by introducing new variables z1=x, z2=y2, z3=y, one can augment the equations in (7.127) as follows:

dzˆ1,t=zˆ2,tdt,   zˆ1,t=xz1,dzˆ2,t=ε22κzˆ2,t+2χzˆ3,t+2εzˆ3,tdZˆt,   zˆ2,t=y2z2,dzˆ3,t=χκzˆ3,tdt+εdZˆt,   zˆ3,t=z3y.(7.175)

These equations are “almost” in the suitable form. The only snag is that one cannot claim that zˆ3,t=zˆ2,t since zˆ3,t is not always positive.

The corresponding Fokker–Planck and Kolmogorov problems can be written as follows:

ϖtˉ2ε2zˉ2ϖzˉ2zˉ22ε2zˉ3ϖzˉ2zˉ312ε2ϖzˉ3zˉ3+zˉ2ϖzˉ1+ε22κzˉ2+2χzˉ3ϖzˉ2+χκzˉ3ϖzˉ3=0,(7.176)
ϖt,x,y2,y,t,zˉ1,zˉ2,zˉ3=δzˉ1xδzˉ2y2δzˉ3y,ϖt+2ε2z2ϖz2z2+2ε2z3ϖz2z3+12ε2ϖz3z3+z2ϖz1+(ε22κz2+2xz3)ϖz2+(χκz3)ϖz3=0,ϖtˉ,z1,z2,z3,tˉ,zˉ1,zˉ32,zˉ3=δz1zˉ1δz2zˉ32δz3zˉ3.(7.177)

As usual, K has the form:

Kt,tˉ,z,m=expαt,tˉ+im1z1zˉ1+iδ2t,tˉz2im2zˉ32+iδ3t,tˉz3im3zˉ3.(7.178)

The corresponding set of ODEs for α,δ2,δ3 is as follows:

αtt,tˉε22δ32t,tˉ+iε2δ2t,tˉ+iχδ3t,tˉ=0,   αtˉ,tˉ=0,iδ2t,tˉ2ε2δ22t,tˉ2iκδ2t,tˉ+im1=0,   δ2tˉ,tˉ=m2,iδ3t,tˉ2ε2δ2t,tˉδ3t,tˉ+2iχδ2viκδ3t,tˉ=0,   δ3tˉ,tˉ=m3.(7.179)

These are matrix Riccati equations.

Once again, the corresponding matrix Riccati equation can be solved explicitly. Since the second equation is separable and hence can be viewed as a scalar Riccati equation, one can start with a familiar ansatz and use Equations (7.111)–(7.113) with γ,l replaced by δ2,m2, and the corresponding characteristic equation is as follows:

λ2+2κλ+2iε2m1=0,(7.180)

and its solutions have the familiar form:

λ±=μ±ζ,μ=κ,   ζ=κ22iε2m1.(7.181)

To linearize the equations in (7.179) as a whole, use the following ansatz:

Ω=E0ω+E++ωE,α=12lnΩ+E0a0+a+E++aEΩ+gtˉt, δ2=E0b+E++bEΩ,   δ3=E0c0+c+E++cEΩ,(7.182)

where a0, a±, b0, b±, and g are constants to be determined. This ansatz is useful since terms proportional to E0,E+,E balance each other, which allows us to find the coefficients explicitly. Initial conditions complete the picture. The actual calculation is omitted for brevity. The result is as follows:

ω±=λ+2iε2m22ζ,   b±=iλ±ω±2ε2,c±=±iχλ±ω±ε2ζ,   c0=m3c+c,    g=χ2λ+λ2ε2ζ2,   a0=iκχc0ζ2,   a±=ω±a0ε2c024ζ+χ2κ2ω+ωε2ζ3,(7.183)

where λ± are given by the equations in (7.181). These expressions can be substituted in the function K to obtain the corresponding t.p.d.f.

7.6.3 Non-Gaussian Augmented OU Process, II

This section studies an affine process of the following form:

dxˆt=yˆtdWˆt,   xˆt=x,dyˆt=χκyˆtdt+εdZˆt,   yˆt=y.(7.184)

The killed process is studied in Section 8 in the context of the Stein–Stein model.

Precisely as before, one can introduce the new variables z1=x, z2=y2, z3=y, and expand the equations in (7.184) as follows:

dzˆ1,t=zˆ3,tdWˆt,   zˆ1,t=xz1,dzˆ2,t=ε22κzˆ2,t+2χzˆ3,t+2εzˆ3,tdZˆt,   zˆ2,t=y2z2,dzˆ3,t=χκzˆ3,tdt+εdZˆt,   zˆ3,t=z3y.(7.185)

It is clear that the equations in (7.185) are affine.

The corresponding Fokker–Planck and Kolmogorov problems can be written as follows:

ϖtˉ12zˉ2ϖzˉ1zˉ12ρεzˉ2ϖzˉ1zˉ2ρεzˉ3ϖzˉ1zˉ32ε2zˉ2ϖzˉ2zˉ22ε2zˉ3ϖzˉ2zˉ312ε2ϖzˉ3zˉ3+ε22κzˉ2+2χzˉ3ϖzˉ2+χκzˉ3ϖzˉ3=0,(7.186)
ϖt,x,y2,y,t,zˉ1,zˉ2,zˉ3=δzˉ1xδzˉ2y2δzˉ3y,ϖt+12z2ϖz1z1+2ρεz2ϖz1z2+ρεz3ϖz1z3+2ε2z2ϖz2z2+2ε2z3ϖz2z3+12ε2ϖz3z3+ε22κz2+2χz3ϖz2+χκz3ϖz3=0,ϖtˉ,z1,z2,z3,tˉ,zˉ1,zˉ32,zˉ3=δz1zˉ1δz2zˉ32δz3zˉ3.(7.187)

One can use K given by (7.178) and write the set of ODEs for α,δ2,δ3 as follows:

αtt,tˉε22δ32t,tˉ+iε2δ2t,tˉ+iχδ3t,tˉ=0,   αtˉ,tˉ=0,iδ2(t,tˉ)2ε2δ22(t,tˉ)2i(κiρεm1)δ2(t,tˉ)12m12=0,   δ2(tˉ,tˉ)=m2,iδ3t,tˉ2ε2δ2t,tˉδ3t,tˉ+2iχδ2t,tˉiκiρεm1δ3t,tˉ=0,δ3tˉ,tˉ=m3.(7.188)

As before, this system can be linearized and solved analytically, which was pointed out by Reference Stein and SteinStein and Stein (1991), Reference Schöbel and ZhuSchöbel and Zhu (1999). One can repeat the result obtained in the previous section verbatim, except for (7.181). The corresponding characteristic equation has the following form:

λ2+2κiρεm1λε2m12=0,(7.189)

and its solutions can be written as follows:

λ±=μ±ζ,μ=κiρεm1,   ζ=ρˉ2ε2m122iρεκm1+κ2.(7.190)

The rest of the formal analysis is the same. But the asymptotic behavior of the t.p.d.f. is, of course, different.

8 Pricing of Financial Instruments

8.1 Background

The formulas derived in Sections 6 and 7 can be used to solve numerous problems of financial engineering within a consistent framework based on Kelvin waves. Here are some representative examples.

Payoffs of European options depend solely on the terminal value of Sˉ=Sˆtˉ of the underlying price at the option’s maturity. The most common European options are calls and puts, but, on occasion, binary options and other types are traded as well. Since the hedging and speculation needs of market participants cannot be satisfied by European options alone, the whole industry emerged to design, price, and hedge the so-called exotic options, with payoffs depending on the entire price trajectory between inception and maturity.

Prices of the fundamental financial instruments, such as forwards and European calls and puts, depend on the underlying prices only at maturity. However, the prices of many other instruments depend on the entire underlying price history between the instrument’s inception and maturity. Typical examples are barrier, American, Asian, lookback, and passport options; see, for example, Reference Lipton-LifschitzLipton-Lifschitz (1999), Reference LiptonLipton (2001), and references therein. Moreover, the prices of bonds also depend on the history of the interest rates and credit spreads throughout their life. This section shows how to price some path-dependent financial instruments using the methodology developed in the previous sections.

8.2 The Underlying Processes

The original approach to modeling financial assets was developed by Bachelier, who assumed that prices Sˆt of such instruments are governed by an arithmetic Brownian motion; see Reference BachelierBachelier (1900):

dSˆt=rSˆtdt+σˆdWˆt,   Sˆt=S.(8.1)

Here, r is the risk-neutralized drift, σˆ is the volatility, and Wˆt is a Wiener process; r, σˆ are dimensional quantities, r=T1, σ=$T1/2. The process for Sˆt given by (8.1) is affine; in fact, it is an OU process with zero mean and mean-repulsion instead of mean-reversion.

Subsequently, the academic community concluded that using a geometric Brownian motion as a driver is more appropriate; see Reference BonessBoness (1964), Reference SamuelsonSamuelson (1965), Reference Black and ScholesBlack and Scholes (1973), and Reference MertonMerton (1973). At present, the basic assumption is that the price Sˆt of an underlying financial instrument follows a geometric Brownian motion process with constant coefficients:

dSˆtSˆt=rdt+σdWˆt,   Sˆt=S.(8.2)

Here, r is the risk-neutralized drift, and σ is the volatility. These are dimensional quantities, r=T1, σ=T1/2.

The choice between using the Bachelier and the Black–Scholes models often depends on the nature of the underlying asset and the market’s specific characteristics. Since the Bachelier model assumes that the underlying asset prices follow a normal distribution, it can be more appropriate for assets whose price changes are additive and can theoretically go below zero, like interest rates, some commodities, or certain types of bonds. Generally, the price movements of the underlying asset are relatively small for short periods, so the Bachelier model provides a good description of these movements. The Bachelier model is often used for pricing commodities, some interest-rate derivatives, and studying the optimal execution. In markets with relatively low volatility, the Bachelier model’s assumption of additive price movements can provide a better fit for pricing and hedging derivatives than the multiplicative approach of the Black–Scholes model.

It was realized, very soon after the seminal paper by Reference Black and ScholesBlack and Scholes (1973) was published, that in practice it provides a rather poor description of reality. Hence, considerable efforts were dedicated to developing more adequate models. Such models include the jump-diffusion, local volatility, path-dependent volatility, stochastic volatility, local-stochastic volatility, rough volatility, and culminate in the universal volatility model; see Reference MertonMerton (1976), Reference Stein and SteinStein and Stein (1991), Reference Bick and ReismanBick and Reisman (1993), Reference HestonHeston (1993), Reference Derman and KaniDerman and Kani (1994), Reference DupireDupire (1994), Reference RubinsteinRubinstein (1994), Reference Hobson and RogersHobson and Rogers (1998), Reference Jex, Henderson and WangJex et al. (1999), Reference LewisLewis (2000), Reference LiptonLipton (2000, Reference Lipton2001); Reference Boyarchenko and LevendorskiiBoyarchenko and Levendorsky (2002), Reference Hagan, Kumar, Lesniewski and WoodwardHagan et al. (2002), Reference LiptonLipton (2002), Reference BergomiBergomi (2015), Reference ReghaiReghai (2015), Reference Gatheral, Jaisson and RosenbaumGatheral et al. (2018), Reference Gershon, Lipton, Rosenbaum and WienerGershon et al. (2022), and references therein.

Replacing constant volatility for a geometric Brownian motion with stochastic volatility driven by a Feller process results in the popular Heston model; see Reference HestonHeston (1993). This model has numerous applications, particularly for pricing equity and foreign exchange derivatives. The governing SDEs are as follows:

dSˆtSˆt=rdt+vˆtdWˆt,   Sˆt=S,dvˆt=χκvˆtdt+εvˆtdZˆt,   vˆt=ν,(8.3)

where dWˆtdZˆt=ρdt. The logarithmic change of variables, given by (8.3), yields the equations of (7.136).

Replacing constant volatility with stochastic volatility driven by an OU process results in the (less popular) Stein–Stein model; see Reference Schöbel and ZhuSchöbel and Zhu (1999); Reference Stein and SteinStein and Stein (1991). The corresponding SDEs have the form:

dSˆtSˆt=rdt+σˆtdWˆt,   Sˆt=S,dσˆt=χκσˆtdt+εdZˆt,   vˆt=ν,(8.4)

Reference Stein and SteinStein and Stein (1991) considered the special case of zero correlation, dWˆtdZˆt=0, while Reference Schöbel and ZhuSchöbel and Zhu (1999) studied the general case of arbitrary correlation, dWˆtdZˆt=ρdt.

Now, it is shown how to use formulas derived in Sections 6 and 7 in the context of financial engineering.

8.3 European Derivatives

8.3.1 Forwards, Calls, Puts, and Covered Calls

The most basic derivatives are forwards. Recall that a forward contract obligates the buyer (seller) to buy (to sell) an underlying asset for an agreed price at a specified future date. These contracts are not standardized and are traded over-the-counter (OTC), not on exchanges. Typical underlying assets are commodities, currencies, and financial instruments. The choice of an asset depends on the needs of the contracting parties. The price agreed upon in a forward contract is called the forward price. This price is derived based on the spot price of the underlying asset, adjusted for factors like time to maturity, interest rates, and dividends. Forward contracts are primarily used for hedging price fluctuations of the underlying asset or speculation. The payoff of a forward contract with maturity tˉ and strike K has the following form:

UFSˉ,K=SˉK,(8.5)

where the strike is chosen in such a way that today’s price of the forward contract is equal to zero. This price can be found without knowing the actual stochastic process Sˆ. The hedging argument shows that the only way to deliver the price of a non-dividend-paying stock at maturity tˉ is to buy it outright at inception t. Similarly, to deliver the strike K at time tˉ, one has to buy a zero coupon bond at time t. Let Zt,tˉ be the price of a bond paying unity at maturity tˉ. Then

Ft,tˉK=SZt,tˉ.(8.6)

In contrast to forwards, a European call option grants the holder the right, but imposes no obligation, to buy an underlying asset at the option maturity for a predetermined strike price. Similarly, a European put option grants the holder the right to sell an underlying asset. Theoretically, buyers utilize calls and puts to hedge future risks; however, they often buy options for speculative purposes. American options can be exercised at any time of the buyer’s choice before the option’s maturity. Bermudan options are exercisable at fixed times between their inception and maturity. A call option is a contract between two parties – a buyer and a seller. Typically, the buyer takes the long position on the underlying (i.e., she expects that at maturity, the underlying price will exceed the strike price) and does not hedge her position. On the other hand, the seller or writer of the option (typically a bank) does hedge and, hence, maintains a market-neutral position. The seller receives cash up-front but incurs potential liabilities at option maturity if the option is exercised. In contrast, the buyer pays money up front in exchange for the potential for future gains. For a put option, the buyer takes a short position, while the seller is still market-neutral.

Payoffs of call and put options with maturity tˉ and strike K have the form

UCSˉ,K=maxSˉK,0,UPSˉ,K=maxKSˉ,0,UC,PSˉ,K=maxϕSˉK,0,(8.7)

where ϕ=1 for a call, and ϕ=1 for a put. Put-call parity implies that their difference is linear in Sˉ and represents a forward contract:

UCSˉ,KUPSˉ,K=SˉK.(8.8)

Several popular models, including Bachelier, Black–Scholes, Heston, and Stein–Stein, are considered below. While the Bachelier model is not scale invariant, all the other models are. A general driver for a scale-invariant model can be written as follows:

dSˆtSˆt=rdt+σtdWˆt+υdΠˆt,   Sˆt=S,(8.9)

where, potentially, the volatility σˆt and the intensity λˆt of the Poisson process Πˆt are driven by SDEs of their own. For such models, it is convenient to decompose call and put payoffs (8.43) into parts, which are easier to study via Kevin waves; see Reference LiptonLipton (2001), Reference Lipton(2002). To this end, introduce the covered call with the payoff of the form

UCCSˉ,K=minSˉ,K.(8.10)

The call and put payoffs can be decomposed as follows:

UCSˉ,K=SˉUCCSˉ,K,   UPSˉ,K=KUCCSˉ,K(8.11)

Thus, the call price is the difference between the forward price and the covered call price, while the put price is the difference between the bond price and the covered call price. In both cases, the covered call is the source of optionality.

8.3.2 Black–Scholes Model

For the standard log-normal process, the backward pricing problem for covered calls can be written as follows:

Ut+12σ2S2USS+rUSrU=0,Utˉ,S=minS,K.(8.12)

It is helpful to rewrite it by using forward rather than spot prices:

Uˆt+12σ2F2UˆFF=0,Uˆtˉ,F=minF,K,(8.13)

where

Fˆt,tˉ=ertˉtSˆt,   Uˆt,F=ertˉtUt,S.(8.14)

Change of variables,

Fˆt,tˉxˆt,tˉ,   Fˆt,tˉ=Kexˆt,tˉ,(8.15)

results in the following process for xˆt:

dxˆt,tˉ=12σ2dt+σdWˆt,   xˆt,tˉ=x=lnFt,tˉK.(8.16)

The t.p.d.f. for this process is Gaussian:

ϖt,x,tˉ,xˉ=12πσ2Texpxˉx+σ2/2T22σ2T.(8.17)

Since the the nondimensional payoff of the covered call has the form

U˜CCx=minex,1,(8.18)

where U˜=Uˆ/K, one obtains the following expression for U˜CC:

U˜CCt,x=exNxσTσT2+NxσTσT2,(8.19)

where N. is the cumulative normal function.

By using (8.19), one can represent call and put prices as follows:

UˆC,Pt,FT=ϕFTNϕd+KNϕd,d±=lnFT/KσT±σT2.(8.20)

See Reference BlackBlack (1976).

Returning to the original variables, write the classical Reference Black and ScholesBlack and Scholes (1973) closed-form formula for the time t prices of calls and puts in its original form:

UC,Pt,S=ϕSNϕd+erTKNϕd,d±=lnerTS/KσT±σT2.(8.21)

Further transforming,

U˜CCt,x=ex/2VCCt,x,(8.22)

yields the following backward problem:

VtCC+12σ2VxxCC18σ2VCC=0,VCCtˉ,x=ex/2,(8.23)

with symmetric “peakon” payoff, which is proportional to the Laplace distribution density. This transform removes the drift in the x direction at the expense of adding killing with intensity σ2/8. Equation (8.19) implies

Vt,x=ex/2NxσTσT2+ex/2NxσTσT2.(8.24)

The Fourier transform of the “peakon” payoff yields

ex/2ikxdx=1k2+1/4.(8.25)

By using this formula, one can derive an alternative expression for UC,P based on Kelvin waves; see Reference LiptonLipton (2002). It is clear that Kelvin waves associated with the killed arithmetic Brownian motion described by (8.16) are the standard Fourier waves of the following form:

Kt,x,k=ek2+1/4σ2T/2+ikx.(8.26)

Equations (8.25) and (8.26) yield the following alternative expression for the price of covered calls given by (8.24):

VCCt,x=12πek2+1/4σ2T/2+ikxk2+1/4dk.(8.27)

See Reference LiptonLipton (2002). Equation (8.27) is central for the subsequent developments. For a single strike, this formula is less efficient than its classical counterpart; however, for a set of strikes, it is faster, because all the prices can be computed in one go, via the Fast Fourier Transform.

As one shall see shortly, these formulas help to handle affine pricing models very naturally.

8.3.3 Heston Model

The transformed forward pricing problem for the Heston model with the “peakon” payoff has the following form:

VtCC+12yVxxCC+2ρεVxyCC+ε2VyyCC+χκˆyVyCCy8VCC=0,VCCtˉ,x,y=ex/2,(8.28)

where κˆ=κρε/2. Thus, one is dealing with the killed stochastic process given by the equations in (7.136). Adapting the corresponding equations to accommodate the updated mean-reversion rate and the presence of the killing term, one gets the following system of ODEs for the corresponding Kelvin wave parameters:

αtt,tˉ+iχγt,tˉ=0,   αtˉ,tˉ=0,iγtt,tˉ12ε2γ2t,tˉκˆiρεkiγt,tˉ12k2+14=0,   γtˉ,tˉ=0.(8.29)

Formulas (7.111)–(7.114) are still applicable. However, the corresponding characteristic equation and its solution are:

λ2+κˆiρεkλε24k2+14=0,λ±=μ±ζ,μ=κˆiρεk2,   ζ=ρˉ2ε2k22iρεk+κˆ2+ε2/42.(8.30)

It is convenient to write α,γ as follows:

αT,k=2χε2μ+ζT+lnμ+ζ+μ+ζe2ζT2ζ,(8.32)
γT,k=k2+14i1e2ζT2μ+ζ+μ+ζe2ζTk2+14iςT,k.(8.33)

Hence, the price of the “peakon” has the following form:

VCCt,x,y=12πeαT,kk2+1/4ςT,ky+ikxk2+1/4dk.(8.34)

Equation (8.34) is frequently called the Lewis–Lipton formula; see, for example, Reference LewisLewis (2000), Reference LiptonLipton (2000), Reference LewisLewis (2001), Reference LiptonLipton (2001), Reference Lipton(2002); Reference SchmelzleSchmelzle (2010), Reference Janek, Kluge, Weron and WystupJanek et al. (2011).

The implied volatility surface generated by a representative Heston model is shown in Figure 13. Recall that the implied volatility ΣT,K is the volatility one must substitute into the Black–Scholes formula to reproduce the market price of a call (or put) option with maturity T and strike K. Thus, the deviation of the volatility surface from the flat surface ΣT,K=Σ0 shows how far a given market (or model) is from the idealized Black–Scholes framework.

Figure 13 A representative implied volatility surface generated by the Heston model. Parameters are the same as in Figure 12. Author’s graphics.

8.3.4 Stein–Stein Model

The transformed forward pricing problem for the Stein–Stein model with the “peakon” payoff has the following form:

VtCC+12z2Vz1z1CC+2ρεz2Vz1z2CC+ρεz3Vz1z3CC+2ε2z2Vz2z2CC+2ε2z3Vz2z3CC+12ε2Vz3z3CC+ε22κˆz2+2χz3Vz2CC+χκˆz3Vz3CCz28VCC=0,VCCtˉ,z1,z2,z3=ez1/2,(8.35)

which corresponds to the killed stochastic process described by the equations in (7.184). By incorporating the killing term, one gets the following set of ODEs for the Kelvin wave parameters

αtt,tˉε22δ32t,tˉ+iε2δ2t,tˉ+iχδ3t,tˉ=0,   αtˉ,tˉ=0,iδ2t,tˉ2ε2δ22t,tˉ2iκˆiρεm1δ2t,tˉ12m12+14=0,   δ2tˉ,tˉ=0,iδ3t,tˉ2ε2δ2t,tˉδ3t,tˉ+2iχδ2t,tˉiκˆiρεm1δ3t,tˉ=0,   δ3tˉ,tˉ=0.(8.36)

The corresponding solution has the form given by the equations in (7.182) with

λ±=μ±ζ,μ=κˆiρεm1,   ζ=ρˉ2ε2m122iρεκm1+κ2+ε2/4,ω±=λ2ζ,   b±=iλ±ω±2ε2,c±=±iχλ±ω±ε2ζ,   c0=c+c,    g=χ2λ+λ2ε2ζ2,   a0=iκχc0ζ2,   a±=ω±a0ε2c024ζ+χ2κ2ω+ωε2ζ3.(8.37)

The generic expression for the price of the “peakon” has the following form:

VCCt,z1,z32,z3=12πeαT,m1+iδ2T,m1z32+iδ3T,m1z3+im1z1m12+1/4dm1.(8.38)

It is clear that this price is a function of tˉ, z1, z3.

8.3.5 Path-Dependent Volatility Model

Reference Hobson and RogersHobson and Rogers (1998) initially proposed path-dependent volatility models; subsequently, they were studied by many authors; see Reference DavisDavis (2004), Reference Di Francesco and PascucciDi Francesco and Pascucci (2004, Reference Di Francesco and Pascucci2005); Guyon Reference Guyon(2014), and Reference Lipton and ReghaiLipton and Reghai (2023), among others. They present a viable alternative to the more popular local volatility models developed by Reference Bick and ReismanBick and Reisman (1993), Reference Derman and KaniDerman and Kani (1994), Reference DupireDupire (1994), and Reference RubinsteinRubinstein (1994).

The main advantage of path-dependent volatility models compared to their local volatility brethren is that the former deal with volatility functions depending on a nondimensional argument, such as Sˆt/Aˆt, where Sˆt is the stock price, and Aˆt is its average, say, σ=σ(Sˆt/Aˆt), while the latter use volatilities depending on a dimensional argument Sˆt, σ=σ(Sˆt), which is conceptually unsound and results in model dynamics deviating from the one observed in the market. The problem with path-dependent models is that building an analytically tractable path-dependent model is exceedingly tricky, so gaining the necessary intuition or benchmarking numerical solutions is complicated. However, this section develops such a model using results derived in Section 7.3.

Here, an original path-dependent model with a semianalytical solution is presented for the first time. The dynamics is adapted from Section 7.3, Equation (7.153) as follows:

Aˆt=expκteκttlnSˆtdt,   Aˆt=A,dSˆtSˆt=c0+c1lnSˆtAtdWˆt,   Sˆt=S.(8.39)

It is not necessary to describe in detail how Sˆt˜, and, hence, Aˆt˜, behave when t˜<t, since it becomes unimportant provided that κT is sufficiently large. For instance, one can assume that Sˆt˜A, when t˜<t, then A=S. Additionally, it is assumed that r=0, so that spot and forward prices coincide, Sˆt=Fˆt,tˉ.

In logarithmic variables xˆt=lnAˆt, yˆt=lnSˆt, the equations in (8.39) assume the form given by the equations in (7.153). Accordingly, the pricing equation for the path-dependent model with the symmetric “peakon” payoff can be written as follows:

VtCC+12a0+a1yxVyyCC14VCC+κyxVxCC=0,VCCtˉ,x,y=ey/2.(8.40)

The Kelvin wave parameters are governed by the equations of the following form:

αtt,tˉa02γ2t,tˉ+14=0,   αtˉ,tˉ=0,iβtt,tˉ+a12γ2t,tˉ+14iκβt,tˉ=0,   βtˉ,tˉ=0,iγtt,tˉa12γ2t,tˉ+14+iκβt,tˉ=0,   γtˉ,tˉ=l.(8.41)
8.3.6 Bachelier Model

In the Bachelier model, the corresponding discounted t.p.d.f. is given by a modified (6.96):

ϖt,S,tˉ,Sˉ=12πΣ2t,tˉexpSˉFT22Σ2t,tˉ,(8.42)

where

Σ2t,tˉ=σˆ2e2rT12r.(8.43)

By virtue of (8.7), one can price European calls and puts as follows:

Vt,FT=erTϕFTKNϕFTKΣT+ΣTnFTKΣT,(8.44)

or, in spot terms:

Vt,S=ϕSerTKNϕSerTKΣ˜T+Σ˜TnSerTKΣ˜T,(8.45)

where

Σ˜2T=σˆ21e2rT2r.(8.46)

See Reference BachelierBachelier (1900), Reference Schachermayer and TeichmannSchachermayer & Teichmann (2008), and Reference TerakadoTerakado (2019) for further details.

8.4 Asian Options with Arithmetic and Geometric Averaging

The most basic path-dependent options are fixed strike Asian calls and puts, whose payoff depends on the underlying value averaged between the inception and maturity. Such options are popular for commodity and energy trading and in many other circumstances. The average Aˆt,tˉ on the interval t,tˉ can be defined in several ways. The simplest and, as a result, the most popular is an arithmetic average:

Aˆt,tˉ=1TttˉSˆsds.(8.47)

A less frequent, but technically easier to deal with, alternative is a geometric average:

Aˆt,tˉ=exp1TttˉlnSˆsds.(8.48)

The payoff of an Asian option with maturity tˉ and fixed strike K is

UAˉt,tˉ=maxϕAˉt,tˉK,0,(8.49)

as before, ϕ=1 for a call, and ϕ=1 for a put. For the floating strike, the payoff is

USˉtˉ,Aˉt,tˉ=maxϕSˉtˉkAˉt,tˉ,0,(8.50)

where the nondimensional parameter k is called weighting; typically, k=1.

Start with the Bachelier model. Equations for pricing Asian Options with an arithmetic average are as follows:

dAˆt=Sˆtdt,   Aˆt=0,dSˆt=rSˆtdt+σdWˆt,   Sˆt=S.(8.51)

Thus, (6.114) and (6.115) are applicable. All one needs is the marginal distribution for Aˉt,tˉ, which is Gaussian:

ϖAˉNR,Σ2,(8.52)

where

R=BrTS,   Σ2=σ2rB0T2BrT+B2rT.(8.53)

Consider the discounted payoff of the Asian call option (say):

Ut,Aˉ=AˉTK+.(8.54)

The corresponding calculation is straightforward:

Ut,S=erTTKAˉTKeAˉR22Σ22πΣ2dAˉ=erTΣTTKRΣηeη222πdηerTTKRTTKRΣeη222πdη=erTΣTnRTKΣerTTKRTNRTKΣ.(8.55)

Analytical pricing of Asian options with arithmetic averaging for the Black–Scholes model is notoriously tricky; see Reference Geman and EydelandGeman and Eydeland (1995), Reference Rogers and ShiRogers and Shi (1995), and Reference LiptonLipton (1999, Reference Lipton2001). At the same time, pricing Asian options with geometric averaging can be done quickly; see Reference Barucci, Polidoro and VespriBarrucci et al. (2001), Reference LiptonLipton (2001), and Reference Di Francesco and PascucciDi Francesco and Pascucci (2005) , and references therein. Such options can be priced using formula (6.45) derived in Section 6. An alternative approach based on the path integral method is discussed in Reference Devreese, Lemmens and TempereDevreese et al. (2010). Define

xˆt=ttyˆsds,   yˆt=lnSˆt.(8.56)

Then

dxˆt=yˆtdt,   xˆt=0,dyˆt=rσ22dt+σdWˆt,   yˆt=lnSˆty.(8.57)

The value of the option can be written as follows:

Ut,S=erTx*ϕϖxˉexpxˉTexplnKdxˉ,(8.58)

where

x*=TlnK.(8.59)

Since (8.57) is a special case of (6.74), one can use the equations in (6.81) to obtain the marginal distribution for xˉ, which is a Gaussian distribution of the form:

ϖxˉ=expxˉp22σx22πσx2,σx2=σ2T33,   p=lnST+12rσ22T2.(8.60)

Thus,

Ut,S=J1t,SJ2t,S,(8.61)

where

J1t,S=erTx*ϕexpxˉp22σx2+xˉT2πσx2dxˉ=ϕe12r+σ26TSNϕd+,J2t,S=erTx*ϕexpxˉp22σx2+lnK2πσx2dxˉ=ϕerTKNϕd,(8.62)

where

d±=lnS/K+12rσ26±σ23Tσ2T/3.(8.63)

Finally, one obtains a well-known formula for the price of a fixed strike Asian option with geometric averaging:

Ut,S=ϕe12r+σ26TSNϕd+erTKNϕd.(8.64)

Of course, a similar formula holds when r,σ are time-dependent. The derivation, although very simple, seems to be new.

8.5 Volatility and Variance Swaps and Swaptions

8.5.1 Volatility Swaps and Swaptions

Recall that the Stein–Stein stochastic volatility model assumes that the volatility is driven by an OU process; see Reference Stein and SteinStein and Stein (1991). One needs to find Green’s function associated with the following augmented SDEs:

dxˆt=yˆtdt,   xˆt=0,dyˆt=χVolκVolyˆtdt+εVoldWˆt,   yˆt=yVol,(8.65)

or, equivalently,

dxˆt=yˆtdt,   xˆt=0,dyˆt=κVolθVolyˆtdt+εVoldWˆt,   yˆt=yVol,(8.66)

which describe the evolution of the volatility σˆtyˆt and its integral xˆt; the equations of (8.65) are identical to the equations of (6.98).

It can be shown that the pair xˉ,yˉ has the bivariate Gaussian distribution with the covariance matrix H given by (6.113), and mean p,q given by (6.114):

pq=TθVol+BˉκVolTyVolθVolθVol+AκVolTyVolθVol.(8.67)

Since the marginal distribution of xˆt given by (6.115 ) is Gaussian, the fair strike of a volatility swap with maturity t is simply the expected value of xˆt/T:

VolSwap=θVol+yVolθVolBˉκVolTT.(8.68)

Here

VolSwap=χVolκVol=θVol=y=1T1/2.(8.69)

Of course, one can calculate the expected value of xˆt/T via more straightforward means. To this end, (8.68) can be derived directly by taking expectations of SDE (8.65). However, as we shall see in the following subsection, (6.115) for the marginal distribution ϖxt,yVol,tˉ,xˉ allows one to solve more interesting problems, such as calculating prices of bonds and bond options; see the discussion that follows.

Moreover, by using this equation, one can price volatility swaptions with payoffs of the form:

Utˉ,xˉ=maxϕxˉx*,0.(8.70)

The price Ut,yVol becomes:

Ut,yVol=erTϕx*ϕxˉx*ϖxt,yVol,tˉ,xˉdxˉ=erTϕ2πh0t,tˉx*ϕxˉx*expxˉp22h0dxˉ=erTϕpx*Nϕpx*h0+h0npx*h0.(8.71)

It is clear that formula (8.71) is a variant of the Bachelier formula (8.44).

8.5.2 Variance Swaps and Swaptions

In contrast to volatility, which, despite common misconceptions, can be negative, variance must be nonnegative since it is a square of a real-valued quantity. Accordingly, the easiest way to model it is by using the augmented Feller process with ϑ>0; see (7.99).

Using (7.127), one can immediately obtain the following expression for the fair value of a variance swap for the Feller process:

VarSwap=θVar+yVarθVarBˉκVarTT,(8.72)

where θVar=χVar/κVar. Here

VarSwap=θVar=yVar=1T.(8.73)

While formulas (8.68) and (8.72) look the same but deal with the volatility and variance, respectively, the corresponding parameters have different meanings.

Alternatively, one can use the degenerate augmented OU process, see the equations of (7.174). Averaging away stochastic terms, one gets the following formula for the fair price of the variance swap:

VarSwap=θVol2+yVol2θVol2BˉκVolTT.(8.74)

It is clear that Equations (8.72) and (8.74) provide different fair values for a variance swap, although these values asymptotically agree. This fact reflects the so-called model risk – by using different models, one gets different answers to the same question.

Equation (7.123) can be used to calculate the price of a variance swaption:

Ut,yVar=12πx*ϕϕxˉx*ϝtˉ,keikxˉdkdxˉ=12πϝtˉ,kϕx*ϕxˉx*eikxˉdxˉdk=12πlimϵ0ϝtˉ,keikx*ϵx*ϕeikϕϵxˉdxˉ=12πlimϵ0ϝtˉ,keikx*ikϕϵ2dk,(8.75)

where ϝtˉ,k is given by (7.124).

8.6 Automated Market Makers

Variance and volatility swaps had long occupied a specific niche within the financial product landscape. Recently, they experienced an unexpected surge in interest due to the influence of cryptocurrency trading. These swaps have proven effective in hedging impermanent loss, a phenomenon generated by automated market makers; see Reference Lipton and HardjonoLipton and Hardjono (2021), Reference Lipton and TreccaniLipton and Treccani (2021), Reference Lipton and SeppLipton and Sepp (2022), Reference Cartea, Drissi and MongaCartea et al. (2023), Reference Fukasawa, Maire and WunschFukasawa et. al (2023), and others. This section closely follows Reference Lipton and HardjonoLipton and Hardjono (2021).

Let us consider a smart contract (SC), called an automated market maker (AMM) designed to facilitate exchanges of two tokens, TN1 and TN2. The analytical formula for the price of the second token in terms of the first defines the nature of the contract. AMMs have gained significant traction in recent years. Initially, anyone can participate as a market maker and liquidity provider by depositing TN1 and TN2 simultaneously and in the correct ratio into the collateral pool. Subsequently, participants can withdraw one token from the pool by delivering the other token according to the rules established by the underlying SC. While AMMs excel in facilitating stablecoin swaps, they can easily accommodate the exchange of various tokens, such as swapping a stablecoin, say USDT, for ethereum (ETH).

The actual exchange rate is determined by rules that rely on prior agreement. The available options are the constant sum, constant product, and mixture rules. Sources including Reference Angeris, Kao, Chiang, Noyes and ChitraAngeris et al. (2019), Reference EgorovEgorov (2019), Reference Zhang, Chen and ParkZhang et al. (2018), Reference Lipton and HardjonoLipton and Hardjono (2021), Reference Lipton and SeppLipton and Sepp (2022), and references therein offer detailed coverage of AMMs and comprehensive insights into their mechanisms.

Assuming that initially tokens TN1, TN2 are equal in value, one can define a constant sum AMM:

X+Y=Σ0,   X0=Y0=N,   Σ0=2N.(8.76)

Here X,Y are the quantities of TN1, TN2 in the pool. Equation (8.76) yields

Y=Σ0X,   dYdX=1.(8.77)

As per (8.77), the pool reaches depletion at X=Σ0, as it becomes advantageous for an arbitrageur to increase X from N to 2N when TN2 surpasses TN1 in value. The marginal price of TN2 relative to TN1, as expressed in the second equation (8.77), remains consistent and equal to one. A constant price is optimal for a constant sum AMM, particularly when dealing with stablecoins like TN1 and TN2, whose prices fluctuate mildly around their equilibrium values. Depleting the pool is rational in scenarios where transaction fees are nonexistent, even with a minimal deviation from equilibrium. However, under more realistic conditions with nonzero transaction fees, arbitrage becomes profitable only if the deviation surpasses a certain threshold.

The constant product rule defines more intricate and, importantly, practical AMMs:

XY=Π0,   X0=Y0=N,   Π0=N2.(8.78)

It is clear that

Y=Π0X,   dYdX=Π0X2.(8.79)

Consequently, an arbitrageur is unable to deplete such a pool, allowing it to persist indefinitely. In this scenario, it becomes evident that the price of TN2 relative to TN1 is no longer steady; instead, it rises (or falls) as X decreases (or increases).

To make liquidity provision more attractive to potential market makers, one can generalize the constant sum and constant product rules. Expressions (8.76) and (8.78) representing these rules can be formulated as follows:

ΣΣ01=0,   X0=Y0=N,   Σ0=2N,Π0Π1=0,   X0=Y0=N,   Π0=N2.(8.80)

where Σ=X+Y, Π=XY are the current sum and product, respectively. These rules can be combined as follows:

Π0Π1+αΣΣ01=0,X0=Y0=N,   Σ0=2N,   Π0=N2.(8.81)

Here, α>0 is an adaptive parameter, characterizing the transition from the constant product to the constant sum rule. The product Π is in the denominator to avoid the possibility of exhausting the entire pool and ensuring that

YXX0,     XYY0.(8.82)

Certainly, when AMM liquidity providers are exposed to arbitragers, they face potential losses stemming from a decline in collateral value below its buy-and-hold threshold. In financial terms, an AMM liquidity provider is an option seller experiencing negative convexity, so that they must impose transaction fees to offset these losses. The losses incurred by AMMs are (somewhat misleadingly) termed “impermanent” because they tend to vanish under the assumption of mean reversion. However, the validity of the mean-reversion assumption in real-world scenarios can vary. Introducing variables x and y where X=Nx and Y=Ny, one can express the constant sum rule described by Equations (8.76) and (8.77) as follows:

x+y=2,   x0=y0=1,(8.83)
yx=2x,  dydx=1.(8.84)

In terms of x and y, the constant product rule given by Equations (8.78) and (8.79) can be written in the following form:

xy=1,    x0=y0=1,(8.85)
yx=1x,  dydx=1x2.(8.86)

Finally, the mixed-rule equations of (8.81) written in terms of x and y become

1xy1+αx+y21=0,   x0=y0=1.(8.87)

Straightforward algebra yields

yα=12α21α+αx+21α+αx2+8αx1/2,dyαdx=121+21α+αx4x221α+αx2+8αx1/2,d2yαdx2=12α+8x321α+αx2+8αx1/2α21α+αx4x221α+αx2+8αx3/2.(8.88)

Assume that the external exchange price S of TN2 expressed in terms of TN1 moves away from its equilibrium value S0=1. Let S>1. For the constant sum contract, an arbitrageur can choose a number x, 1<x2, and deliver x1 of TN1 tokens to the pool in exchange for getting x1 of TN2 tokens. The profit or loss (P&L) is given by

Ωx=S1x1.(8.89)

Since Ω is a linear function of x, it is rational to exhaust the entire pool by choosing the following optimal values x*,y*,Ω*:

x*=2,   y*=0,   Ω*=S1.(8.90)

Similarly, when S<1:

x*=0,   y*=2,   Ω*=S1.(8.91)

The arbitraged portfolio’s value is π*S, where

π*S=2,   S1,2S,   S<1.,(8.92)

while the buy-and-hold portfolio’s value is S+1. The difference ω has the form

ω=S+1π*S=S1.(8.93)

In the DeFi parlance, ω is termed as impermanent loss. However, this description can be misleading as the loss can swiftly become permanent when S moves away from its assumed “equilibrium” value of one. The percentage loss in the actual portfolio compared to the buy-and-hold portfolio is structured as follows:

λ=1S1S+1.(8.94)

A similar calculation can be performed for the constant product contract. When S deviates from one, an arbitrageur can choose a number x>1 and deliver x1 tokens TN1 to the pool, while taking 1y tokens TN2 from the pool, where y=1/x. The P&L has the form:

Ωx=S11xx1.(8.95)

The optimality condition has the form

Ωx=Sx21=0,(8.96)

so that the corresponding optimal values x*,y*,Ω* are

x*=S,   y*=1S,   Ω*=S12.(8.97)

Hence, a constant product collateral pool remains inexhaustible. Throughout each phase, the ideal quantities of TN1 and TN2 maintained in the portfolio are both S. As both tokens’ values within the portfolio must equate, the suggested optimal value of TN2 in terms of TN1 is S*=x*/y*=S. The value of the arbitrage-driven portfolio stands at π*=2S, whereas the value of the buy-and-hold portfolio amounts to S+1. The difference is given by

ω=S+12S=S12.(8.98)

The corresponding percentage loss is

λ=12SS+1=S12S+1.(8.99)

For the mixed-rule AMM, the arbitrageur’s profit for S>1 has the form

Ωx=S1yαxx1,(8.100)

with the optimum achieved at xα*,yα*,Ωα* of the form

yαxα*=1S,   yα*=yαxα*,   Ωα*=S1yα*xα*1,(8.101)

with the optimal xα* via the Newton–Raphson method starting with a suitable xα0:

xαn+1=xαnyαxαn+1Syαxαn.(8.102)

Here yα, yα are given by the equations of (8.88). Due to quadratic convergence of the Newton–Raphson method, ten iterations provide machine accuracy, so that one can set xα*=xα10. The value of the arbitraged portfolio is

π*=xα*+Syαxα*.(8.103)

Figure 14 shows the constant sum, constant product, and mixed-rule curves, along with the relative prices of TN2 in terms of TN1 and the associated impermanent losses. It demonstrates that deviations from the tokens’ equilibrium values result in losses for the market maker. Impermanent loss is relatively minor for the constant product rule, moderate for the mixed rule, and notably high for the constant sum rule. Even when the price S sways by a factor of five from its equilibrium, the impermanent loss within the constant product rule remains manageable, especially compared to the mixed rule.

Figure 14 The constant sum, constant product, and mixed-rule curves, along with the relative prices of TN2 in terms of TN1 and the associated impermanent losses; α=10. Author’s graphics.

One can use variance swaps to hedge impermanent loss. For brevity, consider the constant product rule. The corresponding impermanent loss, shown in Figure 14, is given by (8.98). It can be viewed as a payoff of a nonstandard European option. The hedging approach is straightforward – one approximates this payoff with payoffs of options, which can be priced explicitly. Specifically, one can use two such options: the log and entropy contracts. The corresponding payoffs are as follows:

ULCS=cLCS1lnS,(8.104)
UECS=cECSlnSS1.(8.105)

The prefactors cLC, cEC are chosen in such a way that the value of the impermanent loss (8.98) and the hypothetical payoffs (8.104) and (8.105) agree at the point S=1 up to the third derivative, so that

cLC=cEC=12.(8.106)

Assuming that S is driven by the geometric Brownian motion with stochastic volatility, one can find the value of the log and entropy contracts at time t at the point S=1, by solving the following problems:

Ut+12vS2USS+2ερSUSv+ε2Uvv+χκvUv=0,(8.107)

supplied with terminal conditions of the form

ULCtˉ,S,v=S1lnS,(8.108)

and

UECtˉ,S,v=SlnSS1,(8.109)

respectively.

The corresponding solutions are well-known and easy to find. One can present ULCt,S as follows:

ULCt,S,v,tˉ=ΦLCt,v,tˉ+S1lnS,(8.110)

where

ΦtLC+12v1+ε2ΦvvLC+χκvΦvLC=0,ΦLCtˉ,v,tˉ=0.(8.111)

Accordingly,

ΦLCt,v,tˉ=αLCt,tˉ+βLCt,tˉv,(8.112)

where

αtLCt,tˉ+χβLCt,tˉ=0,   αLCtˉ,tˉ=0,βtLCt,tˉκβLCt,tˉ+12=0,   βLCtˉ,tˉ=0.(8.113)

Thus,

αLCt,tˉ=χ2κTBˉκT,βLCt,tˉ=BˉκT2,(8.114)

so that

ULCt,S,v,tˉ=12χTκ+vχ2κBˉκT+S1lnS.(8.115)

It is clear that ULCt,1,v,tˉ is in agreement with (8.72).

One can calculate UECt,S,v in a similar fashion by representing it in the form:

UECt,S,v,tˉ=ΦECt,v,tˉS+SlnSS1,(8.116)

where, once the common factor S is omitted,

ΦtEC+12v1+2ερΦvEC+ε2ΦvvEC+χκvΦvEC=0,ΦLCtˉ,v,tˉ=0.(8.117)

As before,

ΦECt,v,tˉ=αECt,tˉ+βECt,tˉv,(8.118)

where

αtECt,tˉ+χβECt,tˉ=0,   αECtˉ,tˉ=0,βtECt,tˉκερβECt,tˉ+12=0,   βECtˉ,tˉ=0.(8.119)

Thus,

αECt,tˉ=χ2κερTBˉκερT,(8.120)
βECt,tˉ=BˉκερT2,UECt,S,v,tˉ=12χTκ1+vχ2κ1Bˉκ1TS+SlnSS1.(8.121)

where κ1=κερ. Equations (8.115) and (8.121) allow us to estimate the amount a liquidity provider needs to collect to cover the expected impermanent loss.

However, it turns out (which comes as a surprise, at least to the present author) that one can solve the pricing problem (8.107) with the exact terminal condition (8.98) explicitly, since the impermanent loss does not have any optionality and is a linear combination of the so-called power contracts with payoffs of the form S,S,1.Footnote 10

Thus, by using an appropriate Kelvin wave, one can solve the problem (8.107) with the power terminal condition:

UνS=Sν.(8.122)

Of course, for ν=0,1, the solution is trivial; for other values of ν, additional efforts are needed. To be concrete, it is assumed that 0<ν<1; for other values of ν, the solution can blow up in finite time. The price of the power contract with the payoff Sν (even when the interest rate r0) is given by a Kelvin wave:

Vt,S,tˉ=eαt,tˉ+βt,tˉvSν,(8.123)

where αt,βt solve the following system of ODEs:

αtt,tˉ+χβt,tˉ=0,   αtˉ,tˉ=0,βtt,tˉ+ε22β2t,tˉκνρεβt,tˉ+νν12=0,   βtˉ,tˉ=0,(8.124)

which has an explicit solution given by Equations (7.111)–(7.114) with

λ±2+κνρελ±+ε2νν14=0,λ±=μ±ζ,(8.125)
μ=κνρε2,   ζ=κνρε2ε2νν12.(8.126)

Thus, both μ and ζ are real. Accordingly, one can represent α and β as follows:

αT=2χε2κνρεT2+lnζcoshTμsinhTζ,βT=νν1sinhT2ζcoshTμsinhT.(8.127)

The exact impermanent loss and its approximations are shown in Figure 15. This figure shows that maxULC,UEC strictly dominates the exact solution UEX, but, as time of liquidity provision grows, the corresponding upper bound becomes inaccurate.

Figure 15 The exact impermanent loss and its approximations via log and enthropy contracts. The corresponding parameters are T=3, dt=0.01, χ=0.2, κ=2.0, ε=0.2, ρ=0.5, v=0.15. Author’s graphics.

The calculation of the mixed-rule impermanent loss and its approximations is left to the reader as a difficult exercise.

In P&L modeling for AMMs, the primary aim is to ensure that the liquidity provider makes a profit or, at least, does not incur a loss. This profit stems from transaction fees charged by the pool, which must exceed the impermanent loss caused by collateral value dropping below its buy-and-hold threshold. These fees must exceed the impermanent loss. An arbitrageur needs to add more tokens to the pool than the rule dictates to account for transaction fees. In the presence of nonzero transaction costs, the actual composition of the pool is time- and path-dependent. Given the stochastic nature of the log price, the analysis of P&L can only be conducted probabilistically through Monte Carlo simulations; see Reference Lipton and HardjonoLipton and Hardjono (2021) and Reference Lipton and SeppLipton and Sepp (2022). For the parameter selection used by these authors, automated liquidity provision is profitable on average. This profitability arises because the AMM accumulates more tokens by the process’s conclusion than initially possessed.

8.7 Bonds and Bond Options

8.7.1 Background

We now use the machinery developed in Sections 6 and 7 for pricing bonds and bond options in some popular fixed-income models, including Vasicek–Hull–White and Cox–Ingersoll–Ross.

8.7.2 Vasicek Model

One can use formulas derived in the previous subsection to price bonds and bond options in the popular Vasicek and Hull–White models; see Reference Hull and WhiteHull and White (1990); Reference VasicekVasicek (1977). Recall that Vasicek postulated the following dynamics for the short interest rate yˆt:

dyˆt=χκyˆtdt+εdWˆt,   yˆt=y,(8.128)

or, alternatively,

dyˆt=κθyˆtdt+εdWˆt,   yˆt=y,(8.129)

where κθ=χ.

At time t, the price of a bond maturing at time tˉ, which is denoted by Zt,y,tˉ, boils down to solving the following classical backward problem:

Ztt,y,tˉ+χκyZyt,y,tˉ+12ε2Zyyt,y,tˉyZt,y,tˉ=0,Ztˉ,y,tˉ=1.(8.130)

The standard affine ansatz yields

Zt,y,tˉ=expCBκyC=θε22κ2BκTε24κBκ2=BκTθ+h02,(8.131)

where h0 is given by (6.114).

One can use formulae derived in the previous section to come up with an alternative derivation. Introduce xˆt=ttyˆsds. The distribution of xˆt,yˆt is given by (6.45) with the covariance matrix H, given by (6.114) and the expected value r given by (6.115). Accordingly, the price of a bond can be written as follows:

Zt,y,tˉ=Eexˉ=12πh0exˉxˉp22h0dxˉ=ep+h02=expCTBκTy,(8.132)

so that Equations (8.131) and (8.132) are in agreement.

Knowing the joint Gaussian distribution for xˆt,yˆt, one can price an option on zero coupon bond maturing at time t˘>t, t˘t=T˘. The payoff of a European option with strike K has the form:

Utˉ,yˉ=maxϕexpCT˘BκT˘yˉexplnK,0.(8.133)

At maturity tˉ, the payoff is independent of xˉ; however, at inception it does depend on the realized value of xˉ. By using Equations (6.45), (6.114), and (6.115), one can write Ut,y (recall that here x=0) as follows:

Ut,y=J1t,yJ2t,y,(8.134)

where

J1t,y=12πdetH1/2ϕy*expΛxˉ,yˉxˉ+CT˘BκT˘yˉdxˉdyˉ,(8.135)
J2t,y=12πdetH1/2ϕy*expΛxˉ,yˉxˉ+lnKdxˉdyˉ,(8.136)
Λxˉ,yˉ=h2xˉp22h1xˉpyˉq+h0yˉq22detH,(8.137)

with hi given by (6.114), detH=h0h2h12. Here y* is defined as follows:

y*=CT˘lnKBκT˘.(8.138)

First, consider J1. Completing the square, one gets

Λxˉ,yˉxˉ+CT˘BκT˘yˉ=h2xˉpΞyh22Ξ2yˉ+h0yˉq22detHBκT˘yˉqp+CT˘BκT˘q,(8.139)

where

Ξyˉ=h1yˉqdetHh2.(8.140)

Integrating over xˉ, one obtains the following expression for J1:

J1t,y=ep+CT˘BκT˘q2πh2ϕy*expΞ2yˉ+h0yˉq2+2detHBκT˘yˉq2detHdyˉ.(8.141)

Completing the square one more time, one gets:

Ξ2+h0yˉq2+2detHBκT˘yˉq2detH=yˉq+h1+BκT˘h222h2+h02+BκT˘h1+Bκ2T˘h22,(8.142)

so that

J1t,y=ep+CT˘BκT˘q+h02+BκT˘h1+Bκ2T˘h222πh2t,tϕy*expyˉq+h1+BκT˘h222h2dyˉ=ϕep+CT˘BκT˘q+h02+BκT˘h1+Bκ2T˘h22Nϕy*q+h1+BκT˘h2h2.(8.143)

It is easy to see that Zt,y,t˘ is given by (8.143) with ϕ=1 and y*=, so that

Zt,y,t˘=ep+CT˘BκT˘q+h02+BκT˘h1+Bκ2T˘h22.(8.144)

Thus,

J1t,y=ϕZt,y,t˘NϕCT˘lnKBκT˘q+BκT˘h1+Bκ2T˘h2h2BκT˘.(8.145)

Direct verification of (8.143) is left to the reader as a useful exercise. By using this equation, it is easy but tedious to show that

J1t,y=ϕZt,y,t˘Nϕd+,d+=lnZt,y,t˘Zt,y,tˉKΣt,tˉ,t˘+Σt,tˉ,t˘2,(8.146)

where

Σt,tˉ,t˘=h2BκT˘.(8.147)

Second, consider J2, proceed in the same way as before, and represent J2t,y in the following form:

J2t,y=ϕZt,y,tˉNϕd,d=lnZt,y,t˘Zt,y,tˉKΣt,tˉ,t˘Σt,tˉ,t˘2.(8.148)

Finally, one arrives at the following familiar expression for the bond option price:

Ut,y=ϕZt,y,t˘Nϕd+Zt,y,tˉKNϕd.(8.149)
8.7.3 CIR Model

The CIR model postulates that the short rate follows the Feller process; see Reference Cox, Ingersoll Jr. and RossCox et al. (1985). Accordingly, the bond price can be calculated by using (7.123) with x=0, and k=i:

Zt,y,tˉ=ϖxt,y,tˉ,xˉexˉdxˉ=12πϝt,y,tˉ,keik1xˉdkdxˉ=ϝt,y,tˉ,i,(8.150)

where

ϝt,y,tˉ,i=exp2χμTε22χε2lnλE++λ+E2ζ+2λ+λE+Eyε2λE++λ+E,(8.151)

with

λ±=μ±ζ,μ=κ2,   ζ=κ2+2ε22.(8.152)

Thus,

Zt,y,tˉ=eC˜B˜y,(8.153)

where

C˜=χκTε22χε2lnλE++λ+E2ζ,B˜=E+EλE++λ+E,(8.154)

which coincides with the standard expressions given by Reference Cox, Ingersoll Jr. and RossCox et al. (1985).

8.8 European Options with Stochastic Interest Rates

This section shows how to price equity options with stochastic interest rates. While the formulation of this problem may appear straightforward, its solution proves to be tedious. It is assumed that interest rate is governed by the Ornstein–Uhlenbeck–Vasicek processes.

dyˆt=χκyˆtdt+εdZˆt,   yˆt=y,(8.155)

where Zˆt is the standard Wiener processes. The risk-neutral evolution of the foreign exchange is governed by the following equation:

dSˆtSˆt=yˆtdt+σdWˆt,   Sˆt=S,(8.156)

or, equivalently,

dxˆt=yˆt12σ2dt+σdWˆt,   xˆt=x,(8.157)

where xˆ=lnSˆ/K. In general, dZˆt and dWˆt are correlated, so that dZˆtdWˆt=ρdt.

Consider the familiar backward Kolmogorov problem for European calls and puts:

Ut+12ε2Urr+ρεσUrx+12σ2Uxx+χκrUy+y12σ2UxrU=0,Utˉ,y,x=Kϕex1+.(8.158)

As usual, start with the change of the dependent variable:

U=KB1V.(8.159)

where B=expα1β1y is the domestic bond price, given by (8.131), so that

Bt+12ε2Brr+χκrByrB=0.(8.160)

Hence,

Vt+12ε2Vrr+ρεσVrx+12σ2Vxx+x12ε2β1κrVy+y12σ2ρεσβ1Vx=0,Vtˉ,y,x=ϕex1+.(8.161)

Now, change independent variables t,y,xt,η1,η2, where

η1=y,   η2=α1+β1y+x.(8.162)

Thus,

t=t+α1+β1η1η2,y=η1+β1η2,   x=η2.(8.163)

so that

Vt+α1+β1η1Vη2+12ε2Vη1η1+2β1Vη1η2+β12Vη2η2+ρεσVη1η2+β1Vη2η2+12σ2Vη2η2+χ12ε2β1κη1Vη1+β1Vη2+η112σ2ρεσβ1Vη2=0,Vtˉ,η1,η2=ϕeη21+.(8.164)

Assume that Vt,η1,η2 only depends on t,η2, Vt,η1,η2=Vt,η2, which is consistent with the terminal condition. Thus,

Vt+12ε2β12+ρεσβ1+12σ2Vη2η2+α1+β1η1+b112ε2β1κη1β112σ2ρεσβ1+η1Vη2=0Vtˉ,η1,η2,η2=ϕeη21+.(8.165)

But

α1β1η1+12ε2β12χκη1β1η1=0,(8.166)

so that

Vt+12ε2β12+ρεσβ1+12σ2Vη2η2Vη2=0,Vtˉ,η1,η2=ϕeη21+.(8.167)

This is the classical Black–Scholes problem with time-dependent volatility:

Vt+12Σ2Vη2η2Vη2=0,Vtˉ,η2=ϕeη21+,(8.168)

where

Σ2=ε2Bκ2+2ρεσBκ+σ2.(8.169)

Thus, the price is

U=B1UC,PB2SB1;T,K,Σ2dsT,(8.170)

where UC,P are given by (8.20).

A similar technique can be used for the Heston model and the Stein–Stein model with stochastic interest rates. However, there is one significant difference between these two models - the former model works only when volatility and rate innovations are uncorrelated, while the latter model can handle arbitrary correlations.

9 Conclusions

Due to the space constraints, the discussion must be concluded here. It is left to the reader to explore further the application of mathematical tools and techniques based on Kelvin waves in financial engineering. Three particularly compelling problems are

  • the pricing and risk management of credit derivatives;

  • the exploration of mean-reverting trading strategies, such as pairs trading;

  • the examination of affine jump-diffusion and pseudo-differential processes.

References such as Reference Lipton and SheltonLipton and Shelton (2012), Reference Lipton and Lopez de PradoLipton and Lopez de Prado (2020), and others provide additional insights into these problems.

This Element has established a unified methodology for determining t.p.d.fs and expectations for affine processes through integral representations based on Kelvin waves. This approach has bridged various disciplines, uncovering profound connections between hydrodynamics, molecular physics, stochastic processes, and financial engineering. Both degenerate problems, which possess more independent variables than sources of uncertainty, and their nondegenerate counterparts are covered, showcasing the versatility of the method.

A surprising link is established between the Langevin equation for underdamped Brownian motion and the vorticity equation for two-dimensional flows in viscous incompressible fluids. Utilizing Kelvin wave expansions, the book solves several relevant financial problems, including the deriving convenient formulas for t.p.d.fs and expectations for processes with stochastic volatility, developing an analytically solvable model for path-dependent volatility, pricing of Asian options with geometric averaging, and pricing bonds and bond options by augmenting the short-rate process with its integral process.

The methodology introduced in this book can address a wide spectrum of complex problems, significantly enhancing the comprehension and modeling of stochastic systems across diverse fields.

Acknowledgments

I am grateful to my ADIA colleagues Majed Alromaithi, Marcos Lopez de Prado, Koushik Balasubramanian, Andrey Itkin, Oleksiy Kondratiev, Arthur Maghakian, Dmitry Muravey, Adil Reghai, other Q-team colleagues, my ADIA Lab colleague Horst Simon, and a former Bank of America colleague, Artur Sepp, for their encouragement and council. The kind invitation by Riccardo Rebonato to contribute to Cambridge Elements in Quantitative Finance is much appreciated. I am grateful to Drs. Nicola Ghazi and Piergiorgio Neri from Cleveland Clinic Abu Dhabi for saving the vision in my left eye, thus allowing me to finish this Element. Last but not least, the help of my wife, Marsha Lipton, especially her editorial suggestions and financial insights, has been critical in producing this Element.

Alexander Lipton is a Global Head of Research & Development at Abu Dhabi Investment Authority, an Advisory Board member at ADIA Lab, a Professor of Practice at Khalifa University, and a Connection Science Fellow at MIT. He is a Co-Founder of Sila, a company providing digital wallet & ACH payment services, and an advisory board member at several companies worldwide. From 2006 to 2016, Alexander was Co-Head of the Global Quantitative Group and Quantitative Solutions Executive at Bank of America. Before that, he held senior managerial positions at several leading financial institutions. Additionally, Alexander held visiting professorships at EPFL, NYU, Oxford, and Imperial College. Earlier, Alexander was a Full Professor at the University of Illinois and a Consultant at the Los Alamos National Laboratory. Risk Magazine awarded him the Inaugural Quant of the Year Award in 2000 and the Buy-side Quant of the Year Award in 2021. Alexander has authored/edited thirteen books and over a hundred scientific papers on nuclear fusion, astrophysics, applied mathematics, financial engineering, distributed ledgers, and quantum computing. He holds several US patents.

  • Riccardo Rebonato

  • EDHEC Business School

  • Editor Riccardo Rebonato is Professor of Finance at EDHEC Business School and holds the PIMCO Research Chair for the EDHEC Risk Institute. He has previously held academic positions at Imperial College, London, and Oxford University and has been Global Head of Fixed Income and FX Analytics at PIMCO, and Head of Research, Risk Management and Derivatives Trading at several major international banks. He has previously been on the Board of Directors for ISDA and GARP, and he is currently on the Board of the Nine Dot Prize. He is the author of several books and articles in finance and risk management, including Bond Pricing and Yield Curve Modelling (2017, Cambridge University Press).

About the series

  • Cambridge Elements in Quantitative Finance aims for broad coverage of all major topics within the field. Written at a level appropriate for advanced undergraduate or graduate students and practitioners, Elements combines reports on original research covering an author’s personal area of expertise, tutorials and masterclasses on emerging methodologies, and reviews of the most important literature.

Footnotes

1 Data aequatione quotcunque fluentes quantitae involvente fluxiones invenire et vice versa. – “Given an equation involving any number of fluent quantities to find the fluxions, and vice versa.” Vladimir Arnold paraphrased the statement as follows: “It is useful to solve differential equations.” Cambridge University Library, Department of Manuscripts and University Archives. ItemReference Code: GBR/0012/MS Add.9597/2/18/56.

2 Newton’s anagram is an early example of a one-way hash function. An anagram is easy to calculate, provided the message is known, but not vice versa. Hash functions are indispensable in modern cryptography, including its applications to cryptocurrencies such as Bitcoin; see, for example, Reference Lipton and TreccaniLipton and Treccani (2021).

3 Thus, even the greatest minds occasionally can be myopic – it took eighty years for fluid dynamists to connect the dots and observe that v˜,p˜/ρ solve the nonlinear Euler equations.

4 Kolmogorov published his seminal papers in German where his name appears under the transliteration of Kolmogoroff, the spelling used in the original articles has been retained.

5 A Feller process might or might not be able to reach zero, which depends on the magnitude of the ratio 2χ/ε2.

6 We are grateful to Andrey Itkin for pointing this out.

7 In a special case a=1, b=0, ξ=0, θ=0, Reference He, Chen, Fang and HeHe et al. (2021) attempted to solve the problem considered previously. However, the authors made a severe error in transitioning from (2.2) to (2.3). In contrast to Kolmogorov’s minor error, their error cannot be repaired. Dimensional analysis shows that the proposed solution is completely incorrect. In our notation, it has the following form:

ϖxˉ,yˉdxˉdyˉ=723π3T61+4yˉ24T+3xˉ+yˉT/22T37/2dxˉdyˉ.

Introducing rescaled variables, ζ=xˉT−3/2, η=yˉT−1/2, we get

ϖζ,ηdζdη=723π3T41+4η24+3ζ+η227/2dζdη.

This expression explicitly depends on T, which is impossible, since its integral must be equal to unity.

8 Of course, it is not surprising that such explosions exist since Riccati equations are well-known to have solutions exploding in finite time. Consider the following Riccati initial-value problem

f′+af2+bf+c=0,   fτ=d,

and assume that b2−4ac<0. The solution of this initial-value problem has the form

fτ,t=−ζatanζt−τ−arctan2ad+b2ζ−b2a,

where ζ=4ac−b22. The corresponding blow-up time t* has the form

t*=τ+π2ζ+1ζarctan2ad+b2ζ.

9 Despite the fact that these formulas were originally derived by Reference LiptonLipton (2001), they are frequently mistakenly attributed to Reference Dragulescu and YakovenkoDragulescu and Yakovenko (2002).

10 This fact is true even for nonzero interest rates.

References

Abramowitz, M. & Stegun, I. A., eds. 1964 Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables. Washington, DC: US Government Printing Office.Google Scholar
Aksenov, A. V. 1995 Symmetries of linear partial differential equations and fundamental solutions. Doklady Mathematics 51(3), 329331.Google Scholar
Andersen, L. B. G. & Piterbarg, V. 2007 Moment explosions in stochastic volatility models. Finance and Stochastics 11, 29100.CrossRefGoogle Scholar
Angeris, G., Kao, H. T., Chiang, R., Noyes, C. & Chitra, T. 2019 An analysis of Uniswap markets. Cryptoeconomic Systems Journal. 1(1). DOI: 10.21428/2F58320208.C9738E64Google Scholar
Arnold, L. 1974 Stochastic Equations: Theory and Applications. New York: John Wiley.Google Scholar
Bachelier, L. 1900 Théorie de la spéculation. Annales de l’Ecole Normale Supérieure 17, 2186.CrossRefGoogle Scholar
Barucci, E., Polidoro, S. & Vespri, V. 2001 Some results on partial differential equations and Asian options. Mathematical Models and Methods in Applied Sciences 11(3), 475497.CrossRefGoogle Scholar
Bayly, B. J. 1986 Three-dimensional instability of elliptical flow. Physical Review Letters 57, 21602163.CrossRefGoogle ScholarPubMed
Bayly, B. J., Holm, D. D. & Lifschitz, A., 1996. Three-dimensional stability of elliptical vortex columns in external strain flows. Philosophical Transactions of the Royal Society A 354(1709), 895926.Google Scholar
Berest, Yu. Yu. 1993 Group analysis of linear differential equations in distributions and the construction of fundamental solutions. Differential Equations 29(11), 17001711.Google Scholar
Bergomi, L. 2015 Stochastic Volatility Modeling. Boca Raton, FL: CRC Press.CrossRefGoogle Scholar
Bharucha-Reid, A. T. 1960 Elements of the Theory of Markov Processes and Their Applications. New York: McGraw-Hill.Google Scholar
Bick, A. & Reisman, H. 1993 Generalized implied volatility. Preprint. https://www.researchgate.net/publication/228601694_Generalized_implied_volatilityGoogle Scholar
Black, F. 1976 The pricing of commodity contracts. Journal of Financial Economics 81(3), 167179.CrossRefGoogle Scholar
Black, F. & Scholes, M. 1973 The pricing of options and corporate liabilities. Journal of Political Economy 81(3), 637659.CrossRefGoogle Scholar
Bluman, G. & Kumei, S. 1989 Symmetries and Differential Equations. Berlin: Springer.CrossRefGoogle Scholar
Boness, A. J. 1964 Elements of a theory of a stock option value. Journal of Political Economy 72(2), 163175.CrossRefGoogle Scholar
Boyarchenko, S. & Levendorskii, S. 2002 Non-Gaussian Merton–Black–Scholes Theory. River Edge, NJ: World Scientific.CrossRefGoogle Scholar
Carr, P., Lipton, A. & Madan, D. 2002 The reduction method for valuing derivative securities. Working Paper, New York University. https://www.researchgate.net/publication/2860358_The_Reduction_Method_for_Valuing_Derivative_SecuritiesGoogle Scholar
Cartea, Á. , Drissi, F. & Monga, M. 2023 Predictable losses of liquidity provision in constant function markets and concentrated liquidity markets (August 15). Applied Mathematical Finance. Available at SSRN 4541034: https://ssrn.com/abstract=4541034 or http://dx.doi.org/10.2139/ssrn.4541034CrossRefGoogle Scholar
Chandrasekhar, S. 1943 Stochastic problems in physics and astronomy. Reviews of Modern Physics 15(1), 189.CrossRefGoogle Scholar
Chandrasekhar, S. 1961 Hydrodynamic and Hydromagnetic Stability. Oxford: Clarendon.Google Scholar
Chapman, S. 1928 On the Brownian displacements and thermal diffusion of grains suspended in a nonuniform fluid. Proceedings of the Royal Society of London. Series A 119(781), 3454.Google Scholar
Cordes, H. O. 1995 The Technique of Pseudodifferential Operators. Cambridge: Cambridge University Press.CrossRefGoogle Scholar
Cox, J. C., Ingersoll Jr., J. E. & Ross, S. A. 1985 A theory of the term structure of interest rates. Econometrica 53(2), 385408.CrossRefGoogle Scholar
Craddock, M. 2012 Lie symmetry methods for multi-dimensional parabolic PDEs and diffusions. Journal of Differential Equations 252(1), 5690.CrossRefGoogle Scholar
Craddock, M. & Platen, E. 2004 Symmetry group methods for fundamental solutions. Journal of Differential Equations 207(2), 285302.CrossRefGoogle Scholar
Craik, A. D. D. & Criminale, W. O. 1986 Evolution of wavelike disturbances in shear flows: A class of exact solutions of the Navier–Stokes equations. Proceedings of the Royal Society of London. Series A 406, 1326.Google Scholar
Dai, Q. & Singleton, K. J. 2000 Specification analysis of affine term structure models. Journal of Finance 55(5), 19431978.CrossRefGoogle Scholar
Davis, M. H. A. 2004 Complete-market models of stochastic volatility. Proceedings of the Royal Society of London. Series A 460, 1126.CrossRefGoogle Scholar
Derman, E. & Kani, I. 1994 Riding on a smile. Risk Magazine 7(2), 3239.Google Scholar
Devreese, J. P. A., Lemmens, D. & Tempere, J. 2010 Path integral approach to Asian options in the Black–Scholes model. Physica A 389(4), 780788.CrossRefGoogle Scholar
Di Francesco, M. & Pascucci, A. 2004 On the complete model with stochastic volatility by Hobson and Rogers. Proceedings of the Royal Society of London. Series A 460(2051), 33273338.CrossRefGoogle Scholar
Di Francesco, M. & Pascucci, A. 2005 On a class of degenerate parabolic equations of Kolmogorov type. Applied Mathematics Research eXpress 3, 77116.CrossRefGoogle Scholar
Dragulescu, A. A. & Yakovenko, V. M. 2002 Probability distribution of returns in the Heston model with stochastic volatility. Quantitative Finance 2(6), 443453.CrossRefGoogle Scholar
Duffie, D., Filipovic, D. & Schachermayer, W. 2003 Affine processes and applications in finance. Annals of Applied Probability 13(3), 98410053.CrossRefGoogle Scholar
Duffie, J. D. & Kan, R. 1996 A yield-factor model of interest rates. Mathematical Finance 6(4), 379406.CrossRefGoogle Scholar
Duffie, D., Pan, J. & Singleton, K. 2000 Transform analysis and asset pricing for affine jump-diffusions. Econometrica 68(6), 13431376.CrossRefGoogle Scholar
Duong, M. H. & Tran, H. M. 2018 On the fundamental solution and a variational formulation for a degenerate diffusion of Kolmogorov type. Discrete & Continuous Dynamical Systems: Series A 38(7), 34073438.CrossRefGoogle Scholar
Dupire, B. 1994 Pricing with a smile. Risk Magazine 7(1), 1820.Google Scholar
Ebeling, W., Gudowska-Nowak, E. & Sokolov, I. M. 2008 On stochastic dynamics in physics: Remarks on history and terminology. Acta Physica Polonica B 39(5), 100310017.Google Scholar
Egorov, M. 2019 StableSwap: Efficient mechanism for Stablecoin liquidity. White paper.Google Scholar
Fabijonas, B., Holm, D. D. & Lifschitz, A. 1997 Secondary instabilities of flows with elliptic streamlines. Physical Review Letters 78(10), 19001903.CrossRefGoogle Scholar
Feller, W. 1951 Two singular diffusion problems. Annals of Mathematics 54(1), 173182.CrossRefGoogle Scholar
Feller, W. 1952 The parabolic differential equations and the associated semi-groups of transformations. Annals of Mathematics 55, 468518.CrossRefGoogle Scholar
Feller, W. 1971 An Introduction to Probability Theory and Its Application, vol. 2, 2nd ed. New York: John Wiley.Google Scholar
Filipovic, D. 2009 Term-structure models. Berlin: Springer.CrossRefGoogle Scholar
Fokker, A. D. 1914 Die mittlere Energie rotierender elektrischer Dipole im Strahlungsfeld. Annalen der Physik 348(5), 810820.CrossRefGoogle Scholar
Fourier, J. B. 1822 Théorie Analytique de la Chaleur. Paris: Firmin Didot, Père et Fils.Google Scholar
Friedlander, S. & Lipton-Lifschitz, A. 2003 Localized instabilities in fluids. In Handbook of Mathematical Fluid Dynamics 2, 289354. North-Holland: Amsterdam.CrossRefGoogle Scholar
Friedlander, S. & Vishik, M. 1991 Instability criteria for the flow of an inviscid incompressible fluid. Physical Review Letters 66, 22042206.CrossRefGoogle ScholarPubMed
Friz, P. & Keller-Ressel, M. 2010. Moment explosions in stochastic volatility models. Encyclopedia of Quantitative Finance, 12471253.Google Scholar
Fukasawa, M., Maire, B. & Wunsch, M. 2023 Weighted variance swaps hedge against impermanent loss. Quantitative Finance 23(6), 111.CrossRefGoogle Scholar
Gatheral, J., Jaisson, T. & Rosenbaum, M. 2018 Volatility is rough. Quantitative Finance 18(6), 933949.CrossRefGoogle Scholar
Geman, H. & Eydeland, A. 1995 Asian options revisited: Inverting the Laplace transform. Risk Magazine 8(4), 6567.Google Scholar
Gershon, D., Lipton, A., Rosenbaum, M. & Wiener, Z. eds. 2022 Options-45 Years Since the Publication of the Black–Scholes–Merton Model: The Gershon Fintech Center Conference. Singapore: World Scientific.Google Scholar
Gihman, I. I. & Skorohod, A. V. 1972 Stochastic Differential Equations. New York: Springer.CrossRefGoogle Scholar
Giorno, V. & Nobile, A. G. 2021 Time-inhomogeneous Feller-type diffusion process in population dynamics. Mathematics 9(16), 1879.CrossRefGoogle Scholar
Guyon, J. 2014 Path-dependent volatility. Risk Magazine 27(10).Google Scholar
Hagan, P., Kumar, D., Lesniewski, A. & Woodward, D. 2002 Managing smile risk. Wilmott Magazine 9, 841008.Google Scholar
Hänggi, P., Talkner, P. & Borkovec, M. 1990 Reaction-rate theory: Fifty years after Kramers. Reviews of Modern Physics 62(2), 251341.CrossRefGoogle Scholar
Hanson, F. B. 2007. Applied Stochastic Processes and Control for Jump-Diffusions: Modeling, Analysis and Computation. Philadelphia, PA: Society for Industrial and Applied Mathematics.CrossRefGoogle Scholar
He, C., Chen, J., Fang, H. & He, H. 2021 Fundamental solution of fractional Kolmogorov–Fokker–Planck equation. Examples and Counterexamples 1, 100031.CrossRefGoogle Scholar
Heston, S. L. 1993 A closed-form solution for options with stochastic volatility with applications to bond and currency options. Review of Financial Studies 6, 327343.CrossRefGoogle Scholar
Hobson, D. G. & Rogers, L. C. G. 1998 Complete models with stochastic volatility. Mathematical Finance 8(1), 2748.CrossRefGoogle Scholar
Hörmander, L. 1967 Hypoelliptic second order differential equations. Acta Mathematica 119, 147171.CrossRefGoogle Scholar
Hull, J. & White, A. 1990 Pricing interest rate derivative securities. Review of Financial Studies 3, 573592.CrossRefGoogle Scholar
Ibragimov, N. H. 1985 Transformation Groups Applied to Mathematical Physics. Dordrecht: D. Reidel.CrossRefGoogle Scholar
Ivasishen, S. D. & Medynsky, I. P. 2010 The Fokker–Planck–Kolmogorov equations for some degenerate diffusion processes. Theory of Stochastic Processes 16(1), 5766.Google Scholar
Jacob, N. & Schilling, R. L. 2001 Lévy-type processes and pseudodifferential operators. In Lévy Processes: Theory and Applications, Barndorff-Nielsen, O. E. et al. (eds.), (pp. 139168). Boston, MA: Birkhäuser.CrossRefGoogle Scholar
Janek, A., Kluge, T., Weron, R. & Wystup, U. 2011 FX smile in the Heston model. In Statistical Tools for Finance and Insurance (pp. 133162). Berlin, Heidelberg: Springer.CrossRefGoogle Scholar
Jex, M., Henderson, R. & Wang, D. 1999 Pricing exotics under the smile. Risk Magazine 12(11), 7275.Google Scholar
Kelvin, Lord 1887 Stability of fluid motion: Rectilinear motion of viscous fluid between two parallel planes. Philosophical Magazine 24, 188196.Google Scholar
Klein, O. 1921 Zur statistischen Theorie der Suspension und Lösungen. Inaugural-Dissertation. Uppsala: Almqvist & Wiksells.Google Scholar
Kolmogoroff, A. 1931 Über die analytischen Methoden in der Wahrscheinlichkeitsrechnung. Mathematische Annalen 104(1), 415458.CrossRefGoogle Scholar
Kolmogoroff, A. 1933 Zur Theorie der stetigen zufälligen Prozesse. Mathematische Annalen 108, 149160CrossRefGoogle Scholar
Kolmogoroff, A. 1934 Zufallige Bewegungen (Zur Theorie der Brownschen Bewegung), Annals of Mathematics 35(1), 116117.CrossRefGoogle Scholar
Kovalenko, S., Stogniy, V. & Tertychnyi, M. 2014 Lie symmetries of fundamental solutions of one (2+ 1)-dimensional ultra-parabolic Fokker–Planck–Kolmogorov equation. arXiv preprint arXiv:1408.0166.Google Scholar
Kramers, H. A. 1940 Brownian motion in a field of force and the diffusion model of chemical reactions. Physica 7(4), 284304.CrossRefGoogle Scholar
Kuptsov, L. P. 1972 The fundamental solutions of a certain class of elliptic-parabolic second order equations. Differential Equations 8, 16491660.Google Scholar
Lanconelli, E., Pascucci, A. & Polidoro, S. 2002 Linear and nonlinear ultraparabolic equations of Kolmogorov type arising in diffusion theory and in finance. In Sh, Michael. Birman, Stefan Hildebrandt, Solonnikov, Vsevolod A., and Uraltseva, Nina N., eds., Nonlinear Problems in Mathematical Physics and Related Topics, vol. II, 243265. New York: Kluwer/Plenum.CrossRefGoogle Scholar
Langevin, P. 1908 Sur la théorie du mouvement brownien. Comptes rendus de l’Académie des Sciences 146, 530533.Google Scholar
Lewis, A. 2000 Option Valuation under Stochastic Volatility with Mathematica Code. Newport Beach, CA: Finance Press.Google Scholar
Lewis, A. 2001 A simple option formula for general jump-diffusion and other exponential Lévy processes. Available at SSRN: https://ssrn.com/abstract=282110 or http://dx.doi.org/10.2139/ssrn.282110.CrossRefGoogle Scholar
Lifschitz, A. 1991 Short wavelength instabilities of incompressible three-dimensional flows and generation of vorticity. Physics Letters A 157(8–9), 481487.CrossRefGoogle Scholar
Lifschitz, A. 1995 Exact description of the spectrum of elliptical vortices in hydrodynamics and magnetohydrodynamics. Physics of Fluids A 7, 16261636.CrossRefGoogle Scholar
Lifschitz, A. & Hameiri, E. 1991a A universal instability in fluid dynamics. In 1991 International Sherwood Fusion Theory Conference.Google Scholar
Lifschitz, A. & Hameiri, E. 1991b Local stability conditions in fluid dynamics. Physics of Fluids A 3, 26442651.CrossRefGoogle Scholar
Lipton, A. 1999 Similarities via self-similarities. Risk Magazine 12(9), 1011005.Google Scholar
Lipton, A. 2000 Pricing and Risk-Managing Exotics on Assets with Stochastic Volatility. Presentation, Risk Minds, Geneva.Google Scholar
Lipton, A. 2001 Mathematical Methods for Foreign Exchange: A Financial Engineer’s Approach. Singapore: World Scientific.CrossRefGoogle Scholar
Lipton, A. 2002 The volatility smile problem. Risk Magazine 15(2), 6165.Google Scholar
Lipton, A. 2018 Financial Engineering: Selected Works of Alexander Lipton. Singapore: World Scientific.CrossRefGoogle Scholar
Lipton, A. 2023 Kelvin Waves, Klein–Kramers and Kolmogorov Equations, Path-Dependent Financial Instruments: Survey and New Results. arXiv preprint arXiv:2309.04547.Google Scholar
Lipton, A., Gal, A. & Lasis, A. 2014 Pricing of vanilla and first-generation exotic options in the local stochastic volatility framework: Survey and new results. Quantitative Finance 14(11), 18991922.CrossRefGoogle Scholar
Lipton, A. & Hardjono, T. 2021 Blockchain intra- and interoperability. In Innovative Technology at the Interface of Finance and Operations: Volume II (pp. 130). Cham: Springer International Publishing.Google Scholar
Lipton, A. & Lopez de Prado, M. 2020 A closed-form solution for optimal Ornstein–Uhlenbeck driven trading strategies. International Journal of Theoretical and Applied Finance 23(08), 2050056.CrossRefGoogle Scholar
Lipton, A. & Reghai, A. 2023 SPX, VIX and scale-invariant LSV. Wilmott Magazine 2023(126), 7884.CrossRefGoogle Scholar
Lipton, A. & Sepp, A. 2008 Stochastic volatility models and Kelvin waves. Journal of Physics A: Mathematical and Theoretical 41(34), 344012 (23pp).CrossRefGoogle Scholar
Lipton, A. & Sepp, A. 2022 Automated market-making for fiat currencies. Risk Magazine 35(5). https://www.risk.net/cutting-edge/investments/7948106/automated-market-making-for-fiat-currencies.Google Scholar
Lipton, A. & Shelton, D. 2012 Credit default swaps with and without counterparty and collateral adjustments. Stochastics: An International Journal of Probability and Stochastic Processes 84(5–6), 603624.CrossRefGoogle Scholar
Lipton, A. & Treccani, A. 2021 Blockchain and Distributed Ledgers: Mathematics, Technology, and Economics. Singapore: World Scientific.CrossRefGoogle Scholar
Lipton-Lifschitz, A. 1999 Predictability and unpredictability in financial markets. Physica D: Nonlinear Phenomena 133(1–4), 321347.CrossRefGoogle Scholar
Masoliver, J. 2016 Nonstationary Feller process with time-varying coefficients. Physical Review E 93(1), 012122.CrossRefGoogle ScholarPubMed
Merton, R. C. 1973 Theory of rational option pricing. Bell Journal of Economics and Management Science 4(1), 141183.Google Scholar
Merton, R. C. 1976 Option pricing when underlying stock returns are discontinuous. Journal of Financial Economics 3(1–2), 125144.CrossRefGoogle Scholar
Morse, P. M. & Feshbach, H. 1953 Methods of Theoretical Physics, Part I. New York: McGraw-Hill.Google Scholar
Olver, P. J. 1986 Applications of Lie Groups to Differential Equations, 1st ed. New York: Springer.CrossRefGoogle Scholar
Orr, W. McF. 1907 The stability or instability of the steady motions of a perfect fluid. Proceedings of the Royal Irish Academy A 27, 969.Google Scholar
Ovsiannikov, L. V. 1982 Group Analysis of Differential Equations. New York: Academic Press.Google Scholar
Pascucci, A. 2005 Kolmogorov equations in physics and in finance. In , Elliptic and Parabolic Problems, C. Bandle, et al., eds. Progress in Nonlinear Differential Equations and Their Applications, 63, 313324. Basel: Birkhäuser.Google Scholar
Piessens, R. 2000 The Hankel Transform. In Poularikas, A. D., ed., The Transforms and Applications Handbook: Second Edition. Boca Raton, FL: CRC Press.Google Scholar
Planck, M. 1917 Über einen Satz der statistischen Dynamik und seine Erweiterung in der Quantentheorie. Sitzungsberichte der Königlich Preussischen Akademie der Wissenschaften, 324341.Google Scholar
Reghai, A. 2015 Quantitative Finance: Back to Basics. New York: Palgrave MacMillan.Google Scholar
Risken, H. 1989 The Fokker–Planck Equation: Method of Solution and Applications. New York: Springer.Google Scholar
Rogers, L. C. G. & Shi, Z. 1995 The value of an Asian option. Journal of Applied Probability 32(4), 107710088.CrossRefGoogle Scholar
Rubinstein, M. 1994 Implied binomial trees. Journal of Finance 49(3), 771818.CrossRefGoogle Scholar
Samuelson, P. A. 1965 Rational theory of warrant pricing. Industrial Management Review 6, 1332.Google Scholar
Schachermayer, W. & Teichmann, J. 2008 How close are the option pricing formulas of Bachelier and Black–Merton–Scholes? Mathematical Finance: An International Journal of Mathematics, Statistics and Financial Economics 18(1), 155170.CrossRefGoogle Scholar
Schmelzle, M. 2010 Option pricing formulae using Fourier transform: Theory and application. Preprint, https://pfadintegral.com/articles/option-pricing-formulae-using-fourier-transform.Google Scholar
Schöbel, R. & Zhu, J. 1999 Stochastic volatility with an Ornstein–Uhlenbeck process: An extension. Review of Finance 3(1), 2346.CrossRefGoogle Scholar
Sepp, A., 2007. Affine Models in Mathematical Finance: An Analytical Approach. Tartu: University of Tartu Press.Google Scholar
Stein, E. M. & Stein, J. C. 1991 Stock price distributions with stochastic volatility: An analytic approach. Review of Financial Studies 4, 727752.CrossRefGoogle Scholar
Terakado, S., 2019 On the option pricing formula based on the Bachelier model. Available at SSRN 3428994.CrossRefGoogle Scholar
Uhlenbeck, G. E. & Ornstein, L. S. 1930 On the theory of Brownian motion. Physical Review 36, 823841.CrossRefGoogle Scholar
Vasicek, O. A. 1977 An equilibrium characterization of the term structure. Journal of Financial Economics 5, 177-188.CrossRefGoogle Scholar
Weber, M. 1951 The fundamental solution of a degenerate partial differential equation of parabolic type. Transactions of the American Mathematical Society 71, 2437.CrossRefGoogle Scholar
Wong, M. W. 2014 Introduction to Pseudo-Differential Operators. Singapore: World Scientific Publishing Company.CrossRefGoogle Scholar
Zhang, Y., Chen, X. & Park, D. 2018 Formal specification of constant product (x × y = k) market maker model and implementation. White paper. https://pdfcoffee.com/uniswap-formulas-pdf-free.htmlGoogle Scholar
Figure 0

Figure 1 Newton’s letter to Oldenburg, 1676.

Reproduced by kind permission of the Syndics of Cambridge University Library.
Figure 1

Figure 2

Figure 2

Figure 2

Figure 3

Figure 3 Kelvin waves corresponding to two different orientations of the initial wave vector β0 and a0. (a), (b) β(0)=sinπ/4,0,cosπ/4,a0=0,sinπ/4,0; (c), (d) β(0)=sinπ/3,0,cosπ/3,a0=0,sinπ/3,0. Other parameters are as follows: T=100,ω=1,s=0.5. In the first case, at stays bounded, while at explodes in the second case. This explosion means that the underlying elliptic flow is unstable. Author’s graphics.

Figure 4

Figure 4 Kelvin waves in the viscous fluid with viscosity ν=0.07. Other parameters and initial conditions are the same as in Figure 3. Viscosity dampens the instability but, generally, does not suppress it entirely. Author’s graphics.

Figure 5

Figure 5

Figure 6

Figure 5

Reproduced by kind permission of the Editors of Annals of Mathematics.
Figure 7

Figure 6 A thousand trajectories of a typical Kolmogorov process. Parameters are as follows: T=5,dt=0.01,f=0.2,σ=0.8. (a) xt, (b) yt, (c) xˉT,yˉT, (d) contour lines of ϖ0,0,0,T,x˜,y˜. Author’s graphics.

Figure 8

Figure 7 A thousand trajectories of a typical free particle. Parameters are as follows: T=5,dt=0.01,κ=0.8,σ=1.0. (a) xt, (b) yt, (c) xˉT,yˉT, (d) contour lines of ϖ0,0,0,T,x˜,y˜. Author’s graphics.

Figure 9

Figure 8 A thousand trajectories of a harmonically bounded particle. Parameters are as follows: T=5,dt=0.01,κ=0.2,ω=0.5,σ=0.5. (a) xt, (b) yt, (c) xˉT,yˉT, (d) contour lines of ϖ0,0,0,T,x˜,y˜. Author’s graphics.

Figure 10

Figure 9 Contour lines of ϖ0,0,0,T,x˜,y˜ for an anomalous Kolmogorov process with T=1.5,a=2.5,b=1.5 . Author’s graphics.

Figure 11

Figure 10 T.p.d.fs for three Feller processes with different parameters and regularity conditions. (a) χ=0.1,κ=1.2,ε=0.2,y0=0.15,tˉmax=3; (b), (c) χ=0.1,κ=1.2,ε=0.6,y0=0.15,tˉmax=3. For the first and second processes, the probability of yˉ≥0 is equal to one. For the third process, this probability, shown as a function of time in (d), is less than one. Author’s graphics.

Figure 12

Figure 11 A thousand trajectories of a representative t.p.d.f. for the degenerate augmented Feller process. Parameters are T=3,dt=0.01,χ=0.1,κ=1.2,ε=0.2,x=0 ; y0=0.15. (a) xt, (b) yt, (c) xˉT,yˉT, (d) contour lines of ϖ0,0.15,0,T,x˜,y˜. Author’s graphics.

Figure 13

Figure 12 A thousand trajectories of a representative nondegenerate augmented Feller process. Parameters are T=3,dt=0.01,χ=0.2,κ=2.0,ε=0.2,ρ=−0.5,x=0 , y0=0.15. (a) xt, (b) yt, (c) xˉT,yˉT, (d) contour lines of ϖ0,0.15,0,T,x˜,y˜. Author’s graphics.

Figure 14

Figure 13 A representative implied volatility surface generated by the Heston model. Parameters are the same as in Figure 12. Author’s graphics.

Figure 15

Figure 14 The constant sum, constant product, and mixed-rule curves, along with the relative prices of TN2 in terms of TN1 and the associated impermanent losses; α=10. Author’s graphics.

Figure 16

Figure 15 The exact impermanent loss and its approximations via log and enthropy contracts. The corresponding parameters are T=3,dt=0.01,χ=0.2,κ=2.0,ε=0.2,ρ=−0.5,v=0.15. Author’s graphics.

Save element to Kindle

To save this element to your Kindle, first ensure [email protected] is added to your Approved Personal Document E-mail List under your Personal Document Settings on the Manage Your Content and Devices page of your Amazon account. Then enter the ‘name’ part of your Kindle email address below. Find out more about saving to your Kindle.

Note you can select to save to either the @free.kindle.com or @kindle.com variations. ‘@free.kindle.com’ emails are free but can only be saved to your device when it is connected to wi-fi. ‘@kindle.com’ emails can be delivered even when you are not connected to wi-fi, but note that service fees apply.

Find out more about the Kindle Personal Document Service.

Hydrodynamics of Markets
Available formats
×

Save element to Dropbox

To save content items to your account, please confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your account. Find out more about saving content to Dropbox.

Hydrodynamics of Markets
Available formats
×

Save element to Google Drive

To save content items to your account, please confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your account. Find out more about saving content to Google Drive.

Hydrodynamics of Markets
Available formats
×