API reference

ContinuousTimeMarkovProcess{N}

Abstract supertype for continuous-time Markov processes with N tensor-product state-space axes.

Concrete subtypes must define state_space(X) and generator(X). state_space(X) must return a tuple of state-space axes, including for one-dimensional processes. size(X) is derived from those axes, and ndims(X) returns the type parameter N.

ContinuousTimeMarkovChain(states, Q)
ContinuousTimeMarkovChain(Q)

Returns a finite-state continuous-time Markov chain on the discrete state space states, with generator (transition-rate) matrix Q.

The matrix Q is a generator: rows must sum to zero and off-diagonal entries must be nonnegative transition rates. It is not a discrete-time transition probability matrix.

DiffusionProcess(x, μx, σx)

Returns the diffusion process

dx_t = μ(x_t) dt + σ(x_t) dZ_t

discretized on the strictly increasing grid x (possibly non-uniform), where μx and σx are the drift and volatility evaluated on the grid. The drift is discretized by upwinding and the boundaries are reflecting, so generator always returns a valid generator (transition-rate) matrix.

OrnsteinUhlenbeck(; xbar = 0.0, κ = 0.1, σ = 1.0, p = 1e-10, length = 100, xmin, xmax, pow = 1)

Returns the Ornstein–Uhlenbeck process

dx_t = -κ * (x_t - xbar) * dt + σ * dZ_t

discretized as a DiffusionProcess on an automatically chosen grid.

By default the grid has length points spanning the p and 1 - p quantiles of the stationary distribution N(xbar, σ^2 / 2κ); pass xmin and xmax to override the limits. pow controls the grid spacing (when xmin > 0, points are uniform in x^(1/pow), so pow > 1 concentrates points near xmin). The default p = 1e-10 is deliberately extreme: reflecting boundaries distort solutions near the edges of the grid, and a wide grid pushes that distortion where the process never goes — which matters especially for tail_index.

CoxIngersollRoss(; xbar = 0.1, κ = 0.1, σ = 1.0, p = 1e-10, length = 100, xmin, xmax, pow = 2)

Returns the Cox–Ingersoll–Ross process

dx_t = -κ * (x_t - xbar) * dt + σ * sqrt(x_t) * dZ_t

discretized as a DiffusionProcess on an automatically chosen grid. Requires the Feller condition 2κ * xbar / σ^2 > 1, so that 0 is not attainable.

By default the grid has length points spanning the p and 1 - p quantiles of the stationary Gamma distribution; pass xmin and xmax to override the limits. pow controls the grid spacing (points are uniform in x^(1/pow), so the default pow = 2 concentrates points near zero, where the volatility is smallest).

MultivariateDiffusionProcess(grid; drift, variance, covariance = nothing)

Returns a diffusion process on a tensor-product grid.

grid is a NamedTuple mapping state names to strictly increasing grid vectors. drift has the same names as grid; each entry is either a scalar or an array with the same size as the tensor grid. variance has the same names and gives the diagonal entries of the instantaneous covariance matrix. Cross terms are optional and are passed as pair names such as (; xy = axy).

ProductProcess(X, Y, ...)

Returns the independent product of Markov processes in argument order.

SwitchingProcess(Z, Xs)

Returns a Markov process whose dynamics are modulated by the autonomous Markov process Z — any process in the package, e.g. a ContinuousTimeMarkovChain (regime switching) or a DiffusionProcess (continuous modulation). In the i-th state of Z (in column-major order for multivariate Z), the process follows Xs[i]; all Xs must share the same state space.

state_space(X::ContinuousTimeMarkovProcess)

Return the state-space axes of X as a tuple.

The return value is always a tuple, including for one-dimensional processes: state_space(X)[1] is the first axis, and only(state_space(X)) is the grid of a one-dimensional process.

size(X::ContinuousTimeMarkovProcess)

Return the tensor shape of arrays defined on the state space of X.

This is the length of each state-space axis. Operators such as stationary_distribution(X) and feynman_kac(X, ...) use this shape for process-level inputs and outputs.

length(X::ContinuousTimeMarkovProcess)

Return the total number of grid points or finite states in X, equal to prod(size(X)). This is the length of flattened state vectors and the row/column dimension of generator(X).

ndims(X::ContinuousTimeMarkovProcess)

Return the number of tensor-product state-space axes of X.

generator(X::ContinuousTimeMarkovProcess)

Return the generator (transition-rate) matrix 𝔸 of the discretized process — the matrix representation of the infinitesimal generator

𝔸: f ↦ lim_{t→0} (E[f(x_t) | x_0 = x] - f(x)) / t

acting on the state space flattened in column-major order. Rows sum to zero and off-diagonal entries are non-negative transition rates.

check_generator(𝔸::AbstractMatrix; atol = nothing)

Check that 𝔸 is a valid generator (transition-rate) matrix: square, with nonnegative off-diagonal entries and rows summing to zero — that is, a Metzler matrix whose rows sum to zero.

A violation beyond atol emits a warning reporting its size rather than throwing — a matrix that is close to a generator often still yields accurate results. Non-square or empty matrices throw an ArgumentError. Returns 𝔸, so the call can be used inline.

atol defaults to sqrt(eps) * max(1, maximum(abs, diag(𝔸))) (for a valid generator, the largest entry in absolute value is on the diagonal). The off-diagonal check is minimum(𝔸 - Diagonal(diag(𝔸))), which preserves the structure of the matrix — sparse matrices are checked in one pass over their stored entries, tridiagonal ones over their bands.

Note that a tilted generator (tilted_generator(m, ξ)) is intentionally not a generator in this sense: it has nonnegative off-diagonal entries but nonzero row sums — that is why E[e^{ξ mₜ}] grows at rate Λ(ξ) instead of remaining a probability.

stationary_distribution(𝔸::AbstractMatrix; δ = 0.0, rebirth = Ones(size(𝔸, 1)))

Compute the stationary distribution corresponding to the generator matrix 𝔸 — the solution of the Kolmogorov forward equation 𝔸'g = 0, normalized to sum to one.

The result is a vector of probability masses, one per state — sum(g) == 1 with no grid weights — so expectations are unweighted dot products sum(g .* f). To get a density for a discretized diffusion, divide the masses by the cell widths of the grid.

With a death rate δ > 0 and a rebirth distribution ψ given by rebirth, returns instead the resolvent (δI - 𝔸')⁻¹ δψ — the stationary distribution when units die at rate δ and are reborn with distribution ψ.

stationary_distribution(X::ContinuousTimeMarkovProcess; δ = 0.0, rebirth = Ones(length(X)), kwargs...)

Compute the stationary distribution of the Markov process X, returned as an array of probability masses shaped like size(X): the array sums to one with no grid weights, and expectations are unweighted dot products sum(g .* f). To get a density for a diffusion, divide by the cell widths of the grid — essential when the grid is non-uniform.

The keywords δ and rebirth have the same meaning as in the matrix form; rebirth is an array shaped like size(X). Remaining keyword arguments are forwarded to generator(X) — e.g. check = :warn for a MultivariateDiffusionProcess.

feynman_kac(𝔸::AbstractMatrix, ts; f, ψ, v, direction = :backward)

Solve the Feynman–Kac PDE associated with a generator matrix 𝔸 on the increasing time grid ts, using implicit Euler time steps.

With direction = :backward, returns the solution of

u(x, ts[end]) = ψ(x)
0 = uₜ + 𝔸u - v(x, t)u + f(x, t)

or, equivalently, in integral form,

u(x, t) = E[∫ₜᵀ e^{-∫ₜˢ v(x_u) du} f(x_s) ds + e^{-∫ₜᵀ v(x_u) du} ψ(x_T) | xₜ = x].

With direction = :forward, returns the solution of

u(x, ts[1]) = ψ(x)
uₜ = 𝔸u - v(x, t)u + f(x, t)

or, equivalently, in integral form,

u(x, t) = E[∫₀ᵗ e^{-∫₀ˢ v(x_u) du} f(x_s) ds + e^{-∫₀ᵗ v(x_u) du} ψ(xₜ) | x₀ = x].

The inputs f, ψ, and v are vectors of length size(𝔸, 1). The inputs f and v may also be matrices with length(ts) columns (one per date) or with a single column (held constant over time).

feynman_kac(X::ContinuousTimeMarkovProcess, ts; f = nothing, ψ = nothing, v = nothing, direction = :backward, kwargs...)

Solve the Feynman–Kac PDE associated with a Markov process X on the time grid ts, using implicit Euler time steps.

For the process form, f, ψ, and v are arrays shaped like the state space (size(X)); f and v may also carry a trailing time dimension of length length(ts). The result has one slice per date in ts. Remaining keyword arguments are forwarded to generator(X) — e.g. check = :warn for a MultivariateDiffusionProcess.

jointoperator(operators, Q)

Combine one generator matrix per regime with a regime generator (transition-rate) matrix Q into the sparse block generator of the joint process. Most users can use SwitchingProcess instead.

principal_eigenvalue(𝔸; r0 = ones(size(𝔸, 1)), η0 = nothing, maxiter = 100, tol = 1e-12)

Compute the principal eigenvalue and eigenvector of a Metzler matrix 𝔸 (i.e. a matrix with non-negative off-diagonal entries), returned as a tuple (η, r).

By Perron-Frobenius, the eigenvalue η with largest real part is real, and the corresponding eigenvector r is strictly positive. Metzlerity is what makes this hold — generators and tilted generators both satisfy it — so the matrix is checked on entry, with a warning if an off-diagonal entry is negative (results may still be accurate when the violation is small).

Two cases:

1. If rows or columns sum to zero (𝔸 is a generator, or its transpose), then η = 0. The eigenvector is found by solving 𝔸r = 0 with r[1] = 1; the row-sum test is scaled by the size of the diagonal, so it is robust to the rounding of generator assembly. If that solve is singular, the generator is reducible (it has states that do not communicate), the stationary eigenvector is not unique, and an ArgumentError says so.
2. Otherwise, η is found by inverse iteration with Rayleigh quotient updates. At each step, we solve (𝔸 - σI)w = r, normalize w, and update the shift σ via the Rayleigh quotient σ = r'𝔸r. This converges cubically to a nearby eigenvalue. The initial shift is the max row sum — a Gershgorin upper bound on the real part of the spectrum. When 𝔸 has a real spectrum (e.g. the tridiagonal generators of 1-D diffusions, which are similar to symmetric matrices), the eigenvalue nearest that shift is the principal one, so the iteration lands on η. For a general Metzler matrix with complex eigenvalues this is not guaranteed; pass η0 to start from a known bound if needed. If a shift lands exactly on an eigenvalue (e.g. the Gershgorin bound when all row sums are equal, so the bound is attained), it is nudged by tol to keep the solve nonsingular; the Rayleigh quotient then recovers the exact eigenvalue. For tridiagonal 𝔸, each iteration costs O(n).

AdditiveFunctional

Abstract supertype for additive functionals m of a Markov process — cumulative quantities whose increments depend on the state, such as cumulative log growth. The defining verb is tilted_generator(m, ξ); the operators cgf, tail_index, and principal_eigenvalue are computed from it.

tilted_generator(m::AdditiveFunctional, ξ)

Return the tilted generator matrix 𝔸_ξ of the additive functional m, i.e. the matrix representation of the operator

𝔸_ξ: f ↦ lim_{t→0} (E[e^{ξ mₜ} f(xₜ) | x₀ = x] - f(x)) / t.

tilted_generator(m, 0) is the generator of the underlying process. Everything the package computes about additive functionals derives from this matrix: cgf and tail_index use its principal eigenvalue, and it can be passed directly to feynman_kac to compute finite-horizon moments E[e^{ξ mₜ} ψ(xₜ)]. To define a custom additive functional, subtype AdditiveFunctional and define this function.

A tilted generator is not a generator in the sense of check_generator: its off-diagonal entries must be nonnegative (a Metzler matrix — equivalently, a generator plus a diagonal), but its rows need not sum to zero. Metzlerity is what makes the principal eigenvalue real and its eigenvector positive; a custom implementation must preserve it (principal_eigenvalue checks it on entry and warns if it fails).

cgf(m::AdditiveFunctional, ξ; r0 = nothing, η0 = nothing)

Return the long-run scaled cumulant generating function of m evaluated at ξ,

Λ(ξ) = lim_{t→∞} log(E[e^{ξ mₜ}]) / t,

computed as the principal eigenvalue of the tilted generator tilted_generator(m, ξ) by inverse iteration (Hansen and Scheinkman 2009). r0 is an initial guess for the eigenvector (defaults to a vector of ones) and η0 an initial guess for the eigenvalue. Use cgf_eigenvector to also obtain the associated eigenvector.

cgf_eigenvector(m::AdditiveFunctional, ξ, side = :right; r0 = nothing, η0 = nothing)

Return the pair (Λ(ξ), vector): the long-run scaled cumulant generating function of m at ξ (as in cgf) together with the associated principal eigenvector of the tilted generator.

With side = :right, the right eigenvector — the Hansen–Scheinkman eigenfunction. With side = :left, the left eigenvector, normalized to sum to one — the stationary distribution of the twisted process, i.e. the distribution over states from which the long-run growth of E[e^{ξ mₜ}] is achieved.

tail_index(m::AdditiveFunctional; δ = 0, bracket = (1e-5, 1e3), xatol = 1e-4)

Compute the tail index of the stationary distribution of e^m when units die (are reset) at rate δ, i.e. the ζ such that cgf(m, ζ) = δ.

The root is searched for in bracket; if cgf(m, ξ) - δ has the same sign at both ends, an ArgumentError reports the two values so the bracket can be moved. Remaining keyword arguments are passed to Roots.fzero.

tail_index(μ::Number, σ::Number; δ = 0)

Closed form for constant coefficients: the tail index of the stationary distribution of a size w growing as dw/w = μ dt + σ dZ (note that μ is the arithmetic growth rate of w itself, equal to the drift of log w plus σ^2/2) with death rate δ.

AdditiveFunctionalDiffusion(X::DiffusionProcess, μm, σm; ρ = 0.0)

An additive functional m of the diffusion X, defined by

dmₜ = μm(xₜ) dt + σm(xₜ) dZᵐₜ,        corr(dZᵐₜ, dZₜ) = ρ,

where μm and σm are vectors evaluated on the grid of X.

Legacy type: AdditiveFunctional(X, μm, σm; ρ = ρ) constructs an equivalent functional, and is the preferred entry point; you rarely need to construct this type directly.

FirstDerivative(x, f; direction = :forward, bc = (0, 0))
FirstDerivative(grid, F, dim; direction = :forward, bc = (0, 0))

Lazily compute first-order derivatives (using finite-difference scheme) on a grid. The result is an AbstractArray whose entries are computed on demand; use collect to materialize it.

In one dimension, x is a grid vector and f is a vector with length(f) == length(x). In multiple dimensions, grid is a tuple of axis vectors such as (xs, ys), F is an array with size (length(xs), length(ys)), and dim is the dimension number to differentiate. NamedTuple grids with symbolic dimensions are also accepted as a convenience.

direction selects the one-sided difference: :forward (synonyms :up, :upward) or :backward (synonyms :down, :downward).

bc is the value of the first derivative at the lower and upper boundaries.

SecondDerivative(x, f; bc = (0, 0))
SecondDerivative(grid, F, dim; bc = (0, 0))
SecondDerivative(grid, F, dim1, dim2; direction = :up)

Lazily compute second-order derivatives (using finite-difference scheme) on a grid. The result is an AbstractArray whose entries are computed on demand; use collect to materialize it.

In one dimension, x is a grid vector and f is a vector with length(f) == length(x). In multiple dimensions, grid is a tuple of axis vectors such as (xs, ys), and F is an array with size (length(xs), length(ys)). NamedTuple grids with symbolic dimensions are also accepted as a convenience.

Use SecondDerivative(grid, F, 1, 1) for own second derivatives and SecondDerivative(grid, F, 1, 2; direction = :up) for directional cross derivatives. The :up direction (synonyms :forward, :upward) is the main-diagonal stencil; :down (synonyms :backward, :downward) is the anti-diagonal stencil.

For own second derivatives, bc is the value of the first derivative at the lower and upper boundaries. Cross derivatives do not accept boundary conditions (the stencil is clamped at the edges of the grid) and throw if a nonzero bc is passed.