API

WienerProcess(t0, tf, y0)

Construct a Wiener process on the interval t0 to tf, with initial condition y0.

The noise process noise = WienerProcess(t0, tf, y0) returned by the constructor is a subtype of AbstractNoise{Univariate}.

Sample paths are obtained by populating a pre-allocated vector yt with the sample path, via rand!(rng, noise, yt).

The number of steps for the sample path is determined by the length of the given vector yt, and the time steps are uniform and calculated according to dt = (tf - t0) / (length(yt) - 1). The initial condition is yt[1] = y0, corresponding to the value at time t0.

Since, by definition, $\Delta W_t \sim \mathcal{N}(0, t)$, a sample path is constructed recursively by solving the recursive relation

\[W_{t_i} = W_{t_{i-1}} + \sqrt{\mathrm{dt}} z_i, \qquad i = 1, \ldots,\]

where at each step $z_i$ is drawn from a standard Normal distribution.

OrnsteinUhlenbeckProcess(t0, tf, y0, ν, σ)

Construct an Ornstein Uhlenbeck process $O_t$ on the interval t0 to tf, with initial condition y0, drift -ν and diffusion σ, as defined by the equation

\[\mathrm{d}O_t = -\nu O_t \;\mathrm{d}t + \sigma \;\mathrm{d}W_t.\]

The solution is

\[O_t = e^{-\nu t}O_0 + \sigma \int_0^t e^{-\nu (t - s)}\;\mathrm{d}W_s.\]

The noise process noise = OrnsteinUhlenbeckProcess(t0, tf, y0, ν, σ) returned by the constructor is a subtype of AbstractNoise{Univariate}.

Sample paths are obtained by populating a pre-allocated vector yt with the sample path, via rand!(rng, noise, yt).

The number of steps for the sample path is determined by the length of the given vector yt, and the time steps are uniform and calculated according to dt = (tf - t0) / (length(yt) - 1). The initial condition is yt[1] = y0, corresponding to the value at time t0.

Notice the integral term is a Normal random variable with zero mean and variance

\[\mathbb{E}\left[ \left( \sigma \int_0^t e^{-\nu (t - s)} \;\mathrm{d}W_s\right)^2\right] = \frac{\sigma^2}{2\nu}\left( 1 - e^{-2\nu t} \right).\]

Thus, a sample path is constructed with exact distribution by solving the recursion relation

\[O_{t_i} = e^{-\nu \Delta t} O_{t_{i-1}} + \frac{\sigma}{\sqrt{2\nu}} \sqrt{1 - e^{-2\nu \Delta t}} z_i, \qquad i = 1, \ldots,\]

where at each time step $z_i$ is drawn from a standard Normal distribution.

The Ornstein-Uhlenbeck process has mean, variance, and covariance given by

\[ \mathbb{E}[O_t] = O_0 e^{-\nu t}, \mathrm{Var}[O_t] = \frac{\sigma^2}{2\nu}, \quad \mathrm{Cov}[O_tO_s] = \frac{\sigma^2}{2\nu} e^{-\nu |t - s|}.\]

so that $O_t$ and $O_s$ are significantly correlated only when $|t - s| \lesssim \tau$, where $\tau = 1/\nu$ is a characteristic time scale for the process. When $\tau \rightarrow 0$, i.e. $\nu \rightarrow \infty$, with $\sigma / \nu = \tau\sigma \rightarrow 1$, this approximates a Gaussian white noise.

HomogeneousLinearItoProcess(t0, tf, y0, primitive_a, primitive_bsquare)

Construct a homogeneous linear Itô process noise $Y_t$ on the interval t0 to tf, with initial condition y0, as defined by the equation

\[\mathrm{d}Y_t = a(t) Y_t \;\mathrm{d}t + b(t) Y_t \;\mathrm{d}W_t,\]

provided the primitive of $a=a(t)$ and the primitive of $b^2 = b(t)^2$ are given, via primitive_a and primitive_bsquare, respectively.

The noise process noise = HomogeneousLinearItoProcess(t0, tf, y0, primitive_a, primitive_bsquare) returned by the constructor is a subtype of AbstractNoise{Univariate}.

The exact solution has the form

\[Y_t = y_0 e^{\int_0^t (a(s) - \frac{b(s)^2}{2}) \;\mathrm{d}s + \int_0^t b(s) \;\mathrm{d}W_s}.\]

The basic statistics for this process can be computed by first computing the statistics for its logarithm, which satisfies

\[\mathbb{E}\left[ \ln Y_t \right] = \ln y_0 + \int_0^t (a(s) - \frac{b(s)^2}{2}) \;\mathrm{d}s,\]

and

\[\mathrm{Var}\left( \ln Y_t \right) = \int_0^t b(s)^2 \;\mathrm{d}s.\]

Then, since $\ln Y_t$ is Gaussian, $Y_t$ is log-normal with

\[\mathbb{E}\left[ Y_t \right] = y_0 e^{ \int_0^t a(s) \;\mathrm{d}s},\]

and

\[\mathrm{Var}\left( Y_t \right) = y_0^2 e^{\int_0^t 2a(s) \;\mathrm{d}s}\left( e^{\int_0^t b(s)^2 \;\mathrm{d}s} - 1 \right).\]

A distributionally exact solution is computed on the mesh points in a recursive manner by

\[Y_{t_j} = Y_{t_{j-1}} e^{(p_a(t_j) - p_a(t_{j-1})) - (p_{b^2}(t_j) - p_{b^2}(t_{j-1})/2 + Z_j)}, \qquad j = 1, \ldots,\]

with $Y_0 = y_0$, and where $p_a = p_a(t)$ is the given primitive of $a=a(t)$, $p_{b^2} = p_{b^2}(t)$ is the given primitive of $b^2 = b(t)^2$, and $Z_j \sim \mathcal{N}(0, p_{b^2}(t_j) - p_{b^2}(t_{j-1}))$.

Sample paths are obtained by populating a pre-allocated vector yt with the sample path, via rand!(rng, noise, yt).

The number of steps for the sample path is determined by the length of the given vector yt, and the time steps are uniform and calculated according to dt = (tf - t0) / (length(yt) - 1). The initial condition is yt[1] = y0, corresponding to the value at time t0.

GeometricBrownianMotionProcess(t0, tf, y0, μ, σ)

Construct a Geometric Brownian motion process $Y_t$ on the interval t0 to tf, with initial condition y0, drift μ and diffusion σ, as defined by

\[\mathrm{d}Y_t = \mu Y_t \;\mathrm{d}t + \sigma Y_t \;\mathrm{d}W_t.\]

The noise process noise = GeometricBrownianMotionProcess(t0, tf, y0, μ, σ) returned by the constructor is a subtype of AbstractNoise{Univariate}.

The exact solution is given by

\[Y_t = Y_0 e^{ \left(\mu - \frac{\sigma^2}{2}\right) t + \sigma W_t) }.\]

The discretized solution sample is computed recursively via

\[Y_{t_j} = Y_{t_{j-1}} e^{ \left(\mu - \frac{\sigma^2}{2}\right) \Delta t + \sigma \sqrt{\Delta t} Z_j) },\]

where $Z_t \sim \mathcal{N}(0, 1)$.

Sample paths are obtained by populating a pre-allocated vector yt with the sample path, via rand!(rng, noise, yt).

The number of steps for the sample path is determined by the length of the given vector yt, and the time steps are uniform and calculated according to dt = (tf - t0) / (length(yt) - 1). The initial condition is yt[1] = y0, corresponding to the value at time t0.

CompoundPoissonProcess(t0, tf, λ, dylaw)

Construct a Compound Poisson process on the interval t0 to tf, with point Poisson counter with rate parameter λ and increments given by the distribution dylaw.

The noise process noise = CompoundPoissonProcess(t0, tf, λ, dylaw) returned by the constructor is a subtype of AbstractNoise{Univariate}.

Sample paths are obtained by populating a pre-allocated vector yt with the sample path, via rand!(rng, noise, yt).

The noise process returned by the constructor yields a random sample path of

\[Y_t = \sum_{i=1}^{N_t} \;\mathrm{d}Y_i,\]

where $N_t$ is the number of events up to time $t$.

Then, based on the number n of events, the increment is performed by adding n samples of the given increment distribution dylaw.

The number of steps for the sample path is determined by the length of the given vector yt, and the time steps are uniform and calculated according to dt = (tf - t0) / (length(yt) - 1). The initial condition is set to yt[1] = 0, corresponding to the value at time t0.

PoissonStepProcess(t0, tf, λ, steplaw)

Construct a point Poisson process on the interval t0 to tf, with a point Poisson counter with rate parameter λ and step values given by the distribution steplaw.

The noise process noise = PoissonStepProcess(t0, tf, λ, steplaw) returned by the constructor is a subtype of AbstractNoise{Univariate}.

Sample paths are obtained by populating a pre-allocated vector yt with the sample path, via rand!(rng, noise, yt).

The noise returned by the constructor yields a random sample path of $Y_t = Y_{N_t}$ obtained by first drawing the number n of events between consecutive times with interval dt according to the Poisson distribution n = N(t+dt) - N(t) = Poisson(λdt).

Then, based on the number n of events, the next state is repeated from the previous value, if n is zero, or set to a new sample value of Y, if n is positive. Since it is not cumulative and it has the Markov property, it doesn't make any difference, for the discretized sample, whether n is larger than 1 or not.

The number of steps for the sample path is determined by the length of the given vector yt, and the time steps are uniform and calculated according to dt = (tf - t0) / (length(yt) - 1). The initial condition is yt[1] = y0, corresponding to the value at time t0.

ExponentialHawkesProcess(t0, tf, λ₀, a, δ, dylaw)

Construct an Exponentially Decaying Hawkes process on the interval t0 to tf, with point Poisson counter with rate parameter λ, jump increments given by the distribution dylaw, and exponential decay with rate δ.

An exponentially decaying Hawkes process is a self-exciting point process $\lambda_t$, representing a time-dependent intensity rate for an inhomogenous Poisson counter with an initial intensity $\lambda_0$, a reversion level $a$ with $\lambda_0 \geq a \geq 0$, an exponential decay with rate $\delta > 0$, and positive stationary random jump increments $S_k$, at each arrival time $T_k$. The process is define by

\[ \lambda_t = a + (\lambda_0 - a) e^{-\delta (t-t_0)} + \sum_{t_0 \leq T_k < t} S_k e^{-\delta (t - T_k)}, \quad t \geq t_0.\]

The noise process noise = ExponentialHawkesProcess(t0, tf, λ, δ, dylaw) returned by the constructor is a subtype of AbstractNoise{Univariate}.

Sample paths are obtained by populating a pre-allocated vector yt with the sample path, via rand!(rng, noise, yt).

The noise returned by the constructor yields a random sample path by first drawing the interarrival times, along with the increments given by dylaw, during each mesh time interval, and then applying the exponential decay.

This implementation of the Hawkes process follows A. Dassius and H. Zhao, Exact simulation of Hawkes process with exponentially decaying intensity, Electron. Commun. Probab. 18 (2013), no. 62, 1-13.

TransportProcess(t0, tf, ylaw, f, d)

Construct a transport process on the time interval t0 to tf, with function f=f(t, y) where y is a random vector with dimension d and distribution law for each coordinate given by ylaw.

The noise process noise = TransportProcess(t0, tf, ylaw, f, d) returned by the constructor is a subtype of AbstractNoise{Univariate}.

Sample paths are obtained by populating a pre-allocated vector yt with the sample path, via rand!(rng, noise, yt).

Each random sample path is obtained by first drawing d realizations of the distribution ylaw to build the sample value y and then defining the sample path by Y_t = f(t, y) for each t in the time mesh.

The number of steps for the sample path is determined by the length of the given vector yt, and the time steps are uniform and calculated according to dt = (tf - t0) / (length(yt) - 1). The initial condition is yt[1] = y0, corresponding to the value at time t0.

FractionalBrownianMotionProcess(t0, tf, y0, hurst, d; flags=FFTW.MEASURE)

Construct a fractional Brownian motion process on the interval t0 to tf, with initial condition y0, Hurst parameter hurst and length up to d.

The noise process noise = FractionalBrownianMotionProcess(t0, tf, y0, hurst, d; flags=FFTW.MEASURE) returned by the constructor is a subtype of AbstractNoise{Univariate}.

Sample paths are obtained by populating a pre-allocated vector yt with the sample path, via rand!(rng, noise, yt).

The number of steps for the sample path is determined by the length of the given vector yt, and the time steps are uniform and calculated according to dt = (tf - t0) / (length(yt) - 1). The initial condition is yt[1] = y0, corresponding to the value at time t0. The length of yt must be smaller than or equal to the length d given in the constructor and used for the pre-allocation of the auxiliary vectors.

The method implemented is the one developed by Davies and Harte and uses an FFT transform to drop the order of complexity to O(N log N). For the transform, we use FFTW.jl, and use the flag flags=FFTW.MEASURE for generating the plans. Other common flags can be passed instead.

This implementation of fractional Brownian motion via Davies-Harte method follows Dieker, T. (2004) Simulation of Fractional Brownian Motion. MSc Theses, University of Twente, Amsterdam and A. B. Dieker and M. Mandjes, On spectral simulation of fractional Brownian motion, Probability in the Engineering and Informational Sciences, 17 (2003), 417-434

ProductProcess(noises...)

Construct a multivariate process from independent univariate processes.

The noise process noise = ProductProcess(noises...) returned by the constructor is a subtype of AbstractNoise{Multivariate}.

Sample paths are obtained by populating a pre-allocated matrix yt with the sample path, via rand!(rng, noise, yt).

The number of steps for the sample path is determined by the number of rows of the given matrix yt, and the time steps are uniform and calculated according to dt = (tf - t0) / (size(yt, 1) - 1).

Each columns of yt is populated with a sample path from each univariate process in noise.

UnivariateProcess{T}

Supertype for univariate noise processes.

Alias for AbstractProcess{T, Univariate}.

MultivariateProcess{T}

Supertype for multivariate noise processes.

Alias for AbstractProcess{T, Multivariate}.

rand!(rng::AbstractRNG, noise::AbstractProcess{T}, yt::VecOrMat{T})

Generate sample paths of the noise process.

Populate the vector or matrix yt with a sample path of the process noise, with random numbers generated from rng. See each noise type for details.

RandomEuler(T::DataType=Float64, n::Int=0)

Instantiate a RandomEuler method including a cache vector of length n for a non-allocating solver via the Euler method, solved by solve!(xt, t0, tf, x0, f, yt, params, ::RandomEuler))]

Set n to 0 when solving a scalar equation and set n to the length of the initial vector when solving a system of equations.

RandomHeun(T::DataType=Float64, n::Int=0)

Instantiate a RandomHeun method including three cache vectors of length n for a non-allocating solver via the Heun method, solved by solve!(xt, t0, tf, x0, f, yt, params, ::RandomHeun)).

Set n to 0 when solving a scalar equation and set n to the length of the initial vector when solving a system of equations.

CustomMethod{F, P, N} <: RODEMethod{N}

Custom method for solving a Random ODE. It has two fields:

• solver: a function solver(xt, t0, tf, x0, f, params, yt, solver_params) to solve, on the interval t0 to tf, a Random ODE with right hand side f, "noise" sample path yt, initial condition x0, parameters params for the function f, and extra parameters solver_params for the custom solver;
• solver_params: any argument or series of arguments necessary for the custom solver.

Aliases:

• CustomUnivariateMethod{F, P} = CustomMethod{F, P, Univariate}
• CustomMultivariateMethod{F, P} = CustomMethod{F, P, Multivariate}

ConvergenceSuite(t0, tf, x0law, f, noise, params, target, method, ntgt, ns, m, ks)

Gather the data needed for computing the convergence error for a given RODE of the form

\[ \begin{cases} \displaystyle\frac{\mathrm{d}X_t}{\mathrm{d}t} = f(t, X_t, Y_t), & t_0 \leq t \leq t_f \\ X_{t_0} = X_0. \end{cases}\]

The data comprises of the following:

• the initial and final times t0 and tf;
• the univariate or multivariate distribution x0law for the initial condition $X_0$;
• the right-hand-side term f for the equation, either in the out-of-place form f=f(t, x, y, params), for a scalar equation (i.e. with a univariate initial condition x0law), or in the in-place form f=f(dx, t, x, y, params), for a system of equations (i.e. with a multivariate initial condition x0law);
• the univariate or multivariate process noise for the noise term $Y_t$;
• parameters params for the function f;
• the method target to compute the target solution for the error calculation via solve!(xt, t0, tf, x0, f, yt, target), typically EulerMethod() with a much finer resolution with ntgt mesh points or the order of the square of the highest number of mesh points in ns (see below) or a lower resolution CustomMethod() with an exact distribution of the possible solutions conditioned to the already computed noise points;
• the method to approximate the solution, typically the EulerMethod(), also in the form solve!(xt, t0, tf, x0, f, yt, method);
• the number ntgt of mesh points in the fine mesh on which the target solution will be computed;
• the vector ns with a list of numbers of mesh points to compute the approximate solutions;
• the number m of sample paths to be computed for estimating the strong error via Monte Carlo method.
• the range of steps ks to be used in case one approximates a Random PDE with an increasing number of spatial discretization points, so for each n in ns, one uses a range begin:k:end for the points in the spatial discretization, which defaults to k=[1] in the case of a scalar or of a genuine system of RODEs;

Besides these data obtained from the supplied arguments, a few cache vectors or matrices are created:

• a vector or matrix yt to hold the sample paths of the noise on the finest mesh, with length or row-length being ntgt and the shape depending on whether the noise is univariate or multivariate;
• a vector or matrix xt to hold the sample paths of the target solution, on the finest mesh, with length or row-length being ntgt and the shape depending on whether the law for the initial condition being univariate or multivariate;
• a vector or matrix xnt to hold the sample paths of the approximate solution, with length or row-length being the maximum of those in ns and the shape depending on whether the law for the initial condition being univariate or multivariate.

The actual error is obtained by solving a ConvergenceSuite via solve(rng::AbstractRNG, suite::ConvergenceSuite{T}) with a given RNG.

ConvergenceResult{T, S}(deltas::Vector{T}, trajerrors::Matrix{T}, trajstderrs::Matrix{T}, errors::Vector{T}, stderrs::Vector{T}, lc::T, p::T, pmin::T, pmax::T) where {T, S}

Stores the result of solve(rng, suite) with fields

• deltas: the time steps associated with the number of mesh points in the vector suite.ns;
• trajerrors: a matrix where each column corresponds to the strong error, along the trajectory, at each mesh resolution determined by suite.ns, i.e. trajerrors[i, k] is the error at time $t_0 + i \Delta t$, for the time step $\Delta t = (t_f - t_0) / (n - 1)$ associated with the kth element n = suite.ns[k];
• trajstderrs: a matrix with the corresponding standard error for each entry in trajerrors;
• errors: the maximum, along the trajectory, of the trajerrors;
• stderrs: the corresponding standard error for the Monte-Carlo estimate of the strong errors;
• lc: the logarithm $\log(C)$ of the multiplicative constant in the fitted error CΔtᵖ;
• p: the estimated order of the strong convergence;
• pmin and pmax: the 95% confidence interval for p;

solve(rng, suite::ConvergenceSuite)

Compute the strong errors and the order of convergence of the given ConvergenceSuite suite, with the provided RNG seed rng.

The result is returned in the form of a ConvergenceResult.

generate_error_table(result, suite, info)

Generate the markdown table with the data for the strong errors.

This is obteined from result.errors, with time steps result.deltas and lengths suite.ns, and the provided info for the problem, where info is given as a namedtuple with String fields info.equation, info.ic, and info.noise, some of which may be taken from suite.

plot(results::ConvergenceResult)

Plot the convergence estimate in a log-log scale (time step vs strong error).

It is based on the values provided in results, as computed by solve(::ConvergenceSuite).

The plot consists of a scatter plot for the results.errors and a line plot from the fitted errors ≈ C Δtᵖ, where C = exp(lc), with Δt in results.deltas, lc = results.lc, and p = results.p, with appropriate legends.

plot(suite::ConvergenceSuite; ns = suite.ns, xshow=true, yshow=false, noisealpha=nothing, resolution=2^9)

Plot the target solution, the noise and a few sample paths.

Plot the target solution in suite.xt, the noise in suite.yt, and a few sample paths in the interval t0 to tf, with different time steps as defined by the number of mesh points in suite.ns or as given by the keyword ns as a vector of integers with the desired numbers of mesh points.

The noise, the target solution, and the approximations can be displayed or not, according to the keywords xshow, yshow and ns. If any of them is set to false or nothing, then the corresponding series is not showed.

The linealpha for plotting the noise can be changed via keyword noiselpha.

If suite refers to a system of equations (i.e. with x0law as a MultivariateDistribution instead of a UnivariateDistribution, one can choose to display one or more specific coordinates by specifying the keyword xshow in several possible ways, e.g. xshow=2 (for the second coordinate), or xshow=1:3 (for the first to third coordinates as separate series), or even the sum of all the coordinates, with either xshow=:sum, or the Euclidian norm, if xshow=:norm, or in any other way if xshow is a Function acting on each element of x, when x is a scalar, or on each column of x, when x is vector-valued, such as xshow=sqrt, xshow=sum or xshow=x->2x[1] + x[2]/2, etc.).

Similary, if noise is a ProductProcess, onde can choose to display one or more specific noise contituents, or combinations of them, by specifying the keyword yshow in the same way as for xshow just described.

The resolution keyword can be used to set the maximum resolution for the display of paths, which can be originally about a million points, so reducing it to say the default value of 2^9=512 is plenty enough for the visualization and reduces considerably the size of the generated SVG images. When combining multiple plots into one, it helps reducing even more this value, say resolution=2^7 or 2^8 (128 or 256 points). The resolution must be a factor of suite.ntgt and of all the values of suite.ns.

AbstractProcess{T, N}

Abstract super type for every noise process, with parameter N being either Univariate or Multivariate and T being the eltype of the process.

The following aliases are also defined:

• UnivariateProcess{T} = AbstractProcess{T, Univariate}
• MultivariateProcess{T} = AbstractProcess{T, Multivariate}

The parameter types are borrowed from Distributions.Univariate and Distributions.Multivariate.

RODEMethod{N}

Abstract supertype for the methods for solving a Random ODE.

solve!(xt, t0, tf, x0, f, params, yt, method)

Solve a random ODE with the provided method.

More precisely, sove, inplace, a sample path of the (R)ODE

\[ \begin{cases} \displaystyle\frac{\mathrm{d}X_t}{\mathrm{d}t} = f(t, X_t, Y_t), & t_0 \leq t \leq t_f, \\ X_{t_0} = X_0, \end{cases}\]

with the following arguments:

• a function f(t, x, y, params), if x is a scalar, or f(dx, t, x, y, params), if x is a vector;
• a scalar or vector sample initial condition x0;
• a time interval t0 to tf;
• a sample path yt of a "noise", either a vector (for scalar noise) or a matrix (for vectorial noise);
• parameters params for the function f;
• a numerical method, either RandomEuler() for a scalar equation, RandomEuler(n) for an n-dimensional system of equations, or RandomHeun() for a scalar equation.

The values of xt are updated with the computed solution values.

The time step is obtained from the length of the vector xt via dt = (tf - t0) / (lenght(xt) - 1).

The noise yt should be of the same (row) length as xt.

calculate_trajerrors!(rng, trajerrors, trajstderrs, suite)

Calculate the strong error and the standard error of the suite at each time step along the trajectory.

The strong errors are stored in the provided trajerrors matrix, while the standar errors are stored in the provided trajstderrs. All the info is given in the ConvergenceSuite suite. The RNG seed is given in rng.

This method is used when solving a ConvergenceSuite.