Draft:Pearson diffusion

Pearson diffusions, also called Pearson diffusion processes, are a family of one-dimensional diffusion processes whose drift coefficient is a polynomial of degree at most one and whose squared diffusion coefficient is a polynomial of degree at most two. When an ergodic stationary distribution exists, its density belongs to the Pearson system of probability distributions.

The family contains several classical stochastic processes, including the Ornstein–Uhlenbeck process, the Cox–Ingersoll–Ross (CIR) process, Jacobi diffusion, Fisher–Snedecor diffusion, reciprocal gamma diffusion and Student diffusion. These processes provide dynamic counterparts of Gaussian, gamma, beta and several heavy-tailed probability distributions.

A central feature of Pearson diffusions is their polynomial structure. Their infinitesimal generators map polynomials to polynomials of no higher degree. As a consequence, conditional moments satisfy finite triangular systems of differential equations, allowing explicit formulas for quantities such as conditional expectation and conditional variance.

History

The statistical origin of Pearson diffusions goes back to Karl Pearson's work on frequency distributions in the late nineteenth century. In 1895, Pearson studied asymmetric frequency curves and introduced a differential relation for a broad family of probability densities extending beyond the normal distribution. This line of work developed into what became known as the Pearson system of distributions.

A separate development occurred in the theory of stochastic processes. In 1930, George Uhlenbeck and Leonard Ornstein studied the motion of a Brownian particle subject to friction, leading to the process now known as the Ornstein–Uhlenbeck process. Although it was not introduced as a Pearson diffusion, its linear mean-reverting drift, constant diffusion coefficient and Gaussian stationary distribution later placed it naturally in the Pearson family.

In 1931, Andrey Kolmogorov developed analytical methods for continuous-time Markov processes and differential equations governing their transition probabilities.

The direct connection between Pearson probability distributions and stationary diffusion processes became explicit in work by Eugene Wong in 1964. Wong considered the forward, or Fokker–Planck equation, together with probability densities satisfying the Pearson differential equation, thereby constructing stationary Markov processes having Pearson stationary distributions.

Other members of the family were introduced independently in applications. In 1985, John C. Cox, Jonathan E. Ingersoll and Stephen Ross introduced the square-root diffusion now known as the CIR process in their theory of the term structure of interest rates.

A systematic modern statistical treatment was given by Julie Lyng Forman and Michael Sørensen in 2008. They classified the ergodic Pearson diffusions, related their invariant laws to the complete Pearson system and derived polynomial conditional moments and polynomial eigenfunctions.

Later work extended the family beyond time-homogeneous diffusions. Leonenko, Meerschaert and Sikorskii introduced fractional Pearson diffusions through inverse-stable time changes, while Sutthimat and Mekchay studied time-inhomogeneous Pearson diffusions and obtained closed-form conditional moments without requiring explicit transition probability densities or eigenfunction expansions.

Mathematical definition

A time-homogeneous Pearson diffusion satisfies

$$dX_{t}=\theta (\mu -X_{t})\,dt+{\sqrt {2\theta \left(aX_{t}^{2}+bX_{t}+c\right)}}\,dW_{t},$$

where $\theta >0$, $W_{t}$ is a standard Wiener process, and the coefficients are chosen so that

$$aX_{t}^{2}+bX_{t}+c\geq 0$$

on the state space of the process.

The drift coefficient is linear,

$$\mu _{X}(x)=\theta (\mu -x),$$

while the squared diffusion coefficient is quadratic,

$$\sigma _{X}^{2}(x)=2\theta (ax^{2}+bx+c).$$

The infinitesimal generator is

$${\mathcal {G}}f(x)=\theta (\mu -x)f'(x)+\theta (ax^{2}+bx+c)f''(x).$$

Stationary distribution and Pearson equation

Let

$$D(x)=\theta (ax^{2}+bx+c).$$

If the process admits a stationary density $m(x)$ with zero stationary probability current, then

$$\theta (\mu -x)m(x)-{\frac {d}{dx}}\left[D(x)m(x)\right]=0.$$

Hence

$${\frac {m'(x)}{m(x)}}={\frac {\theta (\mu -x)-D'(x)}{D(x)}}.$$

Since $D(x)$ is quadratic, this expression has the form

$${\frac {m'(x)}{m(x)}}={\frac {\text{linear polynomial}}{\text{quadratic polynomial}}},$$

which is the characteristic differential structure of the Pearson system.

The stationary density may equivalently be written formally as

$$m(x)=C{\frac {1}{D(x)}}\exp \left[\int ^{x}{\frac {\theta (\mu -y)}{D(y)}}\,dy\right],$$

where $C$ is determined by normalization.

Polynomial structure and conditional moments

For a monomial $x^{n}$,

$${\mathcal {G}}x^{n}=-n\theta \left[1-(n-1)a\right]x^{n}+n\theta \left[\mu +(n-1)b\right]x^{n-1}+n(n-1)\theta c\,x^{n-2}.$$

Thus, the generator maps a polynomial of degree $n$ into another polynomial of degree at most $n$.

Define the conditional moment

$$M_{n}(t;s,x)=\mathbb {E} \left[X_{t}^{n}\mid X_{s}=x\right],\qquad t\geq s.$$

The moments satisfy the triangular system

$${\frac {dM_{n}}{dt}}=-n\theta \left[1-(n-1)a\right]M_{n}+n\theta \left[\mu +(n-1)b\right]M_{n-1}+n(n-1)\theta c\,M_{n-2}.$$

This polynomial closure is the basis of explicit conditional-moment calculations for Pearson diffusions.

Conditional expectation and variance

Sutthimat and Mekchay derived closed-form formulas for the first and second conditional moments of Pearson diffusion processes and used them to obtain conditional variances and related statistical quantities.

For the time-homogeneous Pearson diffusion, let

$$\tau =t-s,\qquad r=e^{-\theta \tau },\qquad d=x-\mu .$$

The conditional expectation has the particularly simple form

$$\mathbb {E} \left[X_{t}\mid X_{s}=x\right]=\mu +(x-\mu )e^{-\theta (t-s)}=\mu +dr.$$

Therefore every time-homogeneous Pearson diffusion has the same first-order mean-reversion structure. Differences between the individual Pearson classes first appear clearly in the second moment and conditional variance.

For

$$dX_{t}=\theta (\mu -X_{t})\,dt+{\sqrt {2\theta (aX_{t}^{2}+bX_{t}+c)}}\,dW_{t},$$

define

$$V_{\infty }={\frac {a\mu ^{2}+b\mu +c}{1-a}}.$$

When the second moment exists and $a\neq 1/2$, the conditional variance can be written as

$$\operatorname {Var} \left(X_{t}\mid X_{s}=x\right)=V_{\infty }\left(1-r^{\,2(1-a)}\right)$$

$$\quad +{\frac {2(2a\mu +b)}{1-2a}}(x-\mu )\left(r-r^{\,2(1-a)}\right)+(x-\mu )^{2}\left(r^{\,2(1-a)}-r^{2}\right).$$

At the exceptional value $a=1/2$, the corresponding expression is obtained by taking the continuous limit of this formula.

The formulas below are specializations of this general Pearson moment structure.

Ornstein–Uhlenbeck process

The Ornstein–Uhlenbeck process satisfies

$$dX_{t}=\theta (\mu -X_{t})\,dt+{\sqrt {2\theta \sigma ^{2}}}\,dW_{t}.$$

The conditional expectation is

$$\mathbb {E} \left[X_{t}\mid X_{s}=x\right]=\mu +(x-\mu )e^{-\theta (t-s)}.$$

Its conditional variance is

$$\operatorname {Var} \left(X_{t}\mid X_{s}=x\right)=\sigma ^{2}\left[1-e^{-2\theta (t-s)}\right].$$

As $t-s\to \infty$,

$$\mathbb {E} \left[X_{t}\mid X_{s}=x\right]\longrightarrow \mu$$

and

$$\operatorname {Var} \left(X_{t}\mid X_{s}=x\right)\longrightarrow \sigma ^{2},$$

which are the mean and variance of the stationary Gaussian distribution.

Cox–Ingersoll–Ross process

A convenient Pearson parameterization of the CIR process is

$$dX_{t}=\theta \left({\frac {\beta }{\alpha }}-X_{t}\right)dt+{\sqrt {{\frac {2\theta }{\alpha }}X_{t}}}\,dW_{t},$$

where $\alpha >0$ and $\beta >0$.

Its stationary law is gamma with mean

$$\mu _{\mathrm {CIR} }={\frac {\beta }{\alpha }}.$$

The conditional expectation is

$$\mathbb {E} \left[X_{t}\mid X_{s}=x\right]={\frac {\beta }{\alpha }}+\left(x-{\frac {\beta }{\alpha }}\right)e^{-\theta (t-s)}.$$

Writing

$$r=e^{-\theta (t-s)},$$

the conditional variance is

$$\operatorname {Var} \left(X_{t}\mid X_{s}=x\right)={\frac {\beta }{\alpha ^{2}}}\left(1-r^{2}\right)+{\frac {2}{\alpha }}\left(x-{\frac {\beta }{\alpha }}\right)\left(r-r^{2}\right).$$

Equivalently,

$$\operatorname {Var} \left(X_{t}\mid X_{s}=x\right)={\frac {2x}{\alpha }}r(1-r)+{\frac {\beta }{\alpha ^{2}}}(1-r)^{2}.$$

The long-run variance is therefore

$${\frac {\beta }{\alpha ^{2}}}.$$

Jacobi diffusion

A Jacobi diffusion on $(-1,1)$ may be written as

$$dX_{t}=\theta (\mu _{J}-X_{t})\,dt+{\sqrt {{\frac {2\theta }{K}}(1-X_{t}^{2})}}\,dW_{t},$$

where

$$K=\alpha +\beta +2$$

and

$$\mu _{J}={\frac {\beta -\alpha }{K}},\qquad \alpha >-1,\qquad \beta >-1.$$

Its conditional expectation is

$$\mathbb {E} \left[X_{t}\mid X_{s}=x\right]=\mu _{J}+(x-\mu _{J})e^{-\theta (t-s)}.$$

Its stationary variance is

$$V_{J}={\frac {4(\alpha +1)(\beta +1)}{K^{2}(K+1)}}.$$

Let

$$p_{J}={\frac {2(K+1)}{K}}$$

and

$$r=e^{-\theta (t-s)}.$$

Then the conditional variance takes the form

$$\operatorname {Var} \left(X_{t}\mid X_{s}=x\right)=V_{J}(1-r^{p_{J}})$$

$$\quad -{\frac {4\mu _{J}(x-\mu _{J})}{K+2}}\left(r-r^{p_{J}}\right)+(x-\mu _{J})^{2}\left(r^{p_{J}}-r^{2}\right).$$

Thus the bounded state space changes the rate at which the second conditional moment approaches its stationary value.

Fisher–Snedecor diffusion

A Fisher–Snedecor diffusion may be written as

$$dX_{t}=\theta \left({\frac {\beta }{\beta -2}}-X_{t}\right)dt+{\sqrt {{\frac {4\theta }{\alpha (\beta -2)}}X_{t}(\alpha X_{t}+\beta )}}\,dW_{t}.$$

Its stationary mean is

$$\mu _{F}={\frac {\beta }{\beta -2}},\qquad \beta >2.$$

The conditional expectation is

$$\mathbb {E} \left[X_{t}\mid X_{s}=x\right]=\mu _{F}+(x-\mu _{F})e^{-\theta (t-s)}.$$

When $\beta >4$, the stationary variance exists and is

$$V_{F}={\frac {2\beta ^{2}(\alpha +\beta -2)}{\alpha (\beta -4)(\beta -2)^{2}}}.$$

Define

$$p_{F}={\frac {2(\beta -4)}{\beta -2}},\qquad r=e^{-\theta (t-s)}.$$

For $\beta \neq 6$, the conditional variance can be written as

$$\operatorname {Var} \left(X_{t}\mid X_{s}=x\right)=V_{F}\left(1-r^{p_{F}}\right)$$

$$\quad +{\frac {4\beta (2\alpha +\beta -2)}{\alpha (\beta -6)(\beta -2)}}(x-\mu _{F})\left(r-r^{p_{F}}\right)+(x-\mu _{F})^{2}\left(r^{p_{F}}-r^{2}\right).$$

The case $\beta =6$ is obtained from the continuous limit of the formula.

Unlike the OU, CIR and Jacobi models, the Fisher–Snedecor diffusion is heavy-tailed, and only moments below parameter-dependent orders are finite.

Reciprocal gamma diffusion

A reciprocal gamma diffusion may be written as

$$dX_{t}=\theta \left({\frac {\alpha }{\beta -1}}-X_{t}\right)dt+{\sqrt {{\frac {2\theta }{\beta -1}}X_{t}^{2}}}\,dW_{t}.$$

Its stationary mean is

$$\mu _{R}={\frac {\alpha }{\beta -1}}.$$

The conditional expectation is

$$\mathbb {E} \left[X_{t}\mid X_{s}=x\right]=\mu _{R}+(x-\mu _{R})e^{-\theta (t-s)}.$$

For $\beta >2$, the stationary variance is

$$V_{R}={\frac {\alpha ^{2}}{(\beta -1)^{2}(\beta -2)}}.$$

Define

$$p_{R}={\frac {2(\beta -2)}{\beta -1}}.$$

For $\beta \neq 3$,

$$\operatorname {Var} \left(X_{t}\mid X_{s}=x\right)=V_{R}\left(1-r^{p_{R}}\right)$$

$$\quad +{\frac {4\mu _{R}(x-\mu _{R})}{\beta -3}}\left(r-r^{p_{R}}\right)+(x-\mu _{R})^{2}\left(r^{p_{R}}-r^{2}\right),$$

where

$$r=e^{-\theta (t-s)}.$$

The expression for $\beta =3$ follows by continuity.

Student diffusion

For a symmetric Student diffusion, a convenient parameterization is

$$dX_{t}=\theta (\mu -X_{t})\,dt+{\sqrt {{\frac {2\theta }{\nu -1}}\left[(X_{t}-\mu )^{2}+\delta ^{2}\right]}}\,dW_{t},$$

where

$$\nu >1,\qquad \delta >0.$$

Its conditional expectation is

$$\mathbb {E} \left[X_{t}\mid X_{s}=x\right]=\mu +(x-\mu )e^{-\theta (t-s)}.$$

When $\nu >2$, its stationary variance is

$$V_{S}={\frac {\delta ^{2}}{\nu -2}}.$$

Define

$$p_{S}={\frac {2(\nu -2)}{\nu -1}}.$$

Then

$$\operatorname {Var} \left(X_{t}\mid X_{s}=x\right)=V_{S}\left(1-r^{p_{S}}\right)+(x-\mu )^{2}\left(r^{p_{S}}-r^{2}\right),$$

where

$$r=e^{-\theta (t-s)}.$$

This form illustrates a basic feature of heavy-tailed Pearson diffusions: the conditional expectation has the same exponential mean-reversion form as in the light-tailed models, but the second moment relaxes at a different parameter-dependent rate.

More general Student-type Pearson diffusions may be asymmetric. In that case the general Pearson conditional-variance formula applies with the corresponding quadratic coefficients $a$, $b$ and $c$.

Interpretation of the conditional moments

The six Pearson classes share the same first-order mean-reversion mechanism,

$$\mathbb {E} \left[X_{t}\mid X_{s}=x\right]-\mu =(x-\mu )e^{-\theta (t-s)}.$$

Therefore $\theta$ controls the speed with which the conditional mean loses information about the initial state.

The conditional variances behave differently because they depend on the quadratic diffusion coefficient. For the OU process the variance is independent of the starting state. For CIR and other state-dependent diffusions, the conditional variance depends explicitly on $x$. For heavy-tailed Pearson diffusions, the rate of relaxation of the second moment is determined by the quadratic coefficient and may differ substantially from the rate $2\theta$ of the OU process.

This distinction can be seen from the exponent

$$2\theta (1-a),$$

which governs part of the second-moment dynamics. The first conditional moment therefore depends only on the linear drift, while higher conditional moments also encode the nonlinear structure of the diffusion coefficient.

Sutthimat and Mekchay used this polynomial moment structure to derive conditional moments directly from the associated differential equations, without requiring explicit transition densities or eigenfunctions. Their formulas also yield conditional variance, mixed moments, covariance and correlation.

A later study by Wu cited this work as providing exact conditional-moment formulas for Pearson diffusion processes.

Polynomial eigenfunctions

The polynomial-preserving property also appears in the spectral structure of Pearson diffusions. Polynomial eigenfunctions satisfy

$${\mathcal {G}}Q_{n}(x)=-\lambda _{n}Q_{n}(x).$$

For the general Pearson generator,

$$\lambda _{n}=n\theta \left[1-(n-1)a\right].$$

When the eigenfunctions belong to the discrete spectrum, $Q_{n}$ is a polynomial of degree $n$. Under appropriate normalization,

$$\int _{E}Q_{n}(x)Q_{k}(x)m(x)\,dx=\delta _{nk}.$$

The principal light-tailed Pearson diffusions are connected with classical orthogonal polynomials:

• the Ornstein–Uhlenbeck process is associated with Hermite polynomials;
• the CIR process is associated with Laguerre polynomials;
• the Jacobi diffusion is associated with Jacobi polynomials.

The heavy-tailed members have finite families of square-integrable polynomial eigenfunctions together with continuous spectral components.

For Pearson diffusions with a purely discrete spectrum, the transition density has the spectral expansion

$$p(t,y\mid x)=m(y)\sum _{n=0}^{\infty }e^{-\lambda _{n}t}Q_{n}(x)Q_{n}(y).$$

Spectral classification

Pearson diffusions may be divided into three spectral categories.

The first category contains the Ornstein–Uhlenbeck, CIR and Jacobi diffusions. Their generators have purely discrete spectra.

The second category contains Fisher–Snedecor and reciprocal gamma diffusions. Their spectra contain a finite discrete component together with an absolutely continuous component.

The third category contains Student diffusions. Their spectra also contain discrete and continuous parts, with the continuous component having multiplicity two.

For the heavy-tailed cases the transition density may therefore be written schematically as

$$p(t,y\mid x)=p_{\mathrm {d} }(t,y\mid x)+p_{\mathrm {c} }(t,y\mid x),$$

where

$$p_{\mathrm {d} }(t,y\mid x)=m(y)\sum _{n=0}^{N}e^{-\lambda _{n}t}Q_{n}(x)Q_{n}(y)$$

is the discrete part and $p_{\mathrm {c} }$ is determined by the continuous spectrum.

Fractional Pearson diffusions

Fractional Pearson diffusions introduce memory into the classical Markov model by replacing the ordinary time derivative with a fractional derivative.

For a Caputo derivative of order

$$0<\gamma <1,$$

the backward equation takes the form

$${}^{C}D_{t}^{\gamma }u(t,x)={\mathcal {G}}u(t,x).$$

For processes with a discrete spectrum,

$$u(t,x)=\sum _{n=0}^{\infty }f_{n}E_{\gamma }\left(-\lambda _{n}t^{\gamma }\right)Q_{n}(x),$$

where $E_{\gamma }$ is the Mittag-Leffler function.

The fractional transition density correspondingly becomes

$$p_{\gamma }(t,y\mid x)=m(y)\sum _{n=0}^{\infty }E_{\gamma }\left(-\lambda _{n}t^{\gamma }\right)Q_{n}(x)Q_{n}(y).$$

When $\gamma =1$,

$$E_{1}(-\lambda t)=e^{-\lambda t},$$

recovering the classical Pearson diffusion.

A probabilistic construction is obtained by an inverse-stable time change,

$$X_{\gamma }(t)=X(L_{\gamma }(t)),$$

where $L_{\gamma }$ is the inverse of a stable subordinator independent of $X$.

Time-non-local Pearson diffusions

The fractional construction can be generalized using a Bernstein function $\Phi$. Let $\sigma _{\Phi }$ be the associated subordinator and

$$L_{\Phi }(t)=\inf \left\{s>0:\sigma _{\Phi }(s)>t\right\}$$

its inverse.

The corresponding time-changed Pearson diffusion is

$$X_{\Phi }(t)=X(L_{\Phi }(t)).$$

If $f_{\Phi }(s;t)$ denotes the density of $L_{\Phi }(t)$, then

$$p_{\Phi }(t,y\mid x)=\int _{0}^{\infty }p(s,y\mid x)f_{\Phi }(s;t)\,ds.$$

For processes in the first spectral category,

$$p_{\Phi }(t,y\mid x)=m(y)\sum _{n=0}^{\infty }{\mathfrak {e}}_{\Phi }(t;-\lambda _{n})Q_{n}(x)Q_{n}(y),$$

where

$${\mathfrak {e}}_{\Phi }(t;-\lambda )=\mathbb {E} \left[e^{-\lambda L_{\Phi }(t)}\right].$$

Time-inhomogeneous Pearson diffusions

Sutthimat and Mekchay considered Pearson diffusions in which the model parameters depend explicitly on time. Their model is

$$dX_{t}=\theta (t)\left[\mu (t)-X_{t}\right]dt+{\sqrt {2\theta (t)\left[a(t)X_{t}^{2}+b(t)X_{t}+c(t)\right]}}\,dW_{t}.$$

The time-dependent generator is

$${\mathcal {G}}_{t}f(x)=\theta (t)\left[\mu (t)-x\right]f'(x)+\theta (t)\left[a(t)x^{2}+b(t)x+c(t)\right]f''(x).$$

The polynomial property is retained:

$${\mathcal {G}}_{t}x^{n}=n\theta (t)\left[(n-1)a(t)-1\right]x^{n}$$

$$\quad +n\theta (t)\left[\mu (t)+(n-1)b(t)\right]x^{n-1}+n(n-1)\theta (t)c(t)x^{n-2}.$$

Consequently,

$$M_{n}(t;s,x)=\mathbb {E} \left[X_{t}^{n}\mid X_{s}=x\right]$$

satisfies

$${\frac {d}{dt}}M_{n}(t;s,x)=n\theta (t)\left[(n-1)a(t)-1\right]M_{n}(t;s,x)$$

$$\quad +n\theta (t)\left[\mu (t)+(n-1)b(t)\right]M_{n-1}(t;s,x)+n(n-1)\theta (t)c(t)M_{n-2}(t;s,x).$$

For the first moment,

$${\frac {d}{dt}}M_{1}(t;s,x)=\theta (t)\left[\mu (t)-M_{1}(t;s,x)\right].$$

Define

$$A(u,t)=\exp \left[-\int _{u}^{t}\theta (r)\,dr\right].$$

Then

$$\mathbb {E} \left[X_{t}\mid X_{s}=x\right]=A(s,t)x+\int _{s}^{t}A(u,t)\theta (u)\mu (u)\,du.$$

This formula reduces to

$$\mu +(x-\mu )e^{-\theta (t-s)}$$

when $\theta$ and $\mu$ are constant.

For the second conditional moment,

$${\frac {d}{dt}}M_{2}=2\theta (t)\left[a(t)-1\right]M_{2}+2\theta (t)\left[\mu (t)+b(t)\right]M_{1}+2\theta (t)c(t).$$

The conditional variance is then

$$\operatorname {Var} \left(X_{t}\mid X_{s}=x\right)=M_{2}(t;s,x)-M_{1}(t;s,x)^{2}.$$

This provides a unified route to conditional expectation, variance, covariance, correlation and higher conditional moments without first obtaining the transition probability density.

Statistical applications

Conditional moments can be used in parameter estimation because, for the correct model,

$$X_{t_{i}}^{k}-\mathbb {E} \left[X_{t_{i}}^{k}\mid X_{t_{i-1}}\right]$$

has conditional mean zero.

This leads to martingale estimating equations such as

$$\sum _{i=1}^{N}h(X_{t_{i-1}})\left\{X_{t_{i}}^{k}-\mathbb {E} \left[X_{t_{i}}^{k}\mid X_{t_{i-1}}\right]\right\}=0.$$

Explicit conditional moments therefore make Pearson diffusions useful for moment-based statistical inference, generalized method-of-moments procedures and related estimation methods.

See also

• Diffusion process
• Itô diffusion
• Stochastic differential equation
• Itô's lemma
• Markov process
• Kolmogorov equations
• Fokker–Planck equation
• Infinitesimal generator
• Stationary distribution
• Pearson distribution
• Orthogonal polynomials
• Ornstein–Uhlenbeck process
• Cox–Ingersoll–Ross model
• Gamma distribution
• Beta distribution
• Student's t-distribution
• Conditional expectation
• Variance
• Fractional calculus
• Mittag-Leffler function