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.

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, Y, W, P) -> Y

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

Currently only timedependent sigma, as Ito correction is necessary

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!([::Bridge.TransitionProb], 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))

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.

solve(method, W, P) -> Y

Integrate with method, where P is a bridge proposal from startpoint(P)toendpoint(P)`.

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())
solve(Euler(), W, P°)

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, P)

As P, but without drift.

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.

One step for 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$. Starting point is specified by (t,y)

f!(t,y,k) is a function that takes (t,y) and writes the result in k. ws contains 4 copies of the type of y the result is written into out which is of type y

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 GammaProcessP 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).

uniform_thinning!(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(ThinningAlg(λmax), T, InhomogPoisson(λ))

InhomogPoisson(λ)

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

Bessel{N}(σ)

N-dimensional Bessel process with dispersion σ. Sample with

u = 0.0
t = 0:0.1:1
σ = 1.0
sample(u, t, Bridge.Bessel{3}(σ))

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 cholupper 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 .

Supremum norm

posterior(Val{:λ}, P::GammaProcess, U::SamplePath, prior = (0.0, 0.0))

Marginal posterior distribution of parameter λ. Interpretation of conjugate prior is "observed time, observed increment".

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.

stats = map(OnlineStat, (x, θ, η))

map(push!, stats, (x, θ, η))

mean.(stats)
cov.(stats)

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.

GuidedBridge(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)

PartialBridge

Guided proposal process for diffusion bridge using backward recursion. PartialBridge(tt, P, Pt, L, v, Σ) Guided proposal process for a partial diffusion bridge of P to v on the time grid tt using guiding term derived from linear process Pt. Simulate with bridge!.

PartialBridgeνH

Guided proposal process for diffusion bridge using backward recursion.

PartialBridgeνH(tt, P, Pt, L, v,ϵ Σ)

Guided proposal process for a partial diffusion bridge of `P` to `v` on
the time grid `tt` using guiding term derived from linear process `Pt`.

PartialBridgeνH(tt, P, Pt, ν, Hend⁺)

Guided proposal process on the time grid `tt` using guiding term derived from
linear process `Pt` with backwards equation initialized at `ν, Hend⁺`.

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.

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, V, 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.

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

MeanCov(itr)

Iterator interface for online mean and covariance estimates. Iterates are triples mean, λ, cov/λ

Example:

c = Channel{Vector}(1)
m = Bridge.MeanCov(c)
put!(c, rand(5))
u = iterate(m)
put!(c, rand(5))
u = iterate(m, u[2])
close(c)
u[1][1]

m = Bridge.MeanCov(Channel{Vector{Float64}}(1))
u = register!(m, rand(5))
u = register!(m, rand(5), u)
close(m)
u[1][1]

upsample(x, td, t)

If x is piecewise constant with jumps at td, return values of x at times t.

viridis

Map s onto the first maxviri viridis colors

rescale(x, a=>b, u=>v)

Linearly map the interval [a,b] to [u,v].

Example:

rescale.(x, Ref(extrema(x)))