90 UDC 004.4 Automatic Code Generation for Stochastic Runge–Kutta Methods Migran N. Gevorkyan* , Anastasiya V. Demidova* , Anna V. Korolkova* , Dmitry S. Kulyabov*† * Department of Applied Probability and Informatics, Peoples’ Friendship University of Russia (RUDN University), 6 Miklukho-Maklaya str., Moscow, 117198, Russia † Laboratory of Information Technologies Joint Institute for Nuclear Research 6 Joliot-Curie, Dubna, Moscow region, 141980, Russia Email: gevorkyan_mn@rudn.university, demidova_av@rudn.university, korolkova_av@rudn.university, kulyabov_ds@rudn.university In this paper we consider in detail the realization of Runge–Kutta stochastic numerical methods with weak and strong convergence for systems of stochastic differential 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 . Key words and phrases: stochastic differential equations; stochastic numeric methods; automatic code generation; Python; Julia; the template engine. 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 Gevorkyan Migran N. et al. 91 1. Introduction 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 [1, 2] (in Python language with use of Jinja2 template engine) 2. Stochastic Wiener process 2.1. The Definition and Properties The stochastic process 𝑊 (𝑡), 𝑡 > 0 is called scalar Wiener process if the following conditions are true [3, 4]: – P{𝑊 (0) = 0} = 1, or, in other words, 𝑊 (0) = 0 is almost certain; −1 – 𝑊 (𝑡) is the process with independent increments, i.e. {∆𝑊𝑖 }𝑁 0 are independent random variables: ∆𝑊𝐼 = 𝑊 (𝑡𝐼+1 ) − 𝑊 (𝑡𝐼 ) and 0 6 𝑡0 < 𝑡1 < 𝑡2 < . . . < 𝑡𝑁 6 𝑇 ; – ∆𝑊𝑖 = 𝑊 (𝑡𝐼+1 ) − 𝑊 (𝑡𝐼 ) ∼ 𝒩 (0, 𝑡𝐼+1 − 𝑡𝐼 ) where 0 6 𝑡𝐼+1 < 𝑡𝐼 < 𝑡, 𝐼 = 0, 1, . . . , 𝑁 − 1 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. Program generation of Wiener process 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 different 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 [5] and Julia [6]. To generate an array of numbers distributed according to the standard normal distribution in the case of Python, we used the function random.normal from the NumPy [7] library and, in the case of Julia, the built-in randn function. Both functions give qualitative pseudorandom sequences, since their work uses generators of uniformly distributed 92 ITTMM—2018 pseudorandom numbers based on an algorithm called Mersenne’s vortex [8] (Mersenne Twister), and generators of pseudorandom normally distributed numbers use the Box– Mueller transformation [?, 9]. 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 [0, 1], which is divided into 1000 parts (N=1000). Thus, by default, a path consisting of 1001 points with step dt equal to 0.001. An example of Wiener trajectories is shown on fig. 1. 10 5 Wt 0 5 0 2 4 6 8 10 t Figure 1. Multiple Wiener trajectories In the case of Julia, the Wiener process generator is implemented as the composite data type struct """Stochastic Wiener process""" struct WienerProcess <: AbstractStochasticProcess "Number of process steps" N::Int64 "Time interval starting point" t_0::Float64 "Time interval end point" t_N::Float64 "Step size" dt::Float64 "Time points" T::Vector{Float64} "Winer process values" Gevorkyan Migran N. et al. 93 X::Vector{Float64} "Winer process increments dX ~ N(0, dt)" dX::Vector{Float64} end With following contractors WienerProcess(N::Int64, t_0::Float64, t_N::Float64) WienerProcess(N::Int64, dt::Float64) WienerProcess(N::Int64) WienerProcess() 3. Calculation and approximation of multiple Ito integrals of special form Here we will not go into the general theory of multiple stochastic Ito integrals, a reader can refer to the book [4] for additional information. Here we will consider multiple special integrals, which are included in the stochastic numerical schemes. 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 the following form: 𝑡∫︁ 𝑛+1 𝐼 𝛼 (𝑡𝑛 , 𝑡𝑛+1 ) = 𝐼 𝛼 (ℎ𝑛 ) = d𝑊 𝛼 (𝜏 ), 𝑡𝑛 𝑡∫︁ 𝑛+1∫︁𝜏1 𝐼 𝛼𝛽 (𝑡𝑛 , 𝑡𝑛+1 ) = 𝐼 𝛼𝛽 (ℎ𝑛 ) = d𝑊 𝛼 (𝜏2 )d𝑊 𝛽 (𝜏1 ), 𝑡𝑛 𝑡𝑛 𝑡∫︁ 𝑛+1∫︁𝜏1 ∫︁𝜏2 𝐼 𝛼𝛽𝛾 (𝑡𝑛 , 𝑡𝑛+1 ) = 𝐼 𝛼𝛽𝛾 (ℎ𝑛 ) = d𝑊 𝛼 (𝜏3 )d𝑊 𝛽 (𝜏2 )d𝑊 𝛾 (𝜏1 ), 𝑡𝑛 𝑡𝑛 𝑡𝑛 where 𝛼, 𝛽, 𝛾 = 0 . . . , 𝑚 and 𝑊 𝛼 , 𝛼 = 1, . . . , 𝑚 are components of multidimensional Wiener process. In the case of 𝛼, 𝛽, 𝛾 = 0, the increment of d𝑊 0 (𝜏 ) is assumed to be d𝜏 . The problem is to get analytical formulas for these integrals with ∆𝑊𝑛𝐼 = 𝑊 𝐼 (𝑡𝑛+1 )− 𝑊 𝐼 (𝑡𝑛 ) in them. Despite its apparent simplicity, this is not achievable for all possible combinations of indices. Let us consider in the beginning those cases when it is possible to obtain an analytical expression, and then turn to those cases when it is necessary to use an approximating formulas. In the case of a single integral, the problem is trivial and the analytic expression can be obtained for any index 𝛼: 𝐼 0 (ℎ𝑛 ) = ∆𝑡𝑛 = ℎ𝑛 , 𝐼 𝛼 (ℎ𝑛 ) = ∆𝑊𝑛𝛼 , 𝛼 = 1, . . . , 𝑚. In the case of a double integral 𝐼 𝛼𝛽 (ℎ𝑛 ), the exact formula takes place only at 𝛼 = 𝛽: 1 1 1 (︀ 𝐼 00 (ℎ𝑛 ) = ∆𝑡𝑛 = ℎ2𝑛 , 𝐼 𝛼𝛼 (ℎ𝑛 ) = (∆𝑊𝑛𝛼 )2 − ∆𝑡𝑛 , 𝛼 = 1, . . . , 𝑚, )︀ 2 2 2 in other cases, when 𝛼 ̸= 𝛽 it is not possible in the final form to express 𝐼 𝛼𝛽 (ℎ𝑛 ) in terms of ∆𝑊𝑛𝛼 and ∆𝑡𝑛 , so we can only use numerical approximation. For the mixed case 𝐼 0𝛼 and 𝐼 𝛼0 in [10], simple formulas of the following form are given: (︂ )︂ 1 1 𝐼 0𝛼 (ℎ𝑛 ) = ℎ𝑛 𝐼 𝛼 (ℎ𝑛 ) − √ 𝜁 𝛼 (ℎ𝑛 ) , 2 3 94 ITTMM—2018 (︂ )︂ 1 1 𝐼 𝛼0 (ℎ𝑛 ) = ℎ𝑛 𝐼 𝛼 (ℎ𝑛 ) + √ 𝜁 𝛼 (ℎ𝑛 ) , 2 3 𝛼 ∼ 𝒩 (0, ℎ ) are multidimensional normal distributed random variables. where 𝜁𝑛 𝑛 For the general case 𝛼, 𝛽 = 1, . . . , 𝑚, the book [4] provides the following formulas for approximating the double Ito integral 𝐼 𝛼𝛽 : ∆𝑊𝑛𝛼 ∆𝑊𝑛𝛽 − ℎ𝑛 𝛿 𝛼𝛽 𝐼 𝛼𝛽 (ℎ𝑛 ) = + 𝐴𝛼𝛽 (ℎ𝑛 ), 2 √︃ √︃ ∞ [︃ (︃ )︃ (︃ )︃]︃ 𝛼𝛽 ℎ ∑︁ 1 𝛼 𝛽 2 𝛽 𝛽 𝛼 2 𝛼 𝐴 (ℎ𝑛 ) = 𝑉𝑘 𝑈𝑘 + ∆𝑊𝑛 − 𝑉𝑘 𝑈𝑘 + ∆𝑊𝑛 , 2𝜋 𝑘=1 𝑘 ℎ𝑛 ℎ𝑛 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 𝐼 𝛼𝛽 , which we mentioned above. In the case of 𝛼 ̸= 𝛽, one has to sum the infinite series 𝑎𝛼𝛽 . This algorithm gives the approximation error of order 𝑂(ℎ2 /𝑛), where 𝑛 is number of left terms of an infinite series 𝑎𝑖𝑗 . In the article [11] the matrix form of approximating formulas is introduced. Let 1𝑚×𝑚 , 0𝑚×𝑚 be the unit and zero matrices 𝑚 × 𝑚, then ∆W𝑛 ∆W𝑛𝑇 −ℎ 1 𝑛 𝑚×𝑚 I(ℎ𝑛 ) = + A(ℎ𝑛 ), 2 ∞ ℎ ∑︁ 1 (︁ √︀ √︀ )︁ A(ℎ𝑛 ) = V𝑘 (U𝑘 + 2/ℎ𝑛 ∆W𝑛 )𝑇 − (U𝑘 + 2/ℎ𝑛 ∆W𝑛 )V𝑘𝑇 , 2𝜋 𝑘=1 𝑘 where ∆W𝑛 , V𝑘 , U𝑘 are independent normally distributed multidimensional random variables: ∆W𝑛 = (∆𝑊𝑛1 , ∆𝑊𝑛2 , . . . , ∆𝑊𝑛𝑚 )𝑇 ∼ 𝒩 (0𝑚×𝑚 , ℎ𝑛 1𝑚×𝑚 ), V𝑘 = (𝑉𝑘1 , 𝑉𝑘2 , . . . , 𝑉𝑘𝑚 )𝑇 ∼ 𝒩 (0𝑚×𝑚 , 1𝑚×𝑚 ), U𝑘 = (𝑈𝑘1 , 𝑈𝑘2 , . . . , 𝑈𝑘𝑚 )𝑇 ∼ 𝒩 (0𝑚×𝑚 , 1𝑚×𝑚 ). 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, [10] gives the following formula: 1 (︀ 𝛼 1 (︀ 𝐼 𝛼𝛼𝛼 (ℎ𝑛 ) = (𝐼 (ℎ𝑛 ))3 − 3𝐼 0 (ℎ𝑛 )𝐼 𝛼 (ℎ𝑛 ) = (∆𝑊𝑛𝛼 )3 − 3ℎ𝑛 ∆𝑊𝑛𝛼 . )︀ )︀ 6 6 4. Stochastic Runge–Kutta methods Consider the random process x(𝑡) = (𝑥1 (𝑡), . . . , 𝑥𝑑 (𝑡))𝑇 , where x(𝑡) belongs to the functional space L2 (Ω) with the norm ‖ · ‖. We assume that the random process x(𝑡) is a solution for the Ito SDE [3, 4] if: x(𝑡) = f (𝑡, x(𝑡))d𝑡 + G(𝑡, x(𝑡))dW, where W = (𝑊 1 , . . . , 𝑊 𝑚 )𝑇 is the multidimensional Wiener process, known as the driving process for SDE. The function f : [𝑡0 , 𝑇 ] × R𝑑 → R𝑑 is called as the drift vector, Gevorkyan Migran N. et al. 95 and the matrix-valued function G : [𝑡0 , 𝑇 ] × R𝑑 × R𝑚 → R𝑑 × R𝑚 is called the diffusion matrix. The same equation can be rewritten in the indexed form 𝑚 ∑︁ 𝑥𝛼 (𝑡) = 𝑓 𝛼 (𝑡, 𝑥𝛾 (𝑡))d𝑡 + 𝑔𝛽𝛼 (𝑡, 𝑥𝛾 (𝑡))d𝑊 𝛽 , 𝛽=1 where 𝛼, 𝛾 = 1, . . . , 𝑑, 𝛽 = 1, . . . , 𝑚, and 𝑓 𝛼 (𝑡, 𝑥𝛾 (𝑡)) = 𝑓 𝛼 (𝑡, 𝑥1 (𝑡), . . . , 𝑥𝑑 (𝑡)). On the interval [𝑡0 , 𝑇 ], we introduce the grid 𝑡0 < 𝑡1 < . . . < 𝑡𝑁 = 𝑇 with step ℎ𝑛 = 𝑡𝑛+1 − 𝑡𝑛 , where 𝑛 = 0, . . . , 𝑁 − 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. Euler–Maruyama numerical method The simplest numerical method for solving scalar equations and systems of SDEs is the Euler–Maruyama method. 𝑑 ∑︁ 𝑥𝛼 𝛼 𝛼 𝛼 𝛼 𝛼 0 = 𝑥 (𝑡0 ), 𝑥𝑛+1 = 𝑥𝑛 + 𝑓 (𝑡𝑛 , 𝑥𝑛 )ℎ𝑛 + 𝐺𝛼 𝛾 𝛽 𝛽 (𝑡𝑛 , 𝑥𝑛 )∆𝑊𝑛 . 𝛾=1 The method has a strong order (𝑝𝑑 , 𝑝𝑠 ) = (1.0, 0.5). The value 𝑝𝑑 denotes the de- terministic accuracy order, and value 𝑝𝑠 denotes the stochastic part approximation order. 4.2. Weak stochastic Runge–Kutta-like method with order 1.5 for a scalar Wiener process In the case of a scalar SDE it is possible to construct a numerical scheme with strong convergence 𝑝 = 1.5 [12–16]: 𝑠 𝑠 𝐼 10 (ℎ𝑛 ) 𝐴𝑖0𝑗 𝑓 (𝑡𝑛 + 𝑐𝑗0 ℎ𝑛 , 𝑋0𝑗 )ℎ𝑛 + 𝑔(𝑡𝑛 + 𝑐𝑗1 ℎ𝑛 , 𝑋1𝑗 ) √ ∑︁ ∑︁ 𝑋0𝑖 = 𝑥𝑛 + 𝑖 𝐵0𝑗 , 𝑗=1 𝑗=1 ℎ𝑛 𝑠 𝑠 𝐴𝑖1𝑗 𝑓 (𝑡𝑛 + 𝑐𝑗0 ℎ𝑛 , 𝑋0𝑗 )ℎ𝑛 + 𝑔(𝑡𝑛 + 𝑐𝑗1 ℎ𝑛 , 𝑋1𝑗 ) ∑︁ ∑︁ √︀ 𝑋1𝑖 = 𝑥𝑛 + 𝑖 𝐵1𝑗 ℎ𝑛 , 𝑗=1 𝑗=1 𝑠 ∑︁ 𝑥𝑛+1 = 𝑥𝑛 + 𝑎𝑖 𝑓 (𝑡𝑛 + 𝑐𝑖0 ℎ𝑛 , 𝑋0𝑖 )ℎ𝑛 + 𝑖=1 𝑠 (︂ 𝐼 11 (ℎ𝑛 ) 𝐼 10 (ℎ𝑛 ) 𝐼 111 (ℎ𝑛 ) ∑︁ )︂ + 𝑏1𝑖 𝐼 1 (ℎ𝑛 ) + 𝑏2𝑖 √ + 𝑏3𝑖 + 𝑏4𝑖 𝑔(𝑡𝑛 + 𝑐𝑖1 ℎ𝑛 , 𝑋1𝑖 ), 𝑖=1 ℎ𝑛 ℎ 𝑛 ℎ𝑛 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" inside the stochastic Ito integrals: 𝐼 10 (ℎ𝑛 ), 𝐼 1 (ℎ𝑛 ), 𝐼 11 (ℎ𝑛 ), 𝐼 111 (ℎ𝑛 ). 96 ITTMM—2018 4.3. Stochastic Runge–Kutta method with strong order 𝑝 = 1.0 for vector Wiener process For SDE system with a multidimensional Wiener process, one can construct a stochastic numerical Runge-Kutta scheme of strong order 𝑝𝑠 = 1.0 [4, 17] using single and double Ito integrals [18]. 𝑠 𝑠 𝑚 ∑︁ 𝐴𝑖0𝑗 𝑓 𝛼 (𝑡𝑛 + 𝑐𝑗0 ℎ𝑛 , 𝑋 0𝑗𝛽 )ℎ𝑛 + 𝑗 ∑︁ ∑︁ 𝑋 0𝑖𝛼 = 𝑥𝛼 𝑛 + 𝑖 𝐵0𝑗 𝐺𝛼 𝑙 (𝑡𝑛 + 𝑐1 ℎ𝑛 , 𝑋 𝑙𝑗𝛽 𝑙 )𝐼 (ℎ𝑛 ), 𝑗=1 𝑙=1 𝑗=1 𝑠 𝑚 ∑︁ 𝑠 𝑙𝑘 𝑙𝑗𝛽 𝐼 (ℎ𝑛 ) 𝐴𝑖1𝑗 𝑓 𝛼 (𝑡𝑛 + 𝑐𝑗0 ℎ𝑛 , 𝑋 0𝑗𝛽 )ℎ𝑛 + 𝑗 ∑︁ ∑︁ 𝑋 𝑘𝑖𝛼 = 𝑥𝛼 𝑛 + 𝑖 𝐵1𝑗 𝐺𝛼 𝑙 (𝑡𝑛 + 𝑐1 ℎ𝑛 , 𝑋 ) √ , 𝑗=1 𝑙=1 𝑗=1 ℎ𝑛 𝑠 ∑︁ 𝑚 ∑︁ ∑︁ 𝑠 √︀ 𝑥𝛼 𝛼 𝑛+1 = 𝑥𝑛 + 𝑎𝑖 𝑓 𝛼 (𝑡𝑛 + 𝑐𝑖0 ℎ𝑛 , 𝑋 0𝑖𝛽 )ℎ𝑛 + (𝑏1𝑖 𝐼 𝑘 (ℎ𝑛 ) + 𝑏2𝑖 ℎ𝑛 )𝐺𝛼 𝑖 𝑘 (𝑡𝑛 + 𝑐1 ℎ𝑛 , 𝑋 𝑘𝑖𝛽 ), 𝑖=1 𝑘=1 𝑖=1 𝑛 = 0, 1, . . . , 𝑁 − 1; 𝑖 = 1, . . . , 𝑠; 𝛽, 𝑘 = 1, . . . , 𝑚; 𝛼 = 1, . . . , 𝑑. 4.4. Stochastic Runge–Kutta method with weak order 𝑝 = 2.0 for vector Wiener process Numerical methods with weak convergence are good for approximation the distri- 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 [4, 19–21]. 𝑠 𝑠 ∑︁ 𝑚 𝐴𝑖0𝑗 𝑓 𝛼 (𝑡𝑛 + 𝑐𝑗0 ℎ𝑛 , 𝑋 0𝑗𝛽 )ℎ𝑛 + 𝑗 ∑︁ ∑︁ 𝑙𝑗𝛽 ˆ𝑙 𝑋 0𝑖𝛼 = 𝑥𝛼 𝑛 + 𝑖 𝐵0𝑗 𝐺𝛼 𝑙 (𝑡𝑛 + 𝑐1 ℎ𝑛 , 𝑋 )𝐼 , 𝑗=1 𝑗=1 𝑙=1 𝑠 𝑠 𝐴𝑖1𝑗 𝑓 𝛼 (𝑡𝑛 + 𝑐𝑗0 ℎ𝑛 , 𝑋 0𝑗𝛽 )ℎ𝑛 + 𝑗 ∑︁ ∑︁ √︀ 𝑋 𝑘𝑖𝛼 = 𝑥𝛼 𝑛 + 𝑖 𝐵1𝑗 𝐺𝛼 𝑘 (𝑡𝑛 + 𝑐1 ℎ𝑛 , 𝑋 𝑘𝑗𝛽 ) ℎ𝑛 , 𝑗=1 𝑗=1 𝑠 𝑠 𝑚 ˆ𝑘𝑙 𝑙𝑗𝛽 𝐼 𝐴𝑖2𝑗 𝑓 𝛼 (𝑡𝑛 + 𝑐𝑗0 ℎ𝑛 , 𝑋 0𝑗𝛽 )ℎ𝑛 + 𝑗 ∑︁ ∑︁ ∑︁ ̂︀ 𝑘𝑖𝛼 = 𝑥𝛼 + 𝑋 𝐵2𝑗𝑖 𝐺𝛼 )√ , 𝑛 𝑙 (𝑡𝑛 + 𝑐1 ℎ𝑛 , 𝑋 𝑗=1 𝑗=1 𝑙=1,𝑙̸=𝑘 ℎ𝑛 𝑠 𝑠 ∑︁𝑚 (︃ )︃ ∑︁ ∑︁ 𝐼ˆ𝑘𝑘 𝑥𝛼 𝛼 𝑛+1 = 𝑥𝑛 + 𝑎𝑖 𝑓 𝛼 (𝑡𝑛 + 𝑐𝑖0 , 𝑋 0𝑖𝛽 )ℎ𝑛 + 𝑏1𝑖 𝐼ˆ𝑘 + 𝑏2𝑖 √ 𝐺𝛼 𝑖 𝑘 (𝑡𝑛 + 𝑐1 ℎ𝑛 , 𝑋 𝑘𝑖𝛽 )+ 𝑖=1 𝑖=1 𝑘=1 ℎ𝑛 𝑠 ∑︁ ∑︁ 𝑚 (︁ √︀ )︁ + 𝑏3𝑖 𝐼ˆ𝑘 + 𝑏4𝑖 ℎ𝑛 𝐺𝛼 𝑖 ̂︀ 𝑘𝑖𝛽 ) 𝑘 (𝑡𝑛 + 𝑐2 ℎ𝑛 , 𝑋 𝑖=1 𝑘=1 In the weak numerical schema 𝐼ˆ𝑘𝑙 are 1 ˆ𝑘 ˆ𝑙 √ ˜𝑘 ⎧ ⎪ ⎪ (𝐼 𝐼 − ℎ𝑛 𝐼 ), 𝑘 < 𝑙, 2 ⎪ ⎪ ⎪ √ ⎪ ⎨1 𝐼ˆ𝑘𝑙 = (𝐼ˆ𝑘 𝐼ˆ𝑙 + ℎ𝑛 𝐼˜𝑙 ), 𝑙 < 𝑘, ⎪ ⎪ ⎪ 2 ⎪ 1 ⎩ ((𝐼ˆ𝑘 )2 − ℎ𝑛 ). 𝑘 = 𝑙. ⎪ ⎪ 2 Gevorkyan Migran N. et al. 97 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 Runge–Kutta type To study the calculation errors and the efficiency of different 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 coefficient 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 coefficient 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 coefficients of the method do not have to be stored. We implemented a code generator for the three stochastic numerical methods men- tioned 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 [22] template engine. This template engine was originally created to generate HTML code, but its syntax is universal and allows you to generate text of any kind without reference to any programming or markup language. Information about the coefficients 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 98 ITTMM—2018 parameters are the coefficients of the method. In this case, we give the coefficients of the scalar method with strong convergence 𝑝𝑠 = 1.5, omitting the coefficients a0 , a1 and B0 to save text space. It is necessary to note that the values of the coefficients can be specified in the form of rational fractions, for which they should be presented as JSON strings and enclosed in double quotes. For internal representation of stochastic numerical methods we created three Python classes: ScalarMethod, StrongVectorMethod and WeakVectorMethod. The implementa- tion 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 coefficients. Each class has a method that generates a coefficient table in LaTeX format. The file stoch_rk_generator.py is a script which handles the jinja2 templates and, based on them, generates a code of python functions. For vector stochastic methods, a code is generated for dimensions up to 6. Functions are named based on the information specified in JSON files, such as strong_srk1w2, strong_srk2w5, weak_srk2w6, and so on. In addition to the code in Python, LaTeX formulas are generated. It allows one to 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 coefficients of the method are as follows: 𝐴201 = 1, 𝐴211 = 1, 𝐴311 = 1, 𝐵11 2 3 = 1, 𝐵11 = −1, 𝑎1 = 1/2, 𝑎2 = 1/2, 𝑐20 = 1, 𝑐21 = 1, 𝑐31 = 1, 𝑏11 = 1, 𝑏22 = 1/2, 𝑏23 = −1/2. The numerical scheme formulas are quite cumbersome, despite the large number of zeros in the coefficient table: 𝑋 01𝛼 = 𝑥𝛼𝑛, 𝑋 11𝛼 = 𝑥𝛼𝑛, 𝑋 21𝛼 = 𝑥𝛼𝑛, [︁ ]︁ 𝑋 02𝛼 = 𝑥𝛼 2 𝛼 𝑛 + ℎ𝑛 𝐴01 𝑓 (𝑡𝑛 , 𝑋 01𝛽 ) , [︁ ]︁ 𝑋 12𝛼 = 𝑥𝛼 2 𝛼 𝑛 + ℎ𝑛 𝐴11 𝑓 (𝑡𝑛 , 𝑋 01𝛽 ) 11 21 2 11𝛽 𝐼 (ℎ𝑛 ) 21𝛽 𝐼 (ℎ𝑛 ) + 𝐵11 𝐺𝛼1 (𝑡𝑛 , 𝑋 ) √ + 𝐵11 2 𝐺𝛼 2 (𝑡𝑛 , 𝑋 ) √ , ℎ𝑛 ℎ𝑛 [︁ ]︁ 𝑋 22𝛼 = 𝑥𝛼 2 𝛼 𝑛 + ℎ𝑛 𝐴11 𝑓 (𝑡𝑛 , 𝑋 01𝛽 ) 12 22 2 11𝛽 𝐼 (ℎ𝑛 ) 21𝛽 𝐼 (ℎ𝑛 ) + 𝐵11 𝐺𝛼1 (𝑡𝑛 , 𝑋 ) √ + 𝐵11 2 𝐺𝛼 2 (𝑡𝑛 , 𝑋 ) √ , ℎ𝑛 ℎ𝑛 [︁ ]︁ 𝑋 13𝛼 = 𝑥𝛼 3 𝛼 𝑛 + ℎ𝑛 𝐴11 𝑓 (𝑡𝑛 , 𝑋 01𝛽 ) 11 21 3 11𝛽 𝐼 (ℎ𝑛 ) 21𝛽 𝐼 (ℎ𝑛 ) + 𝐵11 𝐺𝛼1 (𝑡𝑛 , 𝑋 ) √ + 𝐵11 3 𝐺𝛼 2 (𝑡𝑛 , 𝑋 ) √ , ℎ𝑛 ℎ𝑛 [︁ ]︁ 𝑋 23𝛼 = 𝑥𝛼 3 𝛼 𝑛 + ℎ𝑛 𝐴11 𝑓 (𝑡𝑛 , 𝑋 01𝛽 ) 12 22 3 11𝛽 𝐼 (ℎ𝑛 ) 21𝛽 𝐼 (ℎ𝑛 ) + +𝐵11 𝐺𝛼 1 (𝑡𝑛 , 𝑋 ) √ 3 + 𝐵11 𝐺𝛼 2 (𝑡𝑛 , 𝑋 ) √ , ℎ𝑛 ℎ𝑛 [︁ ]︁ 𝑥𝛼 𝛼 𝛼 𝑛+1 = 𝑥𝑛 + ℎ𝑛 𝑎1 𝑓 (𝑡𝑛 , 𝑋 01𝛽 ) + 𝑎2 𝑓 𝛼 (𝑡𝑛 + 𝑐20 ℎ𝑛 , 𝑋 02𝛽 ) + 𝑏11 𝐼 1 (ℎ𝑛 )𝐺𝛼 1 (𝑡𝑛 , 𝑋 11𝛽 ) √︀ √︀ 2 𝛼 2 12𝛽 2 𝛼 3 13𝛽 + 𝑏2 ℎ𝑛 𝐺1 (𝑡𝑛 + 𝑐1 ℎ𝑛 , 𝑋 ) + 𝑏3 ℎ𝑛 𝐺1 (𝑡𝑛 + 𝑐1 ℎ𝑛 , 𝑋 ) Gevorkyan Migran N. et al. 99 √︀ + 𝑏11 𝐼 2 (ℎ𝑛 )𝐺𝛼 2 (𝑡𝑛 , 𝑋 21𝛽 ) + 𝑏22 ℎ𝑛 𝐺𝛼 2 2 (𝑡𝑛 + 𝑐1 ℎ𝑛 , 𝑋 22𝛽 ) √︀ + 𝑏3 ℎ𝑛 𝐺2 (𝑡𝑛 + 𝑐31 ℎ𝑛 , 𝑋 23𝛽 ). 2 𝛼 5. Conclusion 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 effective implementation of the methods and making it possible to use any table of coefficients. The source code of the described programs is open and available at the link https: //bitbucket.org/mngev/sde_num_generation. Acknowledgments 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”. References 1. M. N. Gevorkyan, A. V. Demidova, T. R. Velieva, A. V. Korol’kova, D. S. Kulyabov, L. A. Sevast’yanov, Implementing a Method for Stochastization of One-Step Pro- cesses in a Computer Algebra System, Programming and Computer Software 44 (2) (2018) 86–93. arXiv:1805.03190, doi:10.1134/S0361768818020044. 2. A. V. Demidova, M. N. Gevorkyan, D. S. Kulyabov, A. V. Korolkova, L. A. Sevastianov, The Automation of Stochastization Algorithm with Use of SymPy Computer Algebra Library, EPJ Web of Conferences 173 (2018) 05006_1–4. doi:10.1051/epjconf/201817305006. 3. B. Øksendal, Stochastic differential equations. An introduction with applications, 6th Edition, Springer, Berlin Heidelberg New York, 2003. 4. P. E. Kloeden, E. Platen, Numerical Solution of Stochastic Differential Equations, 2nd Edition, Springer, Berlin Heidelberg New York, 1995. 5. G. Rossum, Python reference manual, Tech. rep., Amsterdam, The Netherlands, The Netherlands (1995). 6. J. Bezanson, A. Edelman, S. Karpinski, V. B. Shah, Julia: A Fresh Approach to Numerical Computing, SIAM Review 59 (2017) 65–98. 7. E. Jones, T. Oliphant, P. Peterson, et al., SciPy: Open source scientific tools for Python, [Online; accessed 08.10.2017] (2001–). URL http://www.scipy.org/ 8. M. Matsumoto, T. Nishimura, Mersenne twister: A 623-dimensionally Equidis- tributed Uniform Pseudo-random Number Generator, ACM Trans. Model. Comput. Simul. 8 (1) (1998) 3–30. doi:10.1145/272991.272995. 9. G. E. P. Box, M. E. Muller, A note on the generation of random normal deviates, The Annals of Mathematical Statistics 29 (2) (1958) 610–611. 10. A. Rößler, Runge-Kutta Methods for the Numerical Solution of Stochastic Dif- ferential Equations, Ph.D. thesis, Technischen Universität Darmstadt, Darmstadt (februar 2003). 11. M. Wiktorsson, Joint characteristic function and simultaneous simulation of iterated Itô integrals for multiple independent Brownian motions, The Annals of Applied Probability 11 (2) (2001) 470–487. 12. K. Burrage, P. M. Burrage, High strong order explicit Runge-Kutta methods for stochastic ordinary differential equations, Appl. Numer. Math. (22) (1996) 81–101. 100 ITTMM—2018 13. K. Burrage, P. M. Burrage, J. A. Belward, A bound on the maximum strong order of stochastic Runge-Kutta methods for stochastic ordinary differential equations., BIT (37) (1997) 771–780. 14. K. Burrage, P. M. Burrage, General order conditions for stochastic Runge-Kutta methods for both commuting and non-commuting stochastic ordinary differential equation systems, Appl. Numer. Math. (28) (1998) 161–177. 15. P. M. Burrage, Runge-Kutta Methods for Stochastic Differential Equations, Ph.D. thesis, University of Qeensland, Australia (1999). 16. K. Burrage, P. M. Burrage, Order conditions of stochastic Runge-Kutta methods by B-series, SIAM J. Numer. Anal. (38) (2000) 1626–1646. 17. A. R. Soheili, M. Namjoo, Strong approximation of stochastic differential equations with Runge–Kutta methods, World Journal of Modelling and Simulation 4 (2) (2008) 83–93. 18. A. Rößler, Strong and Weak Approximation Methods for Stochastic Differential Equations — Some Recent Developments (2010). 19. Y. Komori, T. Mitsuri, Stable ROW-Type Weak Scheme for Stochastic Differential Equations, RIMS Kokyuroku (932) (1995) 29–45. 20. V. Mackevičius, Second-order weak approximations for stratonovich stochastic differential equations, Lithuanian Mathematical Journal 34 (2) (1994) 183–200. doi:10.1007/BF02333416. URL http://dx.doi.org/10.1007/BF02333416 21. A. Tocino, R. Ardanuy, Runge–Kutta methods for numerical solution of stochastic differential equations, Journal of Computational and Applied Mathematics (138) (2002) 219–241. 22. Jinja2 official site. URL http://http://jinja.pocoo.org