Operators

Once a process is represented by its generator matrix, everything else is linear algebra. This page lists the operators that apply to any ContinuousTimeMarkovProcess{N} — univariate or multivariate — with pointers to the tutorials that derive each one by hand.

Throughout, take a simple process:

using InfinitesimalGenerators

X = OrnsteinUhlenbeck(; xbar = 1.0, κ = 0.1, σ = 0.05)
xs = only(state_space(X))

The generator itself

generator(X) returns the generator (transition-rate) matrix $\mathbb{A}$ of the discretized process — the operator $f \mapsto \lim_{t \downarrow 0} \left(E[f(x_t) \mid x_0 = x] - f(x)\right)/t$ acting on the flattened state space. Rows sum to zero and off-diagonals are non-negative, so the matrix can be used directly for anything the package does not provide: exponentials for transition probabilities, resolvents for present values, transposes for densities. Conversely, check_generator(𝔸) checks a matrix you built by hand (rows sum to zero, non-negative off-diagonals) before using it with the package, warning with the size of any violation.

𝔸 = generator(X)

100×100 LinearAlgebra.Tridiagonal{Float64, Vector{Float64}}:
 -11.005     11.005       ⋅       …     ⋅          ⋅          ⋅ 
   6.05498  -16.96      10.905          ⋅          ⋅          ⋅ 
    ⋅         6.05498  -16.86           ⋅          ⋅          ⋅ 
    ⋅          ⋅         6.05498        ⋅          ⋅          ⋅ 
    ⋅          ⋅          ⋅             ⋅          ⋅          ⋅ 
    ⋅          ⋅          ⋅       …     ⋅          ⋅          ⋅ 
    ⋅          ⋅          ⋅             ⋅          ⋅          ⋅ 
    ⋅          ⋅          ⋅             ⋅          ⋅          ⋅ 
    ⋅          ⋅          ⋅             ⋅          ⋅          ⋅ 
    ⋅          ⋅          ⋅             ⋅          ⋅          ⋅ 
   ⋮                              ⋱                        
    ⋅          ⋅          ⋅             ⋅          ⋅          ⋅ 
    ⋅          ⋅          ⋅             ⋅          ⋅          ⋅ 
    ⋅          ⋅          ⋅             ⋅          ⋅          ⋅ 
    ⋅          ⋅          ⋅             ⋅          ⋅          ⋅ 
    ⋅          ⋅          ⋅       …     ⋅          ⋅          ⋅ 
    ⋅          ⋅          ⋅            6.05498     ⋅          ⋅ 
    ⋅          ⋅          ⋅          -16.86       6.05498     ⋅ 
    ⋅          ⋅          ⋅           10.905    -16.96       6.05498
    ⋅          ⋅          ⋅             ⋅        11.005    -11.005

Stationary distributions

stationary_distribution(X) solves the Kolmogorov forward equation $\mathbb{A}' g = 0$ and returns the probability mass at each grid point, shaped like size(X); to get a density, divide by the cell widths. See the distribution dynamics tutorial for the computation by hand:

g = stationary_distribution(X)
sum(g .* xs)   # ≈ xbar

1.0

The keyword form stationary_distribution(X; δ = δ, rebirth = ψ) computes the stationary distribution when agents die at rate $\delta$ and are reborn with distribution $\psi$ — the resolvent $(\delta I - \mathbb{A}')^{-1} \delta \psi$ — the relevant object in perpetual-youth and firm-entry models:

ψ = zeros(length(xs)); ψ[1] = 1.0    # everyone reborn at the bottom
gδ = stationary_distribution(X; δ = 0.05, rebirth = ψ)
sum(gδ .* xs)

0.7688024920634218

Expectations: Feynman–Kac

feynman_kac(X, ts; f, ψ, v, direction) computes conditional expectations of the general form

\[E\left[\int_0^T e^{-\int_0^t v(x_s) ds} f(x_t) \, dt + e^{-\int_0^T v(x_s) ds} \psi(x_T) \,\Big|\, x_0\right]\]

by implicit Euler time stepping — forecasts (ψ), flow payoffs and present values (f), state-dependent discounting or hazard rates (v). See the expectations tutorial for the computation by hand and the exact semantics of each keyword:

u = feynman_kac(X, 0:0.1:10; ψ = collect(xs))   # E[x_T | x_t], one column per date
maximum(abs, u[:, 1] - (1.0 .+ exp(-0.1 * 10.0) .* (xs .- 1.0)))

0.005733814684626637

Additive functionals: long-run CGFs and tail indices

An AdditiveFunctional(X, μm, σm) represents a cumulative quantity $dm_t = \mu_m(x_t) dt + \sigma_m(x_t) dZ^m_t$ driven by the state — think $m = \log w$ for a size $w$ growing at a state-dependent rate; X can be any process in the package. cgf(m, ξ) returns the long-run scaled cumulant generating function $\Lambda(\xi) = \lim_{t\to\infty} \log E[e^{\xi m_t}]/t$ — the principal eigenvalue of the tilted generator tilted_generator(m, ξ) — and tail_index(m; δ) returns the Pareto exponent of the stationary distribution of $e^m$ under death rate $\delta$. See Computing tail indices:

m = AdditiveFunctional(X, collect(xs .- 1.0), 0.1 .* ones(length(xs)))
tail_index(m; δ = 0.05)

0.5928746121978757

Finite differences

The finite-difference operators the package is built on are also exported, as lazy arrays whose entries are computed on demand. In one dimension, the grid is a vector and the function values a vector of the same length; in multiple dimensions, the grid is a tuple of axis vectors such as (xs, ys), dimensions are selected by number, and the function values are an array of matching size. The direction keyword accepts :forward/:backward and the EconPDEs-style synonyms :up/:down interchangeably:

f = sin.(xs)
FirstDerivative(xs, f; direction = :forward, bc = (0, 0))
FirstDerivative(xs, f; direction = :backward, bc = (0, 0))
SecondDerivative(xs, f; bc = (0, 0))

ys = range(0, 2, length = 40)
F = [x * y for x in xs, y in ys]
grid = (xs, ys)
FirstDerivative(grid, F, 1; direction = :forward)
SecondDerivative(grid, F, 1, 1)
SecondDerivative(grid, F, 1, 2; direction = :up)

The argument bc of FirstDerivative and SecondDerivative is the value of the first derivative at each limit of the grid. It defaults to zero, the right condition for reflecting boundaries — and the hook through which HJB boundary conditions like borrowing constraints enter (see Application to a consumption-saving problem). The convention is deliberately the same as the bc keyword of EconPDEs.jl's pdesolve (the outward first derivative at each boundary), so boundary conditions carry over unchanged between the two packages.