Performing Parallel Monte Carlo and Moment Equations Methods for Itô and Stratonovich Stochastic Differential Systems: R Package Sim.DiffProc

We introduce Sim.DiffProc, an R package for symbolic and numerical computations on scalar and multivariate systems of stochastic differential equations (SDEs). It provides users with a wide range of tools to simulate, estimate, analyze, and visualize the dynamics of these systems in both forms, Ito and Stratonovich. One of Sim.DiffProc key features is to implement the Monte Carlo method for the iterative evaluation and approximation of an interesting quantity at a fixed time on SDEs with parallel computing, on multiple processors on a single machine or a cluster of computers, which is an important tool to improve capacity and speed-up calculations. We also provide an easy-to-use interface for symbolic calculation and numerical approximation of the first and central second-order moments of SDEs (i.e., mean, variance and covariance), by solving a system of ordinary differential equations, which yields insights into the dynamics of stochastic systems. The final result object of Monte Carlo and moment equations can be derived and presented in terms of LATEX math expressions and visualized in terms of LATEX tables. Furthermore, we illustrate various features of the package by proposing a general bivariate nonlinear dynamic system of Haken-Zwanzig, driven by additive, linear and nonlinear multiplicative noises. In addition, we consider the particular case of a scalar SDE driven by three independent Wiener processes. The Monte Carlo simulation thereof is obtained through a transformation to a system of three equations. We also study some important applications of SDEs in different fields.

Before we give a rigorous definition of SDEs, we show how they arise as randomly perturbed ODEs and give a physical interpretation. If X t is a differentiable function defined for t ≥ 0, f (t, x) is a function of time t and state x, the following relation is satisfied for all t (0 ≤ t ≤ T ): Then X t is a solution of the ODE with the initial condition X 0 = x 0 . The above equation can be written in other forms (by continuity of X t ): Although the Wiener process (Brownian motion) is not differentiable with respect to time, we write the White noise ξ t as its fictitious derivative of the Wiener process W t , i.e., ξ t = dW t /dt = W t , just with the purpose to mimic standard ODE calculus. This will allow us to introduce the Itô stochastic integral in the following form: if g(t, x) is the intensity of the White noise, then we define the Itô (1944) integral as: where the middle term T 0 g(t, X t )W t dt has no mathematical meaning. Indeed, SDEs are just a symbolic writing for the mathematically correct stochastic integral representation.
SDEs arise when the coefficients of (1) are perturbed by White noise. For example, if X t denotes the population density, then the population growth can be described by the ODE: the growth is exponential with birth rate a. Random perturbation of the birth rate results in the equation: or the SDE: There are two widely used types of stochastic calculus: Itô (1944Itô ( , 1946 and Stratonovich (1966), differing with respect to the stochastic integral used. The Itô integral takes the left end point of each subinterval as the evaluation point, whereas the Stratonovich integral uses the midpoint of each subinterval. Under certain conditions, the Itô integral has the most convenient property of being a martingale (see Gikhman and Skorokhod 2007), but has a transformation rule, called the Itô's formula, which is more complicated than the chain rule of deterministic calculus. In contrast, the Stratonovich integral meets the transformation rules of classical calculus (see Kloeden and Platen 1991a,b). Practical considerations typically dictate which version is appropriate. The following classical example demonstrates how the Itô integral differs from Stratonovich 1 integral: Itô stochastic calculus: Stratonovich stochastic calculus: Let X(t) = {X 1 (t), X 2 (t), . . . , X d (t)}, with d ≥ 1 be an d-dimensional random process with components satisfying the Itô SDEs: dX j (t) = f j (t, X(t))dt + m k=1 g jk (t, X(t))dW jk (t).
(2) W jk (t) is a standard Wiener process (Gaussian random variables with means zero, initial value zero with probability one and variance equal to t at time t), which we assume to be independent, where j = 1, 2, . . . , d and k = 1, 2, . . . , m. The drift f j (t, x) and diffusion g jk (t, x) are known functions which are assumed to be sufficiently regular (Lipschitz, bounded growth) for the existence and uniqueness of the solution (see, e.g., Øksendal 2003). The same solution process X(t) can be obtained from the Stratonovich SDEs denoted by: dX j (t) = f j (t, X(t))dt + m k=1 g jk (t, X(t)) • dW jk (t).
It is possible to switch between the two stochastic approaches by using a modified jth drift coefficient: (4) is called the drift correction formula. Note that the diffusion coefficients are the same in both representations. Such a Stratonovich SDE (3) is equivalent to an Itô SDE: g ik (t, X(t)) ∂ ∂X i g jk (t, X(t)) dt + m k=1 g jk (t, X(t))dW jk (t).
(5) Thus, we can write the Stratonovich SDE in the form of an Itô SDE. Whichever interpretation of an SDE is deemed appropriate in a particular situation, we can always switch to the corresponding SDE in the other interpretation when this is advantageous. In other words, (2) and (3) have the same solution. Notice that if the noise is additive (g independent of x), then the two representations are equivalent.
The current development of computing tools (software and hardware) has allowed solving many theoretical problems on SDEs that have become the object of practical research. This has also enabled many researchers in different domains to use these equations for modeling and analyzing practical problems. In response to these many uses, there already exist some open source packages that can perform stochastic calculus in the R software environment (R Core Team 2020); see the freely available packages sde (Iacus 2008(Iacus , 2016 and yuima (Brouste et al. 2014;Iacus and Yoshida 2018). They provide functions for simulation and inference for multidimensional SDEs, driven by Wiener process, fractional Brownian motion and Lévy processes, with or without jumps specified, as well as continuous time autoregressive moving average (CARMA), continuous time generalized autoregressive conditional heteroskedasticity (COGARCH) and point processes. The packages DiffusionRgqd, Diffusion-Rimp and DiffusionRjgqd (Pienaar and Varughese 2016a,b,c) constitute a handy collection of tools for performing inference and analysis using a data-imputation scheme and the method of lines on scalar and bivariate time-inhomogeneous diffusion and jump-diffusion processes with quadratic drift and diffusion. The QPot package (Moore, Stieha, Nolting, Cameron, and Abbott 2016) covers advanced techniques for quasi-potential analysis and visualizes the dynamics of two-dimensional systems of SDEs. The diffeqr package (Rackauckas and Nie 2017) is useful to solve SDE problems using the DifferentialEquations.jl package from the Julia programming language (Bezanson, Edelman, Karpinski, and Shah 2017).
In general, simulation allows to study and experiment with a given system where the complex interactions are known, to measure the effects of certain changes in the interactions on the behavior of the system, to experiment with new situations. It should be noted that if one looks for a faithful representation of the phenomena observed, one is quickly confronted with difficulties due to the calculations not being available in an explicit way. Monte Carlo (MC) techniques will allow us to approach these calculations numerically. There exists a well-developed literature on general Monte Carlo methods, we might mention Fishman (1996); Glasserman (2004); Liu (2004); Thomopoulos (2013); Graham and Talay (2013), and for computer vision with the R language see, e.g., Robert and Casella (2010). Methodological advances have led to much more computational demand in statistical computing, such as Markov chain Monte Carlo (MCMC) algorithms, bootstrapping, cross-validation, Monte Carlo simulation, etc. Many areas of statistical applications are experiencing rapid growth in the size of data sets (see, e.g., Mishra, Dehuri, Kim, and Wang 2016). Parallel computing is a common approach to resolve these problems. This is available in R since 2011 (version 2.14.0) through the base package parallel. This provides parallel facilities previously contained in packages multicore for parallel computation on shared memory (unix) platforms (multiprocessor and multi-core computers, Urbanek 2014) and snow for parallel computation on distributed memory platforms, with slight revisions (simple network of workstations, Tierney, Rossini, Li, and Sevcikova 2018). Several approaches to parallel computing are available in R, see Schmidberger, Morgan, Eddelbuettel, Yu, Tierney, and Mansmann (2009) ;Eugster, Knaus, Porzelius, Schmidberger, and Vicedo (2011); McCallum and Weston (2011), and an extensive and updated list of R packages is reported in the Comprehensive R Archive Network (CRAN) Task View on High-Performance and Parallel Computing with R at https://CRAN.R-project.org/view=HighPerformanceComputing (Eddelbuettel 2020).
The determination of response moments in SDEs is one of the fundamental problems in stochastic analysis and is much more complicated in the case of nonlinear equations. The main reason is the fact that the widest class of nonlinear SDEs subjected to additive or multiplicative random excitation do not have exact solutions in analytical forms. In many cases, therefore, it is necessary to adopt approximate solutions. In the moment equations method (MEM) approach (see, e.g., Bover 1978;Socha 2008;Alibrandi and Ricciardi 2012), the response statistical characterization is given by the moments. This method consists in writing a system of ODEs for the response statistical moments of any order, taking advantage of the Fokker-Planck equations or of the Itô differential rule. The MEM can be easily solved for linear systems. Unfortunately, for nonlinear systems in general, these equations constitute an infinite hierarchy. In fact the moment differential equations up to a given order contain moments of higher order and, consequently, this system is not solvable.
This has motivated us to take an interest in this particular field by introducing the R package Sim.DiffProc (simulation of diffusion processes, Guidoum and Boukhetala 2020). For the R package developed in the present paper we focus on two problems: Firstly, we develop a module for performing Monte Carlo methods whose essence is the use of repeated experiments to iteratively evaluate and approximate a quantity of interest (replicates of a statistic(s) to SDEs) at a fixed time t for scalar and multivariate Itô and Stratonovich SDEs, with the techniques for parallel computing on computer clusters, and multi-core systems using the parallel package. Secondly, we develop a module for symbolic and numerical computations by analytical determination of the dynamics of first and central second-order moments of an SDE (i.e., mean, variance and covariance of the dynamical variables). We also provide the built-in functions for random number generation based on the simulated trajectories at a fixed time t, and the evaluation of transition densities for SDEs in both interpretations, with nonparametric techniques. In addition, the simulation and estimation outputs can be graphically represented in 2D and 3D plots (like Figures 10). Our package is structured according to the S3 object-oriented system of classes and methods. Sim.DiffProc is currently the only package on CRAN, available at http://CRAN.R-project.org/package=Sim.DiffProc, that uses parallel processing on multiple processors and a cluster of computers in stochastic computing.
The paper is organized as follows: Section 2 we explain the implementation and describe the core functions of the package, and present the usage and options of the parallel Monte Carlo method, also we briefly review some theoretical concepts concerning the moment equa-tions method for an analytical approximation of first and central second-order moments. In Section 3 we give illustrations with simulations in different domains using scalar and multivariate Itô and Stratonovich SDEs. In Section 4 we give an overview of using the package to approximate the European options pricing for Black-Scholes (1D) and Heston models (2D) with an application using real data. The presentation of the results of a simulation study and moment equations with L A T E X tables/math expressions is demonstrated in Section 5. Concluding remarks discussing the package and its limitations, as well as future goals, are provided in Section 6.
Assumptions 1 and 3 require stationarity and/or ergodicity, and for assumption 2 neither stationarity nor ergodicity is necessary.
Let X δ be a time-discrete approximation of a continuous-time process X, with δ the maximum time increment of the step size.
1. An approximation method is strong order γ of convergence, if at fixed time T we have: 2. An approximation method is said to converge weakly of order β, if for any fixed time T and any 2(β + 1) continuous differentiable function g of polynomial growth, we have: In both cases δ 0 > 0 and C do not depend on δ.
Notice that all numerical methods implemented in this package are of strong order γ of convergence.
To obtain: • A formal display of SDEs, with all inputs. Creating L A T E X tables and mathematical expressions (reproducible research). • Numerical solutions of SDEs.
• Summary descriptive statistics for the solution process based on the simulated Mtrajectories at any fixed time t ∈ [t 0 , T ], i.e., mean, variance, median, mode, quantiles, minimum, maximum, coefficient of variation (relative variability), skewness, kurtosis, and the kth central moment.
• Kernel approximation of the transitional density (uni-, bi-and tri-variate) of SDEs, based on the simulation of M -trajectories at any fixed time t ∈ [t 0 , T ].
• Monte Carlo simulation-approximation of any interesting quantity at any fixed time t ∈ [t 0 , T ], with parallel computing.
• Symbolic and numerical computations by an analytical approximation of the first and central second-order moments (means, variances and covariances) of SDEs.
• Plot(s) routines for various classes of objects in the Sim.DiffProc package.
In Table 1 we give a summary of the key functions introduced in this package (where k ∈ {1, 2, 3}).

Parallel implementation of Monte Carlo method: 'MCM.sde' details
Monte Carlo methods are based on the analogy between probability and volume. They involve massive computation due to the replications in simulations (some problems are: computing column means, bootstrapping, etc.). These structures allow one to develop highly sophisticated Monte Carlo simulation based on the strong law of large numbers. In many cases, these perform extremely well in circumstances where other methods fail. Often, the purpose of the Monte Carlo simulation is to approximate an expectation: where X t is the solution of the SDE (2) (or (3)) at fixed time t, for some given function h by the expectation: where X ∆ t is the corresponding iterate of the numerical scheme with a constant step size ∆. In practice, (7) can only be approximated by the arithmetic mean: The approximation (8) depends on which M realizations of the numerical scheme are used to calculate it. To obtain a statistically reliable result, in the following, we apply the form with R batches with M realizations each. Then their mean is determined, which is the total sample mean of realizations X ∆ t (ω 1,1 ), . . . , The mean (9) is itself only an approximation of the desired value (6), sinceμ t,∆,M,R converges almost surely (i.e., for almost every generated sequence) to E (h (X t )) by the strong law of large numbers whenever E |h(X t )| < ∞. The confidence interval about (9) is then given by: whereσ 2 R is the unbiased sample variance of the batch means, and St R−1 (1 − α/2) is the inverse of Student's St cumulative distribution function with R − 1 degrees of freedom and confidence parameter α ∈]0, 1[, if R is small. For large R, we can use the normal distribution instead of the St-distribution (due to the central limit theorem). The accuracy of the simulation is therefore inversely proportional to the square root of the number of trials. This means that to double the accuracy the number of trials has to be quadrupled. It is clear that how good the approximation of (6) is will depend on the accuracy of the stochastic numerical scheme and step size ∆. A study of Komori, Saito, and Mitsui (1994) indicates that as ∆ decreases, the lack of independence in the samples from a random number generator (RNG) typically degrades the computation before rounding errors become significant. We will discuss this important feature in the following.
Generally, the Monte Carlo method on SDEs is not easy to implement, particularly in a multidimensional situation. Furthermore, the computational times are very significant. Hence computational performance can be increased by the use of parallel computing. Parallel computing is the simultaneous execution of the source code of one or more programs, to improve capacity and speed up the computations (see, e.g., Schmidberger and Mansmann 2010;Lim and Tjhi 2015). R offers a unified way of doing parallelization with the base package parallel, it can be easily run in parallel on several cores. There are several different approaches available to achieve parallelization and not all approaches are available for all platforms (see McCallum and Weston 2011). We will use two main methods to parallelize the Monte Carlo method: implicit parallelism (multiple cores on a single machine, i.e., multicore) with function mclapply(), and explicit parallelism (cluster with a multiple machine, i.e., snow) with function parLapply(). The conceptual differences between these two methods are illustrated in the following paragraphs.
The first method uses the forking 2 techniques (Urbanek 2014) on machines running with POSIX operating systems with multiple cores, and would simply call lapply() on Microsoft Windows operating systems, so that the code works but sees no speed gain. In contrast, on MacOS and Unix operating systems, the most straightforward way to enable parallel processing is by switching from using lapply() to mclapply(). This method will use all cores available to it. If you do not want to (either because you are on a shared system or you just want to save processing power for other purposes) you can set this to a value lower than the number of cores you have. Setting it to 1 disables parallel processing, and setting it higher than the number of available cores has no effect. We can easily check the number of cores we have access to with the command detectCores() available in package parallel.
The second method to parallel processing is more complicated, but works on any system including Windows systems. The computations in parallel can be extended to all cores on a single computer or to networked computers (see Appendix E), it uses the PVM, MPI, NWS standards as well as direct SOCKETS communication technologies for distributed memory computing, for more details on this technique see the original articles of Rossini, Tierney, and Li (2007) and Tierney, Rossini, and Li (2009). The general process we will follow is: i) Start a cluster with n nodes.
ii) Execute any preprocessing code necessary in each node (e.g., loading a package).
Results based on Monte Carlo simulations are not reproducible unless we set the random seed in their R code using function set.seed(). The parallel package contains support for multiple RNG streams based on the algorithms discussed in L'Ecuyer (1999) and L'Ecuyer, Simard, Chen, and Kelton (2002), with support for both multicore and snow clusters, and is the only RNG in R that supports multiple streams to generate random numbers in parallel environments. In addition there are also several packages on CRAN, see, e.g., rlecuyer (Sevcikova, Rossini, and L'Ecuyer 2019), rstream (Leydold 2020) and doRNG (Gaujoux 2020). In R we selected the RNG by calling RNGkind("L'Ecuyer-CMRG") and the clusterSetRNGStream() function (see help("RNGkind") and help("clusterSetRNGStream")). The idea of L'Ecuyer-CMRG RNG is that each separate worker generates random numbers independently from the others, and its period is around 2 191 , that will provide one stream of random numbers for about 264 workers. The seed of this generator can be easily advanced a given number of steps, making it very useful for true parallel RNG support.
Roughly, the best and easiest way to analyze the performance of an application is to measure the execution time. Consequently, an application can be compared with an improved version through the execution times: Speed-up ≈ execution time for a program without enhancements (sequential) execution time for a program using the enhancements (parallel) .
To summarize, the Monte Carlo simulation-approximation in this paper consists of: 1) Exhibiting a probabilistic representation of µ t (any interesting quantity at a fixed time t) of the type (6) such that the numerical solution of the stochastic system (2) is efficiently simulated.
2) Applying the strong law of large numbers in order to approximate (6).
3) To get the system status at fixed time T (i.e., the state of SDEs at this set time T ); we automatically use linear interpolation with the R function approxfun() available in base package stats (i.e., the state of the SDEs is linearly interpolated if T = n∆).
Remark 1. The whole idea might work in general when the simulation step size ∆ is very small (high frequency "assumption 2"). To use the numerical methods mentioned previously (Milstein scheme 1 and 2, Taylor scheme of strong order 1.5) the assumptions on the step size ∆ of simulation can be relaxed, but more regularity on the coefficients of drift and diffusion is needed (i.e., existence of partial derivatives up to order 2).
Parallel computing with multicore (shared memory) can be summarized into the following steps: 1) Detect the number of CPUs cores (ncpus) using detectCores().
3) Parallel computation of (9) with mclapply(), and leave mc.set.seed to TRUE (the mclapply() function can set the random number seeds for each worker for us, when the mc.set.seed argument is set to TRUE; we do not need to call clusterSetRNGStream()).
For parallel computing with snow (distributed memory) the steps are: 1) Start the cluster (cl); we use the function makePSOCKcluster().
2) Send a list of objects to all workers using the function clusterExport().
3) Load packages on all workers using the function clusterEvalQ().
4) Set parallel RNG on the cluster using the function clusterSetRNGStream().
The R function call consists of (see usage in the package help files): MCM.sde(model, statistic, R = 100, time, exact = NULL, names = NULL, level = 0.95, parallel = c("no", "multicore", "snow"), ncpus = getOption("ncpus", 1L), cl = NULL, ...) The input arguments for MCM.sde() are: model: an object of class 'snssde1d', 'snssde2d' and 'snssde3d'. statistic: a function that when applied to the model (SDEs) returns a vector containing the statistic(s) of interest. R: number of Monte Carlo replicates (R batches > 1), this will be a single positive integer. time: fixed time at which the quantity of interest statistic is approximated. exact: a named list giving the exact statistic(s), if it exists the bias calculation will be performed. names: names for the statistic(s) of interest. Default names = c("mu1", "mu2",...). level: confidence level of the required interval(s). parallel: the type of parallel operation to be used. "multicore" does not work on Microsoft Windows operating systems, but on Unix it uses parallel operations. Default parallel = "no". ncpus: an integer value specifying the number of cores to be used in the parallelized procedure. Default is 1 core of the machine. cl: an optional parallel cluster for use if parallel = "snow". Default cl = makePSOCKcluster(rep("localhost", ncpus)). ...: other named arguments for statistic which are passed unchanged each time it is called through the ... argument.
x: an object of class 'MCM.sde'. index: the index of the variable of interest within the output of class 'MCM.sde'. type: type of plots: "hist" computes a histogram, "qqplot" produces a normal QQ plot, "boxplot" produces a box-and-whisker plot, and "CI" plots confidence intervals. Default type = "all". ...: arguments to be passed to par() function to set graphical parameters. Tuckwell (1996, 2000) introduced an elegant and systematic analytical method which consists of replacing the system of stochastic differential equations with a system of deterministic equations representing the dynamics of the means, variances, and covariances of the components (2) (or (3)). This method requires much less computation since it involves only the solution of a system of coupled deterministic ordinary differential equations. Also, it leads directly to solutions for the moments of the state variables, such as their means and variances which are the fundamental properties that we seek for in a model. We report a slight improvement of their method by adding the case of Stratonovich interpretation; hence the approach will be adopted in both cases, Itô and Stratonovich SDEs.

Moment equations method: 'MEM.sde' details
Let us assume that the distribution function of X(t) is concentrated near the mean point X(t) = X 1 (t), X 2 (t), . . . , X d (t) (that is, P X(t) − X(t) < , for some (usually small) positive , is close to 1) and is symmetric around this point (see Rodriguez and Tuckwell 1996). For relatively small noise amplitudes, we take advantage of the fact that the third and higher-order odd central moments are close to zero and that the fourth and higherorder even moments are small relative to the second moment. Therefore, the exact means , respectively, which obey a system of (1/2)d(d + 3) ordinary differential equations.
The vector of d means m(t) for the various components X j , at time t is found to satisfy the following systems of ODEs: whereas the vector of d variances S(t) at time t is determined by: and the (1/2)d(d − 1) covariance C(t) at time t is determined by: with: where ν = 1 and 2 stand, respectively, for the case of Itô and Stratonovich SDEs.
For the general linear system the solutions of these equations give explicit results (as illustrated in Appendix B). Otherwise, the moment equations can be solved numerically. R offers several tools for solving a system of ODEs, many of these tools can be found in the CRAN Task View on Differential Equations at https://CRAN.R-project.org/view=DifferentialEquations (Soetaert and Petzoldt 2020). For an easy transition from a package such as deSolve (Soetaert, Cash, and Mazzia 2012), we include a function containing the symbolic computations of the system (11-13). For this we use the package Deriv (Clausen and Sokol 2020) and the D() function available in the base package stats, and a list of parameters and their values in function MEM.sde(), and we will use function ode() (general solver for ODEs) from package deSolve, to solve numerically the obtained system.
The R function call consists of (see usage in the package help files): MEM.sde(drift, diffusion, corr = NULL, type = c("ito", "str"), solve = FALSE, parms = NULL, init = NULL, time = NULL, ...) The input arguments for MEM.sde() are: drift: an R vector of expressions which contains the drift specification (1D, 2D and 3D). diffusion: an R vector of expressions which contains the diffusion specification (1D, 2D and 3D). corr: an R vector of expressions which contains the correlation structure between the Wiener processes (for any 2D and 3D). The default is corr = NULL. type: type of SDEs to be used; "ito" for Itô form and "str" for Stratonovich form. The default is type = "ito". solve: if solve = FALSE only the symbolic computation of system (11-13) will be made; and if solve = TRUE a numerical approximation of the obtained system will be performed. parms: parameters passed to drift and diffusion expressions. init: initial (state) values of system (11-13). time: time sequence (vector) for which output is sought; the first value of time must be the initial time. ...: arguments to be passed to function ode() available in package deSolve, if solve = TRUE. The returned object can be summarized using: summary(object, at, ...) which provides the summary results at a fixed time, defined by the argument at.
object: an object of the class 'MEM.sde'. at: fixed time between minimum and maximum of time sequence (vector). The default is at = T.
For dynamic plot(s) of outputs, we will use the functions available in the deSolve package, e.g., matplot.0D().

Model specification and functionalities
In the following, the functionality of the package is illustrated by examples. We present the main functions and their applications in univariate and multivariate cases (Itô and Stratonovich SDEs) to demonstrate their flexibility, the potential for parallelization, and customization for stochastic analysis. For each example, we start by setting the seed and modifying the R options using R> set.seed(1234, kind = "L Ecuyer-CMRG") R> options(scipen = 10, digits = 5) so that readers can replicate the results.
For t ≥ 0, the models (14) and (15) have log-normal transition density p θ (t, X t |X t 0 = x 0 ); i.e., the density of the distribution of X t given X t 0 = x 0 , with the mean and variance of its logarithm transform (i.e., the log-mean and log-variance) for each SDE type, are given by: and σ 2 mod1 = θ 2 t µ mod2 = log(x 0 ) + 1 2 θ 2 t and σ 2 mod2 = θ 2 t with mean and variance: Before talking about the approximation of transition density p θ (t, X t |X t 0 = x 0 ) for each type of model, we will illustrate in a one-dimensional linear case that the mean and variance for (14) and (15), coincide exactly with results obtained by MEM. We obtain the symbolic computation of a system of two ordinary differential equations for the mean and variance for each model using the MEM.sde() function (procedure outlined in Section 2.4). We can do this by simply running the following commands.
Moment equations of (14): Moment equations of (15): The exact solution of two symbolic systems of ODEs obtained by the function MEM.sde(), gives exactly the results previously known for mean and variance of both models, with initial conditions m(0) = x 0 and S(0) = 0. We obtain the numerical results by assigning the value TRUE to the argument solve, the initial value by c(m = x0, S = 0) and a vector containing the time points at which the solution is determined by seq(0, 1, by = 0.001). In this case, we have assigned the outputs for each model to objects called mem.mod1 and mem.mod2 for further use in the R session.  (14) and (15) calculated using function MEM.sde(), which exactly corresponds to the known values.
the simulated trajectories of SDEs (i.e., output of function snssde1d()). The dsde1d() function (based on function density() available in base package stats) can be used to show the kernel approximation for p θ (t, X t |X t 0 = x 0 ). This consists of giving the initial condition of the SDE, i.e., the starting value x0, and the starting time t0. These two elements can be defined by assigning values to the arguments of function snssde1d(). The argument at = 1 in function dsde1d() represents the fixed time t where the transition density is to be evaluated. For each type of SDE, the kernel estimation of transition density at ending time T = 1 are reported, see the left panel of Figure 2. Since the density approximation is evaluated over time, it can best be visualized with arranged multiple density plots in a staggered fashion, also known as ridgeline plots. Motivated by this, we can use the ggridges package (Wilke 2020), which provides a convenient way of visualizing changes in distributions over time, and we are reshaping the outputs with function melt() available in the reshape (or reshape2) package (Wickham 2007); see the right panel of Figure 2 (for the R code see Appendix A).
Log mean and log variance of (14) and (15) (14) and (15) and true transitional density at ending time T = 1, with initial condition t 0 = 0 and X 0 = 1. Left: X mod1 t=1 |X 0 = 1 and right: X mod2 R> plot(dens.mod1, dens = function(x) dlnorm(x, meanlog = mu.mod1, + sdlog = sigma.mod1), hist = TRUE, xlim = c(0, 8)) R> plot(dens.mod2, dens = function(x) dlnorm(x, meanlog = mu.mod2, + sdlog = sigma.mod2), hist = TRUE, xlim = c(0, 8)) Now we give an illustration of the approximation of a quantity of interest (e.g., mean and variance for an SDE at a fixed time t) using simulated data by the parallelized Monte Carlo method (procedure outlined in Section 2.3), and we compare the results obtained with the exact values. The exact values of mean and variance for each model are declared in R using: We define the function of a quantity of interest by stat.fun1d in the current R session: We apply the Monte Carlo method at a fixed time T = 1, using function MCM.sde(), with parallel = "snow" and ncpus = 4 on Windows on a single machine (on a Unix platform, we can use parallel = "multicore" and for computer clusters see Appendix E).  Figure 4: QQ plots of the quantiles of the estimator of mean and variance for (14) and (15) from parallel Monte Carlo versus the theoretical quantile values from a normal distribution.
R> set.seed(1234, kind = "L Ecuyer-CMRG") R> mcm.mod1 <-MCM.sde(model = mod1, statistic = stat.fun1d, R = 100,  (1), S = V.mod2(1)), + parallel = "snow", ncpus = 4) R> mcm.mod2 The QQ plot (see Figure 4) indicates that the circles are all very close to the red line (theoretical quantile values from a normal distribution), close enough to say these estimators come from a normal distribution. There is a little random wriggle around the red line, which does not disqualify these estimators from being normal.  R> plot(mcm.mod1, index = 1, type = "qqplot") R> plot(mcm.mod1, index = 2, type = "qqplot") R> plot(mcm.mod2, index = 1, type = "qqplot") R> plot(mcm.mod2, index = 2, type = "qqplot") In Table 2, we give for each type an approximation of the mean and variance obtained by function MCM.sde() with various batches at a fixed time T = 1. Through the results obtained, we find that the parallelized Monte Carlo method is very robust and offers excellent performance. Moreover, the central limit theorem guarantees convergence and asymptotic normality of estimators. Due to the parallel random number generation support, this avoids the problem of lack of independence in the samples from multiple RNG. Hence a replication number of few independent realizations is sufficient to obtain a suitable approximation for the first and central second-order moments to provide information on the shape of the distribution.
The standard error (Std.Error) and root mean square error (RMSE) are chosen as criteria for comparison between the Itô and Stratonovich interpretations. In the case of the Stratonovich form, note that the two criteria are higher than those found in the Itô form (see Figure 5). This is explained by the fact that stochastic calculus with the Stratonovich form increases the variability of a statistic, due to the drift correction formula (4).
The tie between bootstrap and Monte Carlo simulation of a statistic is obvious, both are based on repetitive sampling and then a direct examination of the results. A big difference between the methods, however, is that bootstrapping uses the original, initial sample as the population from which to re-sample, whereas Monte Carlo simulation is based on setting up a data generating process (with known values of the parameters). This strategy could be used to approximate some quantities, like in the bootstrap, but also to theoretically investigate some general characteristic of an estimator which is hard to derive analytically. Table 3 shows speed-up and computing time for mean and variance of (14) and (15) approximated by function MCM.sde() at a fixed time T = 1, with different number of processors. This timing test was done on a 2.20 GHz Intel Core i5-5200U processor running Unix platform (Ubuntu 17.04 operating system), and 8 GB 1600 MHz DDR3 of memory. We notice that the execution time with multicore (shared memory), which is suitable for Unix systems, is faster than snow (distributed memory).
The nonlinear mean reversion model Aït-Sahalia (1996) proposed a nonlinear model including polynomial terms for modeling interest rates. In the following, we study the simplest version of this model which satisfies the nonlinear SDE: Some natural restrictions have to be imposed on the parameter values θ = (β 0 , β 1 , β 2 , β 3 , σ, ρ) to have a meaningful specification of the model; for more details we refer the reader to the Appendix in Aït-Sahalia (1996). In general, there are no exact distributional results. This model has been estimated empirically using various techniques by many authors including Aït-Sahalia (1996), Conley, Hansen, Luttmer, and Scheinkman (1997) and Gallant and Tauchen (1998) just to mention a few references. Later Aït-Sahalia (1999 proposed an approx-imation of transition densities via the Hermite polynomial expansion 3 . We shall report the results of the simulation analysis under assumption 2 (high frequency: T = n∆ fixed, ∆ → 0 as n → ∞) only for the case: β 0 > 0, β 1 < 0, β 2 ≤ 0, β 3 > 0, σ > 0 and ρ > 1.
From the Markovian property of the diffusion, it is also possible to define the stochastic transient of the model (16) by the time evolution of the probability density function P (x; t), which satisfies a Fokker-Planck equation also called Kolmogorov forward equation (e.g., see Risken 1996): An analytic solution of the partial differential equation (17) is unlikely due to the nonlinearity of the model (16). By using moment equations method (11-12), we obtain a set of two equations of moments (i.e., m(t) mean and S(t) variance) in closed form through function MEM.sde(). Furthermore, if the user is interested in exporting the equations displayed by printing 'MEM.sde', to be used in a L A T E X document, see Section 5.2. In R: The stochastic transient with initial conditions (m(0) = x 0 , S(0) = 0) is numerically solved from the system (18). We do this by just running the following command: R> mem.mod <-MEM.sde(drift = f, diffusion = g, type = "ito", solve = TRUE, + init = c(m = 5, S = 0), time = seq(0, 5, by = 0.001)) R> summary(mem.mod, at = 5) Approximation of moment at time 5 | m(5) = 2 | S(5) = 0.06557 In Figure 6 we show the results obtained by the Monte Carlo simulation of model (16) and the system of moment equations (18), where there is an excellent agreement between the two approaches over time t ∈ [0, 5]. The following code has been used to produce Figure 6.
With dcSim(): R> parm <-c(beta0, beta1, beta2, beta3, sigma) R> drift <-function(t, x, parm) parm[1] + parm[2] * x + parm[3] * x^2 + + parm[4] * x^-1 R> diffu <-function(t, x, parm) parm[5] * x^1.5 R> x <-seq(1, 3.5, length = 100) R> dens <-NULL R> for (to in x) { + dens <-c(dens, sde::dcSim(x0 = 5, x = to, t = 5, d = drift, + s = diffu, theta = parm, N = 100)) + } Indeed, the approximate solutions of the transition density with functions dcSim() and dsde1d() are nearly identical, as we see in the left panel of Figure 7. Now, we consider in particular the cumulative distribution function (CDF) F X (x) of a random variable X with the function defined at any point x as: There are two well-known estimators for the CDF. The first estimator is the empirical cumulative distribution function (ECDF), and the second is the kernel estimator of the distribution function (KCDF). In fact, we can quickly estimate (19) for model (16) in a specific time T at any point x by both approaches using the formulation (9):  where 1 {·} is the indicator function and H(x) = x −∞ K(t)dt, with K(·) the kernel function and h the bandwidth parameter. R has a function ecdf() to compute the ECDF of an individual data set. The function returns a 'function' object to compute the value of F n (x) at any point x. Let us write an R function ecdf.fun to compute (20a) as follows: Fn <-stats::ecdf(d) + return(c(Fn(seq(1.5, 3.75, by = 0.25)))) + } To calculate the second estimator (20b), we need to use function kcdf() available in package ks (Duong 2007). Furthermore, we calculate the bandwidth h using the 1-and 2-stage plugin selectors of Polansky and Baker (2000), with function hpi.kcde() (to go further into the theory see the recent article of Duong 2016). In R we define the function below: Fhat <-ks::kcde(d, h = ks::hpi.kcde(d, nstage = stage)) + return(predict(Fhat, x = seq(1.5, 3.75, by = 0.25))) + } Finally, the use of MCM.sde() is quite simple, i.e., as in the previous examples. The numerical results of both methods are reported in Table 4, and shown in the right panel of Figure 7.

Ornstein-Uhlenbeck process and its integral
The Ornstein-Uhlenbeck (OU) process has a long history in physics. Introduced in essence by Langevin (Lemons and Gythiel 1997) in his famous 1908 paper on Brownian motion, the process received a more thorough mathematical examination several decades later by Uhlenbeck and Ornstein (1930) and Wang and Uhlenbeck (1945). The OU process is understood here to be the univariate continuous Markov process X t . In mathematical terms, the equation is written as an Itô equation: In this equation, µ and σ are positive constants called, respectively, the relaxation time and the diffusion constant. θ ∈ R is the asymptotic mean of X t . The time integral of the OU process X t (or indeed of any process X t ) is defined to be the process Y t that satisfies: Y t is not itself a Markov process; however, X t and Y t together comprise a bivariate continuous Markov process (see Gillespie 1996). We wish to find the solutions X t and Y t to the coupled time-evolution (21) and (22): We first take note of the following well-known result, which can be derived from the system (23). For any t > 0 the OU process X t and its integral Y t will have a bivariate normal distribution with means and variances given by (where the approximations of the first and central second-order moments with MEM coincide exactly with the known values, see Appendix B for demonstration): and the covariance is: We now have a complete and exact solution to the problem of the time evolution of X t and Y t . Next, we will use this information to compare with parallelized Monte Carlo simulation. For (23), we simulate 5000-trajectories using function snssde2d(), with µ = 1, θ = 2, σ = 0.5, t ∈ [0, 20] and initial value (x 0 , y 0 ) = (−5, 0). The simulation output of (23) is assigned to an object called OUI. We can plot the system via plot(OUI), and use summary(OUI, at = 20) to obtain the descriptive statistics of R> set.seed(1234, kind = "L Ecuyer-CMRG") R> mu <-1 R> theta <-2 R> sigma <-0.5 R> x0 <--5 R> y0 <-0 R> f <-expression(1/mu * (theta -x), x) R> g <-expression(sqrt(sigma), 0) R> OUI <-snssde2d(drift = f, diffusion = g, M = 5000, T = 20, + x0 = c(x0, y0)) First, we declare the exact values of the first and central second-order moments of (23), for comparison with the estimated results.
Next, we apply the parallelized Monte Carlo method for model (23) (20))) R> mcm.oui Itô Sde 2D: | dX(t) = 1/mu * (theta -X(t)) * dt + sqrt(sigma) * dW1(t) | dY(t) = X(t) * dt + 0 * dW2(t) | t in [0,20]  We notice that a number of 100 independent replications of parallelized Monte Carlo is usually enough to get a good approximation at a fixed time T = 20. The resulting plots of simulation over time t ∈ [0, 20] are given in Figure 8, where there is good agreement between results of Monte Carlo and exact values of the first and central second-order moments of system (23).
We have several packages for creating a 3D plot in R, which can be very useful for visualizing the features of a more complex density. For static plotting, we use the scatterplot3d package (Ligges and Mächler 2003) and interactive plotting is provided by rgl (Adler and Murdoch 2020). The arguments pdf = "Joint" and at = 20 produce a kernel approximation of the bivariate density function (based on the function kde2d() available in the package MASS; Venables and Ripley 2002) of the couple (X t , Y t ) at any fixed time using function dsde2d(), with several forms of display ("contour", "image", "persp", "rgl"), e.g., see the right panel of Figure 9. The argument pdf = "Marginal" produces an approximation of a marginal density function of the process. We use the ggridges package to show the evolution of the marginal density of X t and Y t at t = {0, 2, 4, . . . , 20}, see the left panel of Figure 9. We may generate contour plots of the transition density at fixed time points t = {0, 1, . . . , 20} with display = "contour", the animation plots are made using the animation package (Xie 2013); see the left panel of Figure 10.

Simulated Expectation
Transition Density at t = 0 x y R> plot(denJ, display = "rgl", col2d = "lightblue") The bivariate cumulative distribution function (BCDF) F (X,Y ) (x, y) of two random variables X and Y for a couple of fixed level (x, y), are given by: We can easily approximate (24) at any fixed time T by the empirical BCDF. As in a onedimensional case, using formulation (9) for the EBCDF as: and we compare it with the exact bivariate cumulative distribution of system (23), i.e., where Φ corresponds to the bivariate normal cumulative distribution function, with mean vector µ and variance-covariance matrix Σ at fixed time T . In R, we begin by declaring the exact values of (26) at ending time T = 20 in coupled points (x ∈ (1, 2, 3), y ∈ (28, 33, 38)).

Bivariate nonlinear dynamic system: The general stochastic Haken-Zwanzig equations
The stochastic version 6 of a bivariate Haken-Zwanzig model (Friedrich and Haken 1988;Leung and Lai 1993) is used to demonstrate the slaving principle in Synergetics (Haken 2004) and it is also related to the more complex ABCDE model (see, e.g., Friedrich and Haken 1992). We propose the general stochastic model as Stratonovich equations: where rate parameters (µ, α) are non-negative, (σ 1 , σ 2 ) ≥ 0 describes the amplitude of the fluctuating quantity, k = 0, 1 and 2 respectively for additive noise, linear and nonlinear multiplicative noises. Obtaining the analytical solution of (27) is not possible, numerical methods are required. A complete description of the solution implies assignments of values to the two rate parameters (µ, α) and to the two initial conditions (x 0 , y 0 ). This model is stable and has three fixed points: the first point is e u = (0, 0) describing the extinction of the dynamic process, it is an unstable source; and the second (two) points are e s = ± √ µα, µ representing the coexistence of two-modes (x, y), and they are stable since the eigenvalues of stability are given by: where i denotes the imaginary number √ −1. The flow field, fixed points e u and e s , together with the nullclines are shown in Figure 11, which is made using package phaseR (Grayling 2014); for the R code see Appendix C.1.
The stochastic transient of system (27) is now described by a probability density function P (x, y; t), which is found to satisfy a Fokker-Planck equation (forward), and is expressed in the Stratonovich representation by: Should we be able to solve P (x, y; t) analytically, then statistical moments can be calculated explicitly to provide relevant information (for example, the peak location and width of P (x, y; t) are determined by the first two orders moments). However, since an explicit solution of (28) is unlikely; one can quickly obtain the final results of the moment equations for each type of noise with function MEM.sde(), yet numerical treatment is required to solve these coupled equations with initial conditions: m 1 (0) = x 0 , m 2 (0) = y 0 , S 1 (0) = S 2 (0) = C 12 (0) = 0, and for parameters µ = 1 and α = 3. The system (27) allows us to treat and study the effect of the three types of noise: if k = 0 we have an additive noise, which signifies that the two representation Itô and Stratonovich are equivalent; and if k = 1, 2 we have linear and nonlinear multiplicative noise, respectively. We were assigning appropriate values as σ 1 = 0.1 and σ 2 = 0.05; for a concise presentation, we shall report results only for σ 1 > σ 2 . In the case k = 0, we choose the initial state x 0 = y 0 = 0, and in the other two cases x 0 = y 0 = 1.
The right panel of Figure 12 shows the effect of noise on system dynamics. We clearly see the acceleration of the transition to stable points e s = ± √ µα, µ of (27), which also appears in Figure 13 where the transition of the limit cycle to a fixed point e s was predicted in numerical studies of the steady state probability. We found that in cases k = 1 and 2 influence of noise becomes state dependent since the structures of g(x, y) function (diffusion) which are no longer constant (linear and nonlinear function). Usually, additive noise causes smaller effects than multiplicative noises, and also could slow down or speed up the transient process of a dynamic system. From the simulation of system (27), we find that linear multiplicative noise (i.e., the case k = 1) causes significant effects (that is clearly shown in the right panel of Figure 12), due to the fact that the response of this system to noise is nonlinear. For a qualitative comparison of the effects of noise, we evaluate the relative fluctuation R j (t) for the three cases, which is given by: m j (t) %, for j = 1, 2.
As previously, we give an illustration of the approximation of R j (t) at ending time T = 100 with function MCM.sde(), for case k = 1 only (same procedure for the cases k = 0 and 2).

R> -digamma(1)
[1] 0.57722 In Figure 16 we compare the results obtained from MEM with those of MCM. We notice that the means and variances for the two methods are practically indistinguishable, but there is a bad agreement between the covariance curves (C MEM 12 (t) = C MCM 12 (t)). This is due to the assumptions for the validity of the MEM which are not satisfied, implying that agreement with simulation is poor, i.e., the distribution function of system (32) is not concentrated near the mean point.
We conclude this example by estimating the parameters for the system (32). Since we have explicit theoretical moments (33a-33d), we want to estimate the vector of the parameters (µ 1 , µ 2 , σ 1 , σ 2 ) by the method of moments; this consists of using the empirical moments. For this, the empirical moments x t and y t , which we respectively equal to the theoretical moments (33a) and (33b) as follows: and the same for second-order theoretical moments (33d), that we equal to the sample variance: The estimators for µ 1 , µ 2 , σ 1 , σ 2 denoted byμ 1 ,μ 2 ,σ 1 ,σ 2 are described as the solution to the following system: from (33c), we directly get the estimatorν: where z t is the numerical solution of the ODE Z t (since z t = z t (ω j ) and VAR(z t ) = 0, ∀j = 1 . . . M for all t fixed).

Transformation of a scalar SDE into three dimensions
Next we consider a scalar SDE driven by three independent 8 Wiener processes (W 1,t , W 2,t , W 3,t ): To simulate the numerical solution of (35) with package Sim.DiffProc, we make a transformation to a system of three equations as follow: 8 To simulate stochastic models (3D) driven by correlated Wiener processes, see Appendix D.2.
Moreover, if we want to test the normality of the distribution, i.e., if a sample is normally distributed or not, we use the Shapiro and Wilk (1965) normality test. The null hypothesis of this test follows a normal distribution, so if p-value < 0.1, one conludes that the sample does not follow a normal distribution (Royston 1995). This then allows using some statistical tests if the null hypothesis is retained. In R we use function shapiro.test() available in base package stats (we can also use the Kolmogorov-Smirnov test with function ks.test(), it allows testing if a sample follows any given distribution), and we define the function of statistics of interest (i.e., approximation of W The MC results return a non-significant p-value = 0.42369 ± 0.02906, so we cannot reject the null hypothesis, i.e., the model (35) is normally distributed at a fixed time T = 5.
Now we want to approximate the parameters of the distribution of (35) at fixed time T = 5; we will use the fitdistr() function available in package MASS (maximum likelihood fitting of univariate distributions, see help("fitdistr")). As previously, we define the function of statistics of interest, i.e., the estimate of parameters of the normal distribution by estimate.fun (in fact of any distribution, see help("Distributions")), and we use easily function MCM.sde().

Application to European options pricing
In order to illustrate how the R package Sim.DiffProc can be used in practical situations, we give an overview of using the package to approximate the prices of Black-Scholes (1D) and Heston models (2D) European call and put options, with application to real data sets.

Black-Scholes European options pricing
The price of a European option C(s, x) at time s is the discounted value at the risk-free rate of interest r with maturity T , that is (e.g., see Iacus 2011, Chapter 6): where X s,x T is the price of the underlying at maturity T . The payoff function f (x) for European call and put options with strike price K, is given by: We can easily approximate (38) using the formulation (9).
We will use the function GBSOption() available in package fOptions (Wuertz, Setz, and Chalabi 2017). This function is handy to calculate the exact formulas for the call and put options of the generalized Black-Scholes model. So, for example, we calculate the price of a contract with the following parameters: the price of the underlying S 0 = 100, the strike price K = 150, the time to maturity T = 1, the interest rate r = 0.1, the cost of carrying term µ = 0.75, and the volatility σ = 0.25.
R> set.seed(1234, kind = "L Ecuyer-CMRG") R> fx <-expression(mu * x) R> gx <-expression(sigma * x) R> BS <-snssde1d(drift = fx, diffusion = gx, x0 = S0, M = 5000) Secondly, it is quite easy to implement the quantity of interest (40) in R. Let us write a function EuroOpt.fun to compute the call and put option price of the Black-Scholes model (or indeed of any process X t ).
R> EuroOpt.fun <-function(data, i, r, K, T, s) { + d <-data[i, ] + call <-mean(pmax(0, d -K)) * exp(-r * (T -s)) + put <-mean(pmax(0, K -d)) * exp(-r * (T -s)) + return(c(call, put)) + } Now we compare the Monte Carlo approximation of the price and the exact price of a European call/put options already obtained by function GBSOption(). As we see, the parallelized Monte Carlo method gives a reasonable approximation of the price of a European call/put option to maturity T = 1. The resulting plots of simulation in time t ∈ [0, 1] are given in Figure 19, where there is good agreement between the results of Monte Carlo and exact values. So, when the formula of the price is not available in explicit form, for example, the underlying process is not a Black-Scholes model (for example model (16)), the MCM.sde() function is useful to evaluate the European options pricing of any model.

European options pricing under stochastic volatility: Heston model
Under the Heston (1993) model, the stock price X t follows a Black-Scholes type stochastic process, and its stochastic variance Y t follows a Cox, Ingersoll, and Ross (1985) process.
Hence, the Heston model is described by the bivariate system of SDEs: where the two Wiener processes are correlated, that is E (dB 1,t dB 2,t ) = ρdt. In the second equation, volatility is mean-reverting to its long-run level θ with speed ν; σ is the volatility of volatility. If the condition 2νθ > σ 2 holds, then the drift is sufficiently large for the variance process Y t to be guaranteed positive and not reach zero (known as the Feller condition). The correlation is an important part of the model (41). If the correlation is negative, volatility increases as the underlying asset falls and volatility decreases as the underlying rises. In particular, for the Monte Carlo simulation, it is more convenient to first generate independent Wiener processes; we can easily translate these correlated Wiener processes into uncorrelated ones by Cholesky decomposition which gives the decomposition of (B i,t ) i=1,2 onto the orthogonal basis (W i,t ) i=1,2 : Note that E (dB 1,t dB 2,t ) = ρdt. Moreover, since (W i,t ) i=1,2 are two independent Wiener processes: E (dW 1,t dW 2,t ) = 0 and E dW 2 We can rewrite the Heston model (41) in the orthogonal basis as: Now let us use the functions callHestoncf() and putCallParity() available in R package NMOF (Gilli, Maringer, and Schumann 2019). These functions are useful for calculating the price of a European call and put with the characteristic function of the log price in the Heston model. In R we calculate the value of the call and put under the following parameters: the price of the underlying x 0 = 100, initial variance y 0 = 0.3, the strike price K = 120, the time to maturity T = 1, risk-free rate µ = 0.15, dividend yield q = 0, long run variance θ = 0.25, speed of mean reversion ν = 2, volatility of volatility σ = 0.2 and finally the correlation ρ = 0.75.

European options pricing approximation to real data sets
Now, if we look at the market price of, e.g., a call or put option at a given fixed time, we can compare the price predicted from real data by the exact Black and Scholes formula with (40)   therefore interested in the period ending in February 2019, and we get the 1000 closing levels of S&P500 11 prior to that date. For that purpose, we use function get.hist.quote() from the tseries package (Trapletti and Hornik 2019).

R>
[1] 342.83 125.24 As in the previous example we simulate M paths of the Black-Scholes model, with the two estimated parameters of the modelr andσ, and X 0 = S 0 , ∆t = 1/252. R> set.seed(1234, kind = "L Ecuyer-CMRG") R> fx <-expression(r.hat * x) R> gx <-expression(sigma.hat * x) R> BS <-snssde1d(drift = fx, diffusion = gx, x0 = S0, Dt = 1/252, M = 5000) We apply the parallelized Monte Carlo method on the previous function EuroOpt.fun. Finally, we find the price of the call and put options of 1000 closing levels of S&P500. We obtain an excellent agreement between the exact values and the result obtained by the parallelized Monte Carlo method (see the left panel of Figure 22). Thus, in practice when the explicit formulas for calculating the price of the options do not exist, it is always possible to rely on the formula (40).
Indeed, τ = 2.608 then corresponds to 12 October 2017. Moreover, this is shown in the left panel of Figure 22. The right panel of Figure 22 shows that the strike price K is an increasing function of the time τ at which the price of call equals the price of put. We summarize the numerical results obtained in Table 7.

Miscellaneous tools for the Sim.DiffProc environment
A simple utility writing available in function TEX.sde() produces the related L A T E X codes (tables and math expressions) using the Sim.DiffProc package, which can be copied and pasted in a scientific article. Note that there are already several tools available for the R environment for generating L A T E X code. Many of these tools can be found in the CRAN Task View on Reproducible Research at https://CRAN.R-project.org/view=ReproducibleResearch (Blischak and Hill 2020). In the following we present the new tools for constructing L A T E X tables and mathematical expressions with the Sim.DiffProc package. The R function call consists of: TEX.sde(object, ...) The input parameters for TEX.sde() are: object: an object of class 'MCM.sde' or 'MEM.sde'; or an R vector of expressions of SDEs, i.e., drift and diffusion coefficients. ...: if object of class 'MCM.sde', arguments to be passed to kable() function available in package knitr (Xie 2015).

L A T E X table for objects of class 'MCM.sde'
The Monte Carlo results of the 'MCM.sde' class can be presented in terms of L A T E X tables.
Our goal is to make the TEX.sde() method modular and flexible. We thus leave tasks, such as formatting of the table, possible to be modified by the user. For this, we are using function kable() (for more details see help("kable")) available in package knitr (Xie 2015), which converts an output of 'MCM.sde' class into the corresponding L A T E X code for a table. In Section 3.1 we have assigned the outputs of the Monte Carlo method to an object called mcm.mod. In R, we create the L A T E X table for this object using the following code: R> TEX.sde(object = mcm.mod, booktabs = TRUE, align = "r", caption = + paste("\\LaTeX~table for output results of class \\class{MCM.sde} ", + "generated by \\code{TEX.sde()} method.") %%% LaTeX   Table 8 shows the results of applying our TEX.sde() method to an object of class 'MCM.sde'.
If the results of the simulation become complicated enough, creating L A T E X tables from R significantly reduces the work, especially if the simulation study needs to be repeated due to bug fixes, enhancements, or modified implementation.

L A T E X mathematics for objects of class 'MEM.sde'
Since many reports are produced using R and L A T E X, it is very useful to be able to automatically convert mathematical expressions from one language to another. In the following, our goal is to automatically convert an R expression for an object of class 'MEM.sde' into its appropriate L A T E X representation. We will use the R codes available in the book of Wickham (2015, p. 320-328), with some modifications, improvements and adaptation with the Sim.DiffProc package. To give some concrete examples, we want to automatically generate the L A T E X code appropriate for the moment equations obtained from the model (23) using method TEX.sde().

Discussion
In this work we have presented the most used methods in stochastic dynamics systems, analyzing in detail the moment equation method and Monte Carlo simulation (version parallel). The MCS approach is very general and works in almost all circumstances especially when the simulation mesh is very small (under assumption 2 "high frequency"). For high-dimensional functionals, it is sometimes the only method of obtaining a solution for real engineering problems, and it is a powerful tool to approximate an interesting quantity in an empirical manner, especially when we do not have an analytic form of a solution. Besides being more efficient than analytical-based approaches, it has the advantages that the tools of deterministic analysis can be fully exploited. As a drawback, the computational costs required are often very high. The computational efforts increase with the dimension and the complexity of the nonlinear structural models. The cost of this generality is the large sample sizes (dependence between simulated samples) required to obtain reasonably accurate estimates. The parallel RNG support is perhaps the most interesting and important feature of parallel computing, which allows us to avoid the problem of lack of independence in the samples from multiple RNG.
The moment equation method is advantageous for analyzing the effects of noise in stochastic systems. These moments implicitly include information on relative fluctuations R(t). Sometimes, it provides explicit results (generally in the linear case); otherwise, we can quickly solve these equations by standard numerical methods. This avoidsgreedy calculations of time associated with the need to average over many batches of MCS. The moments obtained compared favorably with those obtained by MCS when the noise was small. Stochastic analysis always requires some approximation treatment. The moment equation method provides here a concise picture of a stochastic transient. Generally in the linear case we find explicit results, but in the nonlinear case it is very difficult, as we have seen in Section 3.3. It furthermore enables an analytical study of the transient characteristics so that noise effects could be explicitly and systematically compared. Results derived from Section 3.2 demonstrate that there is no effect of the type of stochastic model (Itô or Stratonovich) on the relative fluctuation. We have confirmed these results analytically in Section 3.1, where we easily find R mod1 (t) = R mod2 (t). Since the response of a dynamic system to noise is nonlinear, noisy systems could reveal more complex and exciting phenomena than deterministic systems. Further study of the effects of noise on nonlinear systems is needed.
We have demonstrated how the Sim.DiffProc package may be used to simulate, estimate and analyze the dynamics of scalar and multivariate systems of Itô and Stratonovich SDEs driven by uncorrelated and correlated Wiener processes. Moreover, we highlighted the performance of symbolic and numerical computation of the package, through several examples and applications. Two central components of the package are the approximation of a quantity of interest for a stochastic model with parallelized MC, and an analytical approach for the determination of the first and central second-order moments, which are highly useful for studying stochastic dynamics. The R language is also a powerful high-level programming language that is well suited for expressing complex mathematical computations, as stochastic analysis of SDEs. There can be a tension between language features that support interactive use, and features that support the development of reliable, high performance code. L A T E X mathematics are complex. Fortunately, they are well documented. That said, they have a fairly simple structure, which could be used as a quick starting point for the more complex editing of a scientific article using R and L A T E X. We hope that the package presented here and the updated survey on the subject might be of help for practitioners and researchers in the field who might want to implement new methods and ideas using R package Sim.DiffProc. One of the most significant limitations of the latest version of the package (version 4.8), is that jump-type stochastic models are not developed in this version, making modeling limited to a class of SDEs. Future efforts will work toward extending the package to more complex systems, driven by fractional Brownian motion and Lévy processes, with jumps specified.

E. Running 'MCM.sde' quickly on many machines
The snow technique is essentially a very coarse grained remote strategy call. It relies upon delivery copies of code and information to remote procedures and then returning results. It is ill-suited for little assignments but very well suited for a reasonable number of moderate to large tasks. This is the technique utilized by R's parallel package. This strategy may appear to be less productive and less advanced than shared memory strategies, but relying on the transmission of objects, it is in principle very easy to extend the technique of a single machine to several machines. This is the strategy consiered in the following. Do the following: Collect a list of machine addresses that you can use. This is the hard part; it depends on your operating system, and something you should get help with if you have not already tried it. This assumes that you already understand how to configure secure shell (SSH) and for a Windows box use PuTTY (see the vignette of package parallel), for more details see the website https://www.stat.uiowa.edu/~luke/R/cluster/cluster.html. In this case, we use IPV4 addresses and hostnames. In our case the list is: • Master machine: 2.20 GHz Intel Core i5-5200U processor, 8 GB 1600 MHz DDR3 of memory. IPV4: "192.168.1.105" and user: "acguidoum".
Note that we do not collect passwords because we assume that we have configured the appropriate authorized_keys and key pairs in the .ssh configurations of all these machines. We call the machine we use to issue the "master" global calculation. It is essential to try all these addresses with ssh in a terminal shell before trying them with R. Also the machine address you choose for the "master" must be an address that the machines can use to return to the main machine. Now, with the system things resolved, the R part is as follows. Start your cluster with: R> machineAddresses <-list(list(host = "192.168.1.105", user = "acguidoum", + ncore = 4), list(host = "192.168.1.108", user = "TM161", ncore = 4))