Library

ContinuousTimeProcess{T}

Types inheriting from the abstract type ContinuousTimeProcess{T} characterize the properties of a T-valued stochastic process, play a similar role as distribution types like Exponential in the package Distributions.

SamplePath{T} <: AbstractPath{T}

The struct

struct SamplePath{T}
    tt::Vector{Float64}
    yy::Vector{T}
    SamplePath{T}(tt, yy) where {T} = new(tt, yy)
end

serves as container for discretely observed ContinuousTimeProcesses and for the sample path returned by direct and approximate samplers. tt is the vector of the grid points of the observation/simulation and yy is the corresponding vector of states.

It supports getindex, setindex!, length, copy, vcat, start, next, done, endof.

Like VSamplePath, but with assumptions on tt and dimensionality. Planned replacement for VSamplePath

valtype(::ContinuousTimeProcess) -> T

Returns statespace (type) of a ContinuousTimeProcess{T].

outertype(P::ContinuousTimeProcess) -> T

Returns the type of outer(x), where x is a state of P

ODESolver

Abstract (super-)type for solving methods for ordinary differential equations.

solve!(method, F, X::SamplePath, x0, P) -> X, [err]
solve!(method, X::SamplePath, x0, F) -> X, [err]

Solve ordinary differential equation $(d/dx) x(t) = F(t, x(t))$ or $(d/dx) x(t) = F(t, x(t), P)$ on the fixed grid X.tt writing into X.yy .

method::R3 - using a non-adaptive Ralston (1965) update (order 3).

method::BS3 use non-adaptive Bogacki–Shampine method to give error estimate.

Call _solve! to inline. "Pretty fast if x is a bitstype or a StaticArray."

solve!(::EulerMaruyama, Y, u, W, P) -> X

Solve stochastic differential equation $dX_t = b(t,X_t)dt + σ(t,X_t)dW_t$ using the Euler-Maruyama scheme in place.

R3

Ralston (1965) update (order 3 step of the Bogacki–Shampine 1989 method) to solve $y(t + dt) - y(t) = \int_t^{t+dt} F(s, y(s)) ds$.

BS3

Ralston (1965) update (order 3 step of the Bogacki–Shampine 1989 method) to solve $y(t + dt) - y(t) = \int_t^{t+dt} F(s, y(s)) ds$. Uses Bogacki–Shampine method to give error estimate.

LeftRule <: QuadratureRule

Integrate using left Riemann sum approximation.

fundamental_matrix(tt, P)

Compute fundamental solution.

a(t, x, P::ProcessOrCoefficients)

Fallback for a(t, x, P) calling σ(t, x, P)*σ(t, x, P)'.

sample(tt, P, x1=zero(T))

Sample the process P on the grid tt exactly from its transitionprob(-ability) starting in x1.

sample(::Thinning, T, P::InhomogPoisson) -> tt

sample!(X, P, x1=zero(T))

Sample the process P on the grid X.tt exactly from its transitionprob(-ability) starting in x1 writing into X.yy.

quvar(X)

Computes the (realized) quadratic variation of the path X.

bracket(X)
bracket(X,Y)

Computes quadratic variation process of X (of X and Y).

ito(Y, X)

Integrate a stochastic process Y with respect to a stochastic differential dX.

girsanov(X::SamplePath, P::ContinuousTimeProcess, Pt::ContinuousTimeProcess)

Girsanov log likelihood $\mathrm{d}P/\mathrm{d}Pt(X)$

lp(s, x, t, y, P)

Log-transition density, shorthand for logpdf(transitionprob(s,x,t,P),y).

llikelihood(X::SamplePath, P::ContinuousTimeProcess)

Log-likelihood of observations X using transition density lp.

llikelihood(X::SamplePath, Pº::LocalGammaProcess, P::LocalGammaProcess)

Log-likelihood dPº/dP. (Up to proportionality.)

llikelihood(X::SamplePath, P::LocalGammaProcess)

Bridge log-likelihood with respect to reference measure P.P. (Up to proportionality.)

solve(method::SDESolver, u, W::SamplePath, P) -> X
solve(method::SDESolver, u, W::SamplePath, (b, σ)) -> X

Solve stochastic differential equation $dX_t = b(t,X_t)dt + σ(t,X_t)dW_t$ using method in place.

Example

solve(EulerMaruyama(), 1.0, sample(0:0.1:10, Wiener()), ((t,x)->-x, (t,x)->I))

import Bridge: b, σ
struct OU <: ContinuousTimeProcess{Float64}
    μ::Float64
end
Bridge.b(s, x, P::OU) = -P.μ*x
Bridge.σ(s, x, P::OU) = I

solve(EulerMaruyama(), 1.0, sample(0:0.1:10, Wiener()), OU(1.4))

solve(method::SDESolver, u, W::VSamplePath, P) -> X

Solve stochastic differential equation $dX_t = b(t,X_t)dt + σ(t,X_t)dW_t$ using method.

euler(u, W, P) -> X

Solve stochastic differential equation $dX_t = b(t,X_t)dt + σ(t,X_t)dW_t$ using the Euler scheme.

euler!(Y, u, W, P) -> X

Solve stochastic differential equation $dX_t = b(t,X_t)dt + σ(t,X_t)dW_t$ using the Euler scheme in place.

EulerMaruyama() <: SDESolver

Euler-Maruyama scheme. Euler is defined as alias.

EulerMaruyama() <: SDESolver

Euler-Maruyama scheme. Euler is defined as alias.

StochasticRungeKutta() <: SDESolver

Stochastic Runge-Kutta scheme for T<:Number-valued processes.

StochasticHeun() <: SDESolver

Stochastic heun scheme.

NoDrift(tt, B, β, a)

R3!

Inplace Ralston (1965) update (order 3 step of the Bogacki–Shampine 1989 method) to solve $y(t + dt) - y(t) = \int_t^{t+dt} F(s, y(s)) ds$.

σ!(t, y, Δw, tmp2, P)

Compute stochastic increment at y, $σ Δw$, modifying tmp2.

b!(t, y, tmp1, P)

Compute drift $b$ in y (without factor $Δt$, modifying tmp1.

GammaProcess

A GammaProcess with jump rate γ and inverse jump size λ has increments Gamma(t*γ, 1/λ) and Levy measure

Here Gamma(α,θ) is the Gamma distribution in julia's parametrization with shape parameter α and scale θ.

GammaBridge(t, v, P)

A GammaProcess P conditional on htting v at time t.

ExpCounting(λ)

Counting process with arrival times arrival(P) = Exponential(1/λ) and unit jumps.

CompoundPoisson{T} <: LevyProcess{T}

Abstract type. For a compound Poisson process define rjumpsize(P) -> T and arrival(P) -> Distribution.

 nu(k, P)

(Bin-wise) integral of the Levy measure $\nu(B_k)$ (sic).

uniformthinning!(X, P::GammaProcess, γᵒ)

Return a Gamma process Y with new intensity γᵒ, such that X-Y has intensity γ-γᵒ and Y and X-Y are independent. In the limit $dt \to \infty$ the new Gamma process has each of is jump removed with probability γᵒ/γ. Overwrites X with Y.

ThinningAlg(λmax)

Sampling method for InhomogPoisson by the 'thinning' algorithm.

Examples:

sample(Thinning(λmax), T, InhomogPoisson(λ))

InhomogPoisson(λ)

Inhomogenous Poisson process with intensity function λ(t). See also ThinningAlg.

Bessel3Bridge(t, v, σ)

Bessel(3) bridge from below or above to the point v at time t, not crossing v, with dispersion σ.

BesselProp

Bessel type proposal

endpoint!(X::SamplePath, v)

Convenience functions setting the endpoint of X tov`.

inner(x[, y])

Short-hand for quadratic form x'x (or x'y).

cumsum0(x)

Cumulative sum starting at 0 such that cumsum0(diff(x)) ≈ x.

mat(yy::Vector{SVector})

Reinterpret X or yy to an array without change in memory.

outer(x[, y])

Short-hand for quadratic form xx' (or xy').

CSpline(s, t, x, y = x, m0 = (y-x)/(t-s), m1 = (y-x)/(t-s))

Cubic spline parametrized by $f(s) = x$ and $f(t) = y$, $f'(s) = m_0$, $f'(t) = m_1$.

integrate(cs::CSpline, s, t)

Integrate the cubic spline from s to t.

logpdfnormal(x, Σ)

logpdf of centered Gaussian with covariance Σ

runmean(x)

Running mean of the vector x.

PSD{T}

Simple wrapper for the lower triangular Cholesky root of a positive (semi-)definite element σ.

Gaussian(μ, Σ) -> P

Gaussian distribution with mean μand covarianceΣ. Definesrand(P)and(log-)pdf(P, x). Designed to work withNumbers,UniformScalings,StaticArraysandPSD`-matrices.

Implementation details: On Σ the functions logdet, whiten and unwhiten (or chol as fallback for the latter two) are called.

refine(tt, n)

Refine range by decreasing stepsize by a factor n.

quaternion(m::SMatrix{3,3})

Compute the (rotation-) quarternion of a 3x3 rotation matrix. Useful to create isodensity ellipses from spheres in GL visualizations.

_viridis

Color data of the Viridis map by Nathaniel J. Smith, Stefan van Der Walt, Eric Firing from https://github.com/BIDS/colormap/blob/master/colormaps.py .

mcstart(x) -> state

Create state for random chain online statitics. The entries/value of x are ignored

mcnext(state, x) -> state

Update random chain online statistics when new chain value x was observed. Return new state.

mcband(mc)

Compute marginal 95% coverage interval for the chain from normal approximation.

mcmeanband(mc)

Compute marginal confidence interval for the chain mean using normal approximation

mcstats(mc)

Compute mean and covariance estimates.

mcmarginalstats(mcstates) -> mean, std

Compute meanand marginal standard deviationsstd` for 2d plots.

LinPro(B, μ::T, σ)

Linear diffusion $dX = B(X - μ)dt + σdW$.

Ptilde(cs::CSpline, σ)

Affine diffusion $dX = cs(t) dt + σdW$ with cs a cubic spline ::CSpline.

LinearNoiseAppr(tt, P, x, a, direction = forward)

Precursor of the linear noise approximation of P. For now no attempt is taken to add in a linearization around the deterministic path. direction can be one of :forward, :backward or :nothing. The latter corresponds to choosing β == 0.

LinearAppr(tt, B, β, a)

LinProBridge

Bridge process of P::LinPro with μ == 0 conditional on ending in v at time t.

GuidedProp

General bridge proposal process, only assuming that Pt defines H and r in the right way.

GuidedBridge

Guided proposal process for diffusion bridge using backward recursion.

GuidedBridge(tt, P, Pt, v)

Constructor of guided proposal process for diffusion bridge of P to v on the time grid tt using guiding term derived from linear process Pt.

GuidedPBridge(tt, P, Pt, V, H♢)

Guided proposal process for diffusion bridge of P to v on the time grid tt using guiding term derived from linear process Pt. Initialize using Bridge.gpupdate(H♢, V, L, Σ, v)

BridgePre() <: SDESolver

Precomputed Euler-Maruyama scheme for bridges using bi.

BridgeProp(Target::ContinuousTimeProcess, tt, v, a, cs)

Simple bridge proposal derived from a linear process with time dependent drift given by a CSpline and constant diffusion coefficient a.

Mdb() <: SDESolver

Euler scheme with the diffusion coefficient correction of the modified diffusion bridge.

bridge(method, W, P) -> Y

Integrate with method, where P is a bridge proposal.

Examples

cs = Bridge.CSpline(tt[1], tt[end], Bridge.b(tt[1], v[1], P),  Bridge.b(tt[end], v[2], P))
P° = BridgeProp(Pσ, v), Pσ.a, cs)
W = sample(tt, Wiener())
bridge(BridgePre(), W, P°)

bridge!(method, Y, W, P) -> Y

Integrate with method, where P is a bridge proposal overwritingY`.

Vs(s, T1, T2, v, P)

Time changed V for generation of U.

gpV!(K::SamplePath, P, KT=zero(T))

Precompute V from $(d/dt)V = BV + β$, $V_T = v$ for a guided proposal.

r(t, x, T, v, P)

Returns $r(t,x) = \operatorname{grad}_x \log p(t,x; T, v)$ where $p$ is the transition density of the process $P$.

gpHinv!(K::SamplePath, P, KT=zero(T))

Precompute $K = H^{-1}$ from $(d/dt)K = BK + KB' + a$ for a guided proposal.

gpupdate(H♢, V, L, Σ, v)
gpupdate(P, L, Σ, v)

Return updated H♢, V when observation v at time zero with error Σ is observed.

LocalGammaProcess

compensator0(kstart, P::LocalGammaProcess)

Compensator of GammaProcess approximating the LocalGammaProcess. For kstart == 1 (only choice) this is $\nu_0([b_1,\infty)$.

compensator(kstart, P::LocalGammaProcess)

Compensator of LocalGammaProcess For kstart = 1, this is $\sum_{k=1}^N \nu(B_k)$, for kstart = 0, this is $\sum_{k=0}^N \nu(B_k) - C$ (where $C$ is a constant).

θ(x, P::LocalGammaProcess)

Inverse jump size compared to gamma process with same alpha and beta.

soft(t, T1, T2)

Time change mapping s in [T1, T2](U-time) tot`in[T1, T2](X`-time).

tofs(s, T1, T2)

Time change mapping t in [T1, T2] ($X$-time) to s in [T1, T2] (U-time).

dotVs (s, T1, T2, v, P)

Time changed time derivative of V for generation of U.

SDESolver

Abstract (super-)type for solving methods for stochastic differential equations.

Increments{S<:AbstractPath{T}}

Iterator over the increments of an AbstractPath. Iterates over (i, tt[i], tt[i+1]-tt[i], yy[i+1]-y[i]).

sizedtype(x) -> T

Return an extended type which preserves size information. Makes one(T) and zero(T) for vectors possible.

piecewise(X::SamplePath, [endtime]) -> tt, xx

If X is a jump process with piecewise constant paths and jumps in X.tt, piecewise returns coordinates path for plotting purposes. The second argument allows to choose the right endtime of the last interval.

BridgePre!() <: SDESolver

Precomputed, replacing Euler-Maruyama scheme for bridges using bi.

aeuler(u, s:dtmax:t, P, tau=0.5)

Adaptive Euler-Maruyama scheme from https://arxiv.org/pdf/math/0601029.pdf sampling a path from u at s to t with adaptive stepsize of 2.0^(-k)*dtmax