InfinitesimalGenerators.jl

InfinitesimalGenerators.jl computes stationary distributions, conditional expectations, and tail indices of continuous-time Markov processes. It works in two steps. The first step is to build the generator matrix: a Markov-chain approximation of the process on a finite grid, representing transition rates between grid points (in particular, rows must sum to zero and off-diagonal elements are non-negative).

The second step is to use this generator to turn the linear PDEs associated with the process into linear algebra. In particular, the package covers the following problems:

The same generator matrices also appear inside finite-difference methods for nonlinear HJB equations. InfinitesimalGenerators.jl focuses on the linear PDEs once a process is known; for nonlinear HJB problems, see EconPDEs.jl and Application to a consumption-saving problem.

Installation

using Pkg
Pkg.add("InfinitesimalGenerators")

The generator matrix

Every tool in the package is built on one object. For a Markov process $x_t$, the infinitesimal generator is the operator

\[(\mathbb{A}f)(x) = \lim_{t \downarrow 0} \frac{E[f(x_t) \mid x_0 = x] - f(x)}{t}.\]

On a discretized state space with $n$ points, $\mathbb{A}$ becomes an $n \times n$ matrix with non-negative off-diagonal entries and rows summing to zero — the generator (transition-rate) matrix of a continuous-time Markov chain. The package constructs this matrix for process objects, including ContinuousTimeMarkovChain, DiffusionProcess, ProductProcess, SwitchingProcess, and MultivariateDiffusionProcess, and the operators built on it — stationary_distribution, feynman_kac, cgf, tail_index — apply to all of them.

All built-in process types subtype ContinuousTimeMarkovProcess{N}, where N is the number of tensor-product state-space axes. A scalar diffusion or a finite-state chain has N == 1; ProductProcess(X, Z) has N == ndims(X) + ndims(Z). The state shape is size(X), and arrays are flattened in Julia's column-major order when applying generator(X).

using InfinitesimalGenerators

# an Ornstein-Uhlenbeck process dx = -0.03 x dt + 0.01 dZ
x = range(-1, 1, length = 100)
X = DiffusionProcess(x, -0.03 .* x, 0.01 .* ones(100))
generator(X)
100×100 LinearAlgebra.Tridiagonal{Float64, Vector{Float64}}:
 -1.60751    1.60751     ⋅        …    ⋅          ⋅          ⋅ 
  0.122512  -1.70002    1.57751        ⋅          ⋅          ⋅ 
   ⋅         0.122513  -1.67003        ⋅          ⋅          ⋅ 
   ⋅          ⋅         0.122513       ⋅          ⋅          ⋅ 
   ⋅          ⋅          ⋅             ⋅          ⋅          ⋅ 
   ⋅          ⋅          ⋅        …    ⋅          ⋅          ⋅ 
   ⋅          ⋅          ⋅             ⋅          ⋅          ⋅ 
   ⋅          ⋅          ⋅             ⋅          ⋅          ⋅ 
   ⋅          ⋅          ⋅             ⋅          ⋅          ⋅ 
   ⋅          ⋅          ⋅             ⋅          ⋅          ⋅ 
  ⋮                               ⋱                        
   ⋅          ⋅          ⋅             ⋅          ⋅          ⋅ 
   ⋅          ⋅          ⋅             ⋅          ⋅          ⋅ 
   ⋅          ⋅          ⋅             ⋅          ⋅          ⋅ 
   ⋅          ⋅          ⋅             ⋅          ⋅          ⋅ 
   ⋅          ⋅          ⋅        …    ⋅          ⋅          ⋅ 
   ⋅          ⋅          ⋅            0.122513    ⋅          ⋅ 
   ⋅          ⋅          ⋅           -1.67003    0.122513    ⋅ 
   ⋅          ⋅          ⋅            1.57751   -1.70002    0.122512
   ⋅          ⋅          ⋅             ⋅         1.60751   -1.60751

Drift terms are discretized by upwinding (forward differences where the drift is positive, backward where it is negative) and boundaries are reflecting, so the discretized operator is always a valid generator (rows sum to zero and off-diagonal elements are non-negative).

Tutorials

Each tutorial first shows how to compute the object by hand from the generator matrix — they are all a few lines of linear algebra — and then introduces the helper function that packages the computation:

  1. Computing expectations (Kolmogorov backward): the Kolmogorov backward equation and feynman_kac.
  2. Computing distributions (Kolmogorov forward): the Kolmogorov forward equation and stationary_distribution.
  3. Application to a consumption-saving problem: how the same generator matrices enter implicit finite-difference methods for nonlinear HJB equations.
  4. Computing tail indices: principal eigenvalues of tilted generators, cgf, and tail_index.
  • EconPDEs.jl solves nonlinear PDEs from optimization problems (HJB equations). InfinitesimalGenerators.jl works with the linear PDEs generated by a known Markov process, for example after an HJB has been solved and the policy-implied law of motion is known.
  • SimpleDifferentialOperators.jl contains more general tools to define operators with different boundary conditions. In contrast, InfinitesimalGenerators always assumes reflecting boundaries.