Function reference

iterated_integrals(W::AbstractVector, h, eps=h^(3/2);
    ito_correction=true,
    error_norm=MaxL2(),
    alg=optimal_algorithm(length(W),h,eps,error_norm)
)

Simulates an approximation of the iterated stochastic integrals $\int_0^h\int_0^sdW_i(t)dW_j(s)$ for all pairs $1\le i, j \le m$ of the given $m$-dimensional Brownian motion with step size h.

Examples

julia> h = 1/2;

julia> W = [1.0, 0.5]
2-element Vector{Float64}:
 1.0
 0.5

julia> diag(iterated_integrals(W, h, h^(3/2))) ≈ 0.5*W.^2 .- 0.5h
true

iterated_integrals(W::AbstractVector, q_12::AbstractVector, h, eps; 
    ito_correction=true,
    error_norm=FrobeniusL2(),
    alg=optimal_algorithm(length(W),q_12,h,eps,error_norm)
)

Simulates an approximation of the iterated stochastic integrals for finite-dimensional approximations of a Q-Wiener process with covariance matrix $Q = Q^\frac{1}{2}*Q^\frac{1}{2}$. Here q_12 is a vector of the eigenvalues of $Q^\frac{1}{2}$; the square root of the covariance matrix. Equivalently these are the square roots of the eigenvalues of $Q$.

Examples

julia> h = 0.01; dim=10; q = [1/k^2 for k=1:dim];

julia> W = √h * sqrt.(q) .* randn(dim);

julia> diag(iterated_integrals(W,sqrt.(q),h,h^(3/2))) ≈ 0.5*W.^2 .- 0.5*h*q
true

iterated_integrals(W::Real, h::Real, eps::Real=0.0; ito_correction=true, kwargs...)

In the case of a scalar Brownian motion the integral can be explicitly calculated as $\int_0^h\int_0^sdW(t)dW(s) = \frac{1}{2}W(h)^2 - \frac{1}{2}h$.

The parameter eps (as well as all additional keyword arguments) has no effect but is available to provide the same interface as the multidimensional version.

iterated_integrals(W::Real, q_12::Real, h::Real, eps::Real; ito_correction=true, kwargs...)

In the case of a scalar Q-Wiener process with (scalar) covariance Q the integral can be explicitly calculated as $\int_0^h\int_0^sdW(t)dW(s) = \frac{1}{2}W(h)^2 - \frac{1}{2}hQ$.

Note that, as in the multidimensional case, the parameter q_12 denotes the square root of the covariance.

The parameter eps (as well as all additional keyword arguments) has no effect but is available to provide the same interface as the multidimensional version.

levyarea(W, n, alg::Fourier)

Simulates an approximation of the iterated Itô-integrals $\int_0^1W_s\otimes dW_s$ of the given $m$-dimensional increment of a Wiener process with step size 1. The parameter $n$ specifies the number of terms in the approximation and thus determines the accuracy. This algorithm is based on a Fourier expansion of the Wiener process. The algorithm needs approximately $m^2+2\cdot m\cdot n$ Float's and $2\cdot m\cdot n$ random numbers. The time complexity is $\mathcal{O}(m^2\cdot n)$.

levyarea(W, n, alg::Milstein)

Simulates an approximation of the iterated Itô-integrals $\int_0^1W_s\otimes dW_s$ of the given $m$-dimensional increment of a Wiener process with step size 1. The parameter $n$ specifies the number of terms in the approximation and thus determines the accuracy. This is an efficient implementation of the algorithm proposed in Milstein, 1994. It is based on a Fourier expansion of the Wiener process. The algorithm needs approximately $m^2+2\cdot m\cdot n$ Float's. The time complexity is $\mathcal{O}(m^2\cdot n)$.

levyarea(W, n, alg::Wiktorsson)

Simulates an approximation of the iterated Itô-integrals $\int_0^1W_s\otimes dW_s$ of the given $m$-dimensional increment of a Wiener process with step size 1. The parameter $n$ specifies the number of terms in the approximation and thus determines the accuracy. This is an efficient implementation of the algorithm proposed in Wiktorsson, 2001. It is based on the Fourier method from Milstein but incorporates an additional tail sum approximation. The algorithm needs approximately $2\cdot m^2+2\cdot m\cdot n+m$ Float's. The time complexity is $\mathcal{O}(m^2\cdot n)$.

levyarea(W, n, alg::MronRoe)

Simulates an approximation of the iterated Itô-integrals $\int_0^1W_s\otimes dW_s$ of the given $m$-dimensional increment of a Wiener process with step size 1. The parameter $n$ specifies the number of terms in the approximation and thus determines the accuracy. This is an efficient implementation of the algorithm proposed in Mrongowius & Rößler, 2021. It is based on the Fourier method from Milstein but incorporates an improved tail sum approximation. The algorithm needs approximately $m^2+2\cdot m\cdot n$ Float's and $1/2m^2+2\cdot m\cdot n + 1/2m$ random numbers. The time complexity is $\mathcal{O}(m^2\cdot n)$.

abstract type AbstractIteratedIntegralAlgorithm end

Abstract type for algorithms for the simulation of iterated integrals.

julia> subtypes(AbstractIteratedIntegralAlgorithm)
4-element Vector{Any}:
 Fourier
 Milstein
 MronRoe
 Wiktorsson

abstract type AbstractErrorNorm end

Abstract type for different kind of errors one might consider. The most important are MaxL2 and FrobeniusL2. These are actually special cases of the MaxLp{p} and SchattenqLp{p,q} norms:

const MaxL2 = MaxLp{2}
const FrobeniusL2 = SchattenqLp{2,2}

All currently defined norms:

julia> subtypes(LevyArea.AbstractErrorNorm)
3-element Vector{Any}:
 LevyArea.LpMax
 LevyArea.MaxLp
 LevyArea.SchattenqLp

optimal_algorithm(dim, stepsize, eps=stepsize^(3/2), norm=MaxL2())
optimal_algorithm(dim, q_12, stepsize, eps, norm=FrobeniusL2())

Returns the optimal algorithm for the given parameters, i.e. the algorithm that needs to simulate the fewest random numbers to achieve the desired precision eps.

Examples

julia> h = 1/128;

julia> optimal_algorithm(10, h, h^(3/2), MaxL2())
MronRoe()

julia> optimal_algorithm(10, 1.0./(1:10).^2, h, h^(3/2), FrobeniusL2())
Milstein()

terms_needed(dim, stepsize, eps, alg, norm)

Returns the number of terms in the approximating sum that is needed to ensure an error in the given norm of at most eps. This depends on the dimension of the Wiener process dim, the current stepsize and the chosen algorithm.

See also: AbstractIteratedIntegralAlgorithm, AbstractErrorNorm

Examples

julia> h = 1/128;

julia> terms_needed(10, h, h^(3/2), Milstein(), MaxL2())
7

Implementation

New algorithms should only have to implement errcoeff and convorder.

terms_needed(dim, q_12, stepsize, eps, alg, norm)

Used for finite-dimensional approximations of a Q-Wiener process with covariance matrix $Q = Q^\frac{1}{2}*Q^\frac{1}{2}$. Here q_12 is a vector of the eigenvalues of $Q^\frac{1}{2}$; the square root of the covariance matrix. Equivalently these are the square roots of the eigenvalues of $Q$.

Examples

julia> h = 1/128;

julia> dim = 10;

julia> q = [1/k^2 for k=1:dim];

julia> terms_needed(dim, sqrt.(q), h, h^(3/2), Milstein(), FrobeniusL2())
9

convorder(alg::AbstractIteratedIntegralAlgorithm)

Returns the convergence order of the algorithm w.r.t. the truncation parameter.

See also: errcoeff

errcoeff(dim, stepsize, alg, norm)
errcoeff(dim, q_12, stepsize, alg, norm)

Returns the coefficient of the truncation parameter in the error estimate. I.e. the error estimate is of the form

\[\lVert I(h)-\tilde{I}^{(p)}(h) \rVert_* ≤ \mathrm{errcoeff}(m,h) \cdot p^{-γ}\]

where the norm is given by norm, the approximation $\tilde{I}^{(p)}$ is calculated using alg and $γ$ is the order of convergence given by convorder(alg).

norv(dim, n, alg::AbstractIteratedIntegralAlgorithm)

Returns the number of random numbers needed to simulate the iterated integrals for a Wiener process of dimension dim and with truncation parameter n.

effective_cost(dim, stepsize, eps, alg, norm)
effective_cost(dim, q_12, stepsize, eps, alg, norm)

Returns the number of random numbers needed to simulate the iterated integrals with the given parameters.

ito_correction!(I, h=1)

Applies the Itô-correction for iterated integrals to I. This amounts to subtracting $\frac{1}{2}h$ from every element on the diagonal.

Example

julia> M = ones(5,5);

julia> LevyArea.ito_correction!(M, 0.5)


julia> M
5×5 Matrix{Float64}:
 0.75  1.0   1.0   1.0   1.0
 1.0   0.75  1.0   1.0   1.0
 1.0   1.0   0.75  1.0   1.0
 1.0   1.0   1.0   0.75  1.0
 1.0   1.0   1.0   1.0   0.75