In this paper we consider in detail the realization of Runge-Kutta stochastic numerical methods with weak and strong convergence for systems of stochastic diferential equations in Ito form. The algorithm for generating the Wiener stochastic process, the algorithm for approximation of Ito stochastic integrals, and the code generation algorithms for numerical schemes are described. Python and Julia languages are used. The Jinja2 template engine is used for the code generation .
Copyright © 2018 for the individual papers by the papers’ authors. Copying permitted for private and academic purposes. This volume is published and copyrighted by its editors.
In: K. E. Samouylov, L. A. Sevastianov, D. S. Kulyabov (eds.): Selected Papers of the VIII Conference “Information and Telecommunication Technologies and Mathematical Modeling of High-Tech Systems”, Moscow, Russia, 20-Apr-2018, published at http://ceur-ws.org
This article is divided into three sections. In the first section we define the Wiener
stochastic process and describe its implementation in Python and Julia . In the second
section the approximation algorithm of Ito integrals is given. Finally, in the third part
we introduce the details of code generator implementation for stochastic numerical
methods [
The stochastic process (), > 0 is called scalar Wiener process if the following
conditions are true [
The symbol ∼ (0, ) denotes that is normally distributed random variable with expected value E[ ] = = 0 and variance D[ ] = 2 = .
The Wiener process is the model of Brownian motion (random walk). If we consider the process () in time points 0 = 0 < 1 < 2 < . . . < − 1 < when it experiences random additive changes, then directly from the definition of Wiener process follows: (1) = (0) + 0, (2) = (1) + 1, . . . , ( ) = (− 1) + − 1, where ∼ (0, ), ∀ = 0, . . . , − 1.
In the case of a multidimensional stochastic process one has to generate sequences of normally distributed random variables.
2.2.
To simulate a one-dimensional Wiener process, it is necessary to generate normally distributed random numbers 1, . . . , and to construct their cumulative sums 1, 1 + 2, 1 + 2 + 3, . . . 1 + 2 + 3 + . . . − 1, 1 + 2 + 3 + . . . − 1 + .
The result will be the simulated sample path of the Wiener process () or — using a diferent terminology — concrete implementation of the Wiener process.
In the case of a multidimensional random process, sequences of normally distributed random variables should be generated (that is, arrays, each of elements).
We implemented Wiener process generator in Python [
To generate the Wiener process in Python one should use the WienerProcess class. The following code gives an example of this class usage. import stochastic N = 100 T = (0.0, 10.0) W = stochastic.WienerProcess(N=N, time_interval=T) print("Step size: ", W.dt) print("Time points: ", W.t) print("Process iterations: ", W.dx) print("Wiener Process trajectory: ", W.x) print("Intervals numbers: ", len(W.dx)) print("Points number: ", len(W.x))
The class constructor does not have any required arguments. By default, a process is
generated on a time interval [
10 5 t W 0 5 0 2 4
6 t 8 10
Calculation and approximation of multiple Ito integrals of special form the following form:
Here we will not go into the general theory of multiple stochastic Ito integrals, a
reader can refer to the book [
In general, for the construction of numerical schemes with order of convergence greater than = 12 , it is necessary to calculate single, double and triple Ito integrals of (, +1) = (ℎ) =
d ( ), +1 ∫︁ ∫︁ +1 1
∫︁ +1 1 2 ∫︁
∫︁ ∫︁ (, +1) = (ℎ) =
d ( 2)d ( 1), (, +1) =
In the case of a double integral (ℎ), the exact formula takes place only at = : 00(ℎ) =
= 1 2
For the mixed case 0 and 0 in [
ℎ (ℎ) − 1
︂) 1 (ℎ) , where ∼
(0, ℎ) are multidimensional normal distributed random variables.
= 1, . . . , , the book [
(ℎ) = ℎ ∑∞︁ 1 [︃ 2 =1 ︃(
− ℎ + √︃ 2 2 ℎ
+ (ℎ), )︃
− ︃( + √︃ 2 ℎ )︃]︃ , where variables: where ∼ (0, 1), ∼
(0, 1), = 1, . . . , ; = 1, . . . , ∞; = 1, . . . , is the numerical schema number. From the formulas it is seen that in the case of = , we may get the final expression for the ̸= , one has to sum the infinite series , which we mentioned above. In the case of . This algorithm gives the approximation error of order (ℎ2/), where is number of left terms of an infinite series .
In the article [
I(ℎ) =
W
W − ℎ1× + A(ℎ), 2 A(ℎ) = 2 =1 ℎ ∑∞︁ 1 (︁ V(U + √︀2/ℎ
W) − (U + √︀2/ℎ
W)V )︁ , W, V, U are independent normally distributed multidimensional random W = ( 1, 2, . . . , ) ∼ (0× , ℎ1× ),
V = (1, 2, . . . , ) ∼ (0× , 1× ), U = (1 , 2 , . . . , ) ∼ (0× , 1× ). 1 6 If the programming language supports vectored operations with multidimensional arrays, these formulas can provide a benefit to the performance of the program.
Finally, consider a triple integral. In the only numerical scheme in which it occurs, it
is necessary to be able to calculate only the case of identical indexes
=
= . For
this case, [
1 and the maximum grid step ℎ = max {ℎ− 1}1 .
Next, we assume that the grid is uniform, then ℎ = ℎ = const. x is grid function, which approximate the stochastic process x(), so x0 = x(0), x ≈ x() ∀ = 1, . . . , . 4.1.
The simplest numerical method for solving scalar equations and systems of SDEs is the Euler–Maruyama method.
0 = (0), +1 = + (, )ℎ + ∑︁ (, ) . terministic accuracy order, and value denotes the stochastic part approximation 4.2.
Weak stochastic Runge–Kutta-like method with order 1.5 for a scalar 1 = + ∑︁ 1 ( + 0 ℎ, 0 )ℎ + ∑︁ 1 ( + 1 ℎ, 1 )√︀ℎ, +1 = + ∑︁ ( + 0 ℎ, 0 )ℎ+ + ∑︁ (︂ =1 1 1(ℎ) + 2 11(ℎ) √ ℎ + 3 10(ℎ) + 4 111(ℎ) )︂ ℎ ( + 1 ℎ, 1 ), inside the stochastic Ito integrals: 10(ℎ), 1(ℎ), 11(ℎ), 111(ℎ). where , = 1, . . . , ( is the number of method’s stages). In the above numerical scheme, the Wiener stochastic process may be present in implicit way. It is "hidden" 4.3.
Stochastic Runge–Kutta method with strong order = 1.0 for vector
bution characteristics of stochastic process (). The weak numerical method does
not need information about the trajectory of driving Wiener process and random
increments for these methods can be generated on another probability space [
( + 0 ℎ, 0 )ℎ + ∑︁ ∑︁ 0 ( + 1 ℎ, )^, = + ∑︁ 1
( + 0 ℎ, 0 )ℎ + ∑︁ 1 ( + 1 ℎ, )√︀ℎ, = + ∑︁ 2 ̂︀ ( + 0 ℎ, 0 )ℎ + ∑︁ 2 ( + 1 ℎ, ) √ +1 = + ∑︁ ( + 0 , 0 )ℎ + ∑︁ ∑︁ ( + 1 ℎ, )+ In the weak numerical schema ^ are ^ =
(^^ + ⎧ 1 (^^ − ⎪ ⎪ ⎪⎪⎪ 2 ⎪⎨ 1 ⎪ 2 ⎪⎪⎪⎪ 1 ⎪ , Here ^ denotes three point distributed random variable. It means, that ^ may have three values {− √3ℎ, 0, √3ℎ} with probabilities 1/6, 2/3 and 1/6 respectively. ~ denotes two point distributed random variable {− √ℎ, √ℎ} with probabilities 1/2 and 1/2 respectively.
4.5. Automatic code generation for stochastic numerical methods of
To study the calculation errors and the eficiency of diferent stochastic numerical methods, it is necessary to have a universal implementation of such methods. The universality means the possibility to use any stochastic method with a desired strong or weak error by setting its coeficient table. With direct transfer of mathematical formulas to the program code, one need to use about five nested cycles, which extremely reduces performance, since such code does not take into account a large number of zeros in the coeficient tables and arithmetic operations on zero components are still performed, although this is an extra waste of processor time.
One way to achieve versatility and acceptable performance is to generate code for a numerical method step. This approach minimizes the number of arithmetic operations and saves memory, since the zero coeficients of the method do not have to be stored.
We implemented a code generator for the three stochastic numerical methods mentioned above: – scalar method with strong convergence = 1.5, – vector method with strong convergence = 1.0, – vector method with weak convergence of = 2.0.
We use Python to implement the code generator and Jinja2 [
Information about the coeficients of each particular method is stored as a JSON file of the following structure: { "name": "method’s name (the future name of the function)", "description": "method’s short description", "stage": 4, "det_order": "2.0", "stoch_order": "1.5", "A0": [...], "B0": [...], "A1": [...], "B1": [ ["0", "0", "0", "0"], ["1/2", "0", "0", "0"], ["-1", "0", "0", "0"], ["-5", "3", "1/2", "0"] ], "c0": ["0", "3/4", "0", "0"], "c1": ["0", "1/4", "1", "1/4"], "a": ["1/3", "2/3", "0", "0"], "b1": ["-1", "4/3", "2/3", "0"], "b2": ["-1", "4/3", "-1/3", "0"], "b3": ["2", "-4/3", "-2/3", "0"], "b4": ["-2", "5/3", "-2/3", "1"] }
The parameter stage is the number of method’s stages, det_order is the error order of the deterministic part (), stoch_order is the error order of the stochastic part (), name is the name of the method, which will then be used to create the name of the generated function, so it should be written in one word without spaces. All other parameters are the coeficients of the method. In this case, we give the coeficients of the scalar method with strong convergence = 1.5, omitting the coeficients B0 to save text space. It is necessary to note that the values of the coeficients can be specified in the form of rational fractions, for which they should be presented as JSON a0, a1 and strings and enclosed in double quotes.
For internal representation of stochastic numerical methods we created three Python classes: ScalarMethod, StrongVectorMethod and WeakVectorMethod. The implementation of these classes is contained in the file
coefficients_table.py. The constructors of these classes read the JSON file and, based on them, create objects, which can later be used for code generation. The Fraction class from the Python standard library is used to represent rational coeficients. Each class has a method that generates a coeficient check the correctness of the generator. For example, we give below the formula generated automatically based on the data from JSON file for Runge–Kutta method strong_srk1w2 with stages = 3, and 2 dimensioned Wiener process. Nonzero coeficients of the method are as follows:
021 = 1, 121 = 1, 131 = 1, 121 = 1, 131 = − 1, 1 = 1/2, 2 = 1/2, 02 = 1, 12 = 1, 13 = 1, 11 = 1, 22 = 1/2, 32 = − 1/2.
The numerical scheme formulas are quite cumbersome, despite the large number of zeros in the coeficient table: 01 = , 11 = , 21 = , 02 = + ℎ[︁021 (, 01 )]︁, 12 = + ℎ[︁121 (, 01 )]︁ 22 = + ℎ[︁121 (, 01 )]︁ 13 = + ℎ[︁131 (, 01 )]︁ + 1211 (, 11 ) 11(ℎ) √ + 1212 (, 21 ) 21(ℎ) , √ + 1211 (, 11 ) 12(ℎ) √ + 1212 (, 21 ) 22(ℎ) , √ + 1311 (, 11 ) 11(ℎ) √ + 1312 (, 21 ) 21(ℎ) , √ 23 = + ℎ[︁131 (, 01 )]︁ + +1311 (, 11 ) 12(ℎ) √ + 1312 (, 21 ) 22(ℎ) , √ ℎ ℎ ℎ ℎ ℎ ℎ ℎ ℎ [︁ +1 = + ℎ 1 (, 01 ) + 2 ( + 0ℎ, 02 )]︁ + 111(ℎ)1 (, 11 )
2 + 22√︀ℎ1 ( + 1ℎ, 12 ) + 32√︀ℎ1 ( + 1ℎ, 13 ) 2 3
We present the details of the software implementation of the Wiener process generation in Julia and Python, the algorithm of approximation of multiple Ito integrals and the details of the software implementation of the code generator for stochastic numerical methods of Runge–Kutta type. The code generator allows one to get an efective implementation of the methods and making it possible to use any table of coeficients. The source code of the described programs is open and available at the link https: //bitbucket.org/mngev/sde_num_generation.
The work is partially supported by Russian Foundation for Basic Research (RFBR) grants No 16-07-00556. Also the publication was prepared with the support of the “RUDN University Program 5-100”.