Iterated Stochastic Integrals in Matlab
This package implements state-of-the-art methods for the simulation of iterated stochastic integrals. These appear e.g. in higher order algorithms for the solution of stochastic (partial) differential equations.
Install the Matlab toolbox file (LevyArea.mltbx
) from the Releases page or through the Matlab Add-On Explorer.
Alternatively you can download the files into a folder and add it to your Matlab path.
You can do this either manually or by calling the LevyArea_setup.m
script.
The main function of the toolbox is the function
iterated_integrals
. It can be called by prepending the
package name levyarea.iterated_integrals
. However, since
this may be cumbersome one can import the used functions once by
>> import levyarea.iterated_integrals
>> import levyarea.optimal_algorithm
and then one can omit the package name by simply calling
iterated_integrals
or optimal_algorithm
,
respectively.
In the following, we assume that the two functions are imported, so that we can always call them directly without the package name.
First we generate a Wiener increment:
>> m = 5; % dimension of Wiener process
>> h = 0.01; % step size or length of time interval
>> err = 0.05; % error bound
>> W = sqrt(h) * randn(m,1); % increment of Wiener process
Here,
The default call uses h^(3/2) as the precision and chooses the best algorithm automatically:
>> II = iterated_integrals(W,h)
If not stated otherwise, the default error criterion is the II
containing a realisation
of the approximate iterated stochastic integrals that correspond
to the given increment
The desired precision can be optionally provided using a third positional argument:
>> II = iterated_integrals(W,h,err)
Again, the software package automatically chooses the optimal algorithm.
To determine which algorithm is chosen by the package without simulating any iterated
stochastic integrals yet, the function optimal_algorithm
can
be used. The arguments to this function are the dimension of the Wiener
process, the step size and the desired precision:
>> alg = optimal_algorithm(m,h,err); % output: 'Fourier'
It is also possible to choose the algorithm directly using a key-value pair.
The value can be
one of 'Fourier'
, 'Milstein'
, 'Wiktorsson'
and 'MronRoe'
:
>> II = iterated_integrals(W,h,'Algorithm','Milstein')
The desired norm for the prescribed error bound can also be
selected using a key-value pair. The accepted values are
'MaxL2'
and 'FrobeniusL2'
for the
>> II = iterated_integrals(W,h,err,'ErrorNorm','FrobeniusL2')
The simulation of numerical solutions to SPDEs often requires
iterated stochastic integrals based on
>> q = 1./(1:m)'.^2; % eigenvalues of covariance operator
>> QW = sqrt(h) * sqrt(q) .* randn(m,1); % Q-Wiener increment
>> IIQ = iterated_integrals(QW,h,err,'QWiener',sqrt(q))
In this case, the function utilizes a
scaling of the iterated stochastic integrals and also adjusts the error
estimates appropriately such that the error bound holds w.r.t.\ the
iterated stochastic integrals
Note that all discussed keyword arguments (key-value pairs) are
optional and can be combined as needed. Additional information
can be found using the help
function:
>> help iterated_integrals
>> help levyarea.optimal_algorithm
Please cite this package and/or the accompanying paper if you found this package useful. Example BibLaTeX code can be found in the CITATION.bib file.
A Julia version of this package is also available under LevyArea.jl.