Brownian Motion and Iterated Integrals on Mathematica

BrownianUtils.m and RoughPaths.m are two Mathematica packages for investigating SDE driven by (fractional) Brownian motion. They accompany our paper “Parameter Estimation for Rough Differential Equations”.

Click RoughPaths_24Nov2010 here to download (v. 24 Nov 2010) all the files in one go.

This package contains RoughPaths.m, BrownianUtils.m, RoughPathsDemo.nb, BrownianUtilsDemo.nb, IteratedIntegral.tm, IteratedIntegralAllPaths.tm.

BrownianUtils.m

This package offers functions for generating Brownian paths, Fractional Brownian paths and solutions of 1D SDE driven by a (fractional) Brownian path.

Select All Code:
npoints = 500;
T = 0.25;
npaths = 10;
paths = BrownianPaths[T, npoints, npaths];
ListLinePlot[paths, DataRange -> {0, T}]

produces 10 Brownian motion paths:

Select All Code:
h = 0.4;  (* Hurst index *)
res = 250; (* resolution *)
T = 0.25; (* final point *)
numberPaths = 10; (* Number of paths *)
paths = FractionalBrownianPaths[T, res, h, numberPaths];
ListLinePlot[paths, DataRange -> {0, T}]

produces 10 fractional Brownian motion paths, with Hurst parameter 0.4:

(The Hurst parameter relates to the correlation between increments and therefore the roughness of the path. With h=1/2, we end up with the regular Brownian motion.)

Select All Code:
1
2
3
4
5
P[x_] := 1 - x;
Q[x_] := 0.2(1-x)^2;
npaths = 20;
solution = Davie[{P, Q}, incrdistr, T, npoints, npaths, Range[npoints]];
ListLinePlot[solution, DataRange -> {0, T}]

produces 20 simulations of the SDE dY_t=(1-Y_t)dt+0.2(1-Y_t)^2dfB_t:

The algorithm follows the procedure by Davie, described in Deya2010.

RoughPaths.m

This package offers functions for manipulating Iterated Integrals, computing stochastic expansions of solution of SDEs driven by rough paths and the estimation of parameters for SDEs driven by rough paths.

The best is to use our step-by-step demo RoughPathDemo.nb and experiment with the tools, but here are some things you can do with it:

Stochastic expansion of an SDE

Select All Code:
1
2
3
4
5
6
7
Clear[a,b]
P[x_] := a (1 - x);
Q[x_] := b x^2;
expansion = Picard[{P, Q}, 0, 3]
expansion2 = expansion^2;
moment1 = Esp[t, expansion]
moment2 = Esp[t, expansion2]

produces the stochastic expansion of dY_t=a(1-Y_t)dX^0_t+bY_t^2dX^1_t, i.e. a linear combination of iterated integrals capturing the statistics of the drivers. In this case, we have 2 drivers, numbered from 0 to 1. 0 is typically reserved for time.

Select All Code:
1
2
3
a j[{0}]-a^2 j[{0,0}]+a^3 j[{0,0,0}]+2 a^2 b j[{0,0,1}]-6 a^3 b j[{0,0,0,1}]-2 a^3 b j[{0,0,1,0}]+6 a^4 b j[{0,0,0,0,1}]
+12 a^3 b^2 j[{0,0,0,1,1}]+4 a^3 b^2 j[{0,0,1,0,1}]-24 a^4 b^2 j[{0,0,0,0,1,1}]-12 a^4 b^2 j[{0,0,0,1,0,1}]
-4 a^4 b^2 j[{0,0,1,0,0,1}]+48 a^4 b^3 j[{0,0,0,0,1,1,1}]+24 a^4 b^3 j[{0,0,0,1,0,1,1}]+8 a^4 b^3 j[{0,0,1,0,0,1,1}]

Those expansions can be used for simulating solutions of the SDEs, given realisations of the 2 drivers X. They can also be used to derive an approximation of the solutions’ moments.

If one assumes the first driver to be time and the second driver to be Brownian motion, the expectations of the iterated integrals j[…] have an analytic form and we can have a closed expression for the moments of the SDE:

Select All Code:
1
2
3
4
5
6
a t - (a^2 t^2)/2 + (a^3 t^3)/6 + 1/4 a^3 b^2 t^4 - 1/10 a^4 b^2 t^5
a^2 t^2 - a^3 t^3 + (7 a^4 t^4)/12 - (a^5 t^5)/6 + 
 7/10 a^4 b^2 t^5 + (a^6 t^6)/36 - 17/20 a^5 b^2 t^6 + 
 191/420 a^6 b^2 t^7 - 11/105 a^7 b^2 t^8 + 21/80 a^6 b^4 t^8 + 
 1/144 a^8 b^2 t^9 - 43/180 a^7 b^4 t^9 + 33/700 a^8 b^4 t^10 + 
 1/50 a^8 b^6 t^11

The code was heavily influenced by Tocino08, with some improvements for speed and memory management. In particular, the shuffle product is done iteratively rather than recursively and the expectation for Brownian motion is computed more directly. A future publication will detail the changes.

Parameter Estimation of SDEs driven by rough paths

Our approach is fully described in our paper "Parameter Estimation for Rough Differential Equations". The idea is to estimate the empirical moments of the solutions of an SDE given a large number of observations at time T, and to identify those moments with the theoretical moments obtained with the stochastic expansions. If the functions used to define the SDE (e.g. drift and variance) are polynomials, the theoretical expansions at a time T are also polynomials in the parameters. Identifying them reduces to solving a system of polynomials.

The demo gives the example of the SDE dY_t=a(1-Y_t)dt+bY_t^2dW_t, with a=1 and b=2 and observed at T=0.25. We first generate some data using 3000 solutions of the SDE, with Davie[] and get the first two empirical moments m1 and m2:

Select All Code:
1
2
3
4
5
6
7
P[x_] := 1 - x;
Q[x_] := 2 x^2;
npaths = 3000;
T=0.25;
solutions =  Davie[{P, Q}, incrdistr, T, npoints, npaths, Range[npoints]];
m1 = Last@Mean[solutions]
m2 = Last@Mean[solutions^2]

Next we compute the theoretical moments through Picard iterations:

Select All Code:
1
2
3
4
5
Clear[a, b, y0];
P[x_] := a (1 - x);
Q[x_] := b x^2;
expansion = Picard[{P, Q}, 0, 3]
expansion2 = expansion^2;

The expansions are written in terms of Iterated Integrals. The expectations of the expansions (i.e. the theoretical moments) can be obtained in closed form in the case W_t is a Brownian motion:

Select All Code:
1
2
moment1 = Esp[t, expansion]
moment2 = Esp[t, expansion2]

The theoretical moments are polynomials in the parameters a and b. In particular, at T=0.25, they are:

Select All Code:
1
2
3
4
5
6
7
8
9
moment1 = Esp[t, expansion]
moment2 = Esp[t, expansion2]
0.25 a - 0.03125 a^2 + 0.00260417 a^3 + 0.000976563 a^3 b^2 - 
 0.0000976563 a^4 b^2
0.0625 a^2 - 0.015625 a^3 + 0.00227865 a^4 - 0.00016276 a^5 + 
 6.78168*10^-6 a^6 + 0.000683594 a^4 b^2 - 0.00020752 a^5 b^2 + 
 0.0000277565 a^6 b^2 - 1.59854*10^-6 a^7 b^2 + 
 2.6491*10^-8 a^8 b^2 + 4.00543*10^-6 a^6 b^4 - 
 9.11289*10^-7 a^7 b^4 + 4.49589*10^-8 a^8 b^4 + 4.76837*10^-9 a^8 b^6

Identifying the theoretical moments (moment1,moment2) and empirical moments (m1,m2),and limiting to non-complex solutions, we get:

Select All Code:
1
2
0.993803, -2.102
0.993803,  2.102

which is close to the real values (1,2). Note that b is only identifiable up to its sign, since b=2 or b=-2 characterize the same SDE.

In the case where a closed form the expectation of the iterated integrals is not available, we can approximate them by averaging a large number of realisations of the drivers and thus get empirical expectations which we plug into the theoretical expression. This is also done in the demo, for the same SDE, but driven by fractional Brownian motion.

It can also be shown that the parameter estimates follow asymptotic normality. This can be seen by repeating the estimation 100 times and plotting the estimates after a suitable centering and scaling:

This set of points has a covariance matrix close to the identity:

Select All Code:
1
2
0.969113	0.0269188
0.0269188	0.951896

Comments are closed.