Conditional Optimization of Algorithms for Estimating Distributions of Solutions to Stochastic Differential Equations

: This article discusses an alternative method for estimating marginal probability densities of the solution to stochastic differential equations (SDEs). Two algorithms for calculating the numerical– statistical projection estimate for distributions of solutions to SDEs using Legendre polynomials are proposed. The root-mean-square error of this estimate is studied as a function of the projection expansion length, while the step of a numerical method for solving SDE and the sample size for expansion coefficients are fixed. The proposed technique is successfully verified on three one-dimensional SDEs that have stationary solutions with given one-dimensional distributions and exponential correlation functions. A comparative analysis of the proposed method for calculating the numerical–statistical projection estimate and the method for constructing the histogram is carried out.


Introduction
Many mathematical models of dynamical systems in various fields of science (radio engineering, statistical mechanics, automatic control, chemistry, medicine, reliability theory, etc.) can be described by stochastic differential equations (SDEs).Moreover, some boundary value problems in mathematical physics can be reduced to Cauchy problems for SDEs.This approach is based on both a probabilistic representation of the solution to the boundary value problem in the form of a functional of the solution to the corresponding SDE [1] and a numerical method for statistical modeling of this solution [2][3][4][5][6][7].
For analyzing, filtering, predicting, or smoothing random processes in various dynamical systems described by SDEs, the probability density is often estimated [8,9] because it provides maximum information about the random process.
Mathematical models related to phase transitions, e.g., processes of particle association into clusters, can be described by SDEs [10].There are many examples of such processes: formation of aerosols in the atmosphere, condensation in high-speed gas flows, polymerization, crystallization, deposition of metal vapors, etc.Such problems require estimating distributions of solutions to SDEs.
Numerical methods for solving SDEs [2][3][4][5][6][7] make it possible to find approximate solutions and also to estimate their probability densities.As such an estimate, the histogram is usually considered [11][12][13][14][15], which is a statistical analogue for the probability density of the SDE solution.Note that histograms are used not only when solving SDEs but also when analyzing various data [16][17][18].However, this article studies estimating distributions in the context of SDEs.
Consider the following SDE: dX(t) = G(t, X(t))dt + σ(t, X(t))dW(t), t ∈ [0, T]; X(0) = X 0 , where X(t) is an n X -dimensional vector random process; G is a given vector-valued function, σ is a given matrix-valued function and W(t) is an n w -dimensional standard Wiener process with independent components.The distribution of an initial vector X 0 is given.The initial vector X 0 and the Wiener process W(t) are mutually independent.Since SDE (1) can be understood as the Itô SDE or the Stratonovich SDE, we will use different notations for them.The Formula (1) is accepted for the Itô SDE.The Stratonovich SDE has the following formula: where g is a vector-valued function A transition from the Itô SDE to the Stratonovich SDE and vice versa is always possible.Note that if σ is a constant matrix, then the Itô SDE and the corresponding Stratonovich SDE coincide.
This article presents an alternative method for estimating marginal probability densities of solutions to SDEs.The statistical projection estimates in the Monte Carlo method were first proposed by N.N.Chentsov [19].He developed a general technique for optimizing such estimates, which requires improvement in specific problems [20,21].The randomized projection estimates based on Legendre polynomials [22] for marginal probability densities of solutions to SDEs are constructed in this article.The problem of optimizing the root-mean-square error is considered.The problem of the optimal, or consistent, choice is studied for the following parameters of statistical algorithms: the sample size; the step size in a numerical method for solving SDEs; and the projection expansion length.The proposed method has been successfully tested on several examples.A comparative analysis of the proposed method for calculating the numerical-statistical projection estimate and the method for constructing the histogram are carried out.

Numerical Methods for Solving the Cauchy Problem for SDEs and Their Conditional Optimization
To solve problems whose mathematical models are specified by SDEs, various numerical methods can be used.Let us give the necessary definitions [2].

Definition 1.
A numerical method has the mean-square order of convergence q (on the Cauchy problem solution) if where X n is the numerical solution at the time t n = nh; h is a given constant step size, and E denotes the mean or mathematical expectation.We will also say that the method has the qth order of mean-square convergence.
Note that the order of mean-square convergence is often considered to be the square root of q defined above [4,6].Definition 2. A numerical method weakly converges with order q (on the Cauchy problem solution) if for any function s(X) that belongs to a class F with its partial derivatives up to order 2q + 2, where F is the class of functions ν(X), for which there are constants C > 0 and κ ≥ 0 such that the inequality |ν(X)| ≤ C(1+|X| κ ) holds for all X ∈ R n X .We will also say that the method has the qth order of weak convergence.
Let us give some examples of numerical methods for solving the Cauchy problem for SDEs.To solve the Itô SDE (1) numerically, we can use the Euler-Maruyama method [3], which has the first order of mean-square convergence: where {ξ n } is a sequence of mutually Gaussian random vectors with independent compo- nents having standard normal distribution, and ξ n and X n are independent.
To solve the Stratonovich SDE (2) numerically, we can use the generalized Rosenbrocktype method [2]: where I is the identity matrix of size n X × n X , and σ •j is the jth column of the matrix σ.
The numerical method (4) is asymptotically unbiased with any step size h > 0. This means that if one applies it with the given step size h to the one-dimensional linear SDE where λ and σ are the constants, and λ > 0; then the distribution of the numerical solution x n converges as n → ∞ to the normal distribution with zero mean and the variance σ 2 /2λ.And if the numerical method (4) is applied with the given step size h to the multidimensional linear SDE where A and σ are constant matrices, and all eigenvalues of the matrix A have negative real parts; then the distribution of the numerical solution X n converges as n → ∞ to the normal distribution with zero mean and the covariance matrix D satisfying the Lyapunov equation The asymptotic unbiasedness property is a generalization of the A-stability property of numerical methods for solving ordinary differential equations.The asymptotically unbiased method (4) is recommended for solving stiff SDEs.Note that the Euler-Maruyama method (3) does not have the asymptotic unbiasedness property with any step size h > 0.
In the general case, the numerical method (4) has the first order of mean-square and weak convergence.However, for SDEs driven by the one-dimensional Wiener process and for SDEs with a constant matrix σ, the numerical method (4) has the second-order mean-square and weak convergence.
Numerical methods for solving SDEs mentioned above can be used to estimate the required functionals, e.g., moments or cumulants of the SDE solution.Let us introduce notations for the following functionals: where X(t) is the piecewise linear random process obtained from X n .To estimate the functional of the SDE solution, we should simulate N paths of the random process X(t) using some numerical method.Further, J(h) is estimated by the arithmetical mean of the obtained sample values for f (h); i.e., Then, the estimation error J N (h) is determined by the relation where C 1 , C 2 are positive constant independent of N and h, and D denotes the variance.
To reduce the computational cost of the statistical algorithm for calculating the mean of some functional, the problem of the optimal, or consistent, choice of algorithm parameters, e.g., the sample size N and the step size h, arises.The following theorem has been proved in [2]: Theorem 1.Let J N (h) be an estimate of some functional of the SDE solution at time t n ∈ [0, T] obtained by the numerical method for solving SDEs with the step size h.Then, the minimum computational cost for calculating the functional is achieved for parameters where γ is the required calculation accuracy.
The corresponding computational cost of the statistical algorithm is where C is a positive constant independent of N and h.

The Projection Expansion of Marginal Probability Densities of the SDE Solution
Due to the computational complexity of orthogonal expansions with adapted weights, randomized projection estimates of solutions to integral equations using Legendre polynomials have been constructed in [21].Next, we apply a similar technique to construct the numerical-statistical projection estimate of the probability density of the SDE solution using Legendre polynomials.
Legendre polynomials have the following advantages compared to other complete orthonormal systems of functions.Firstly, the expansion coefficients of the probability density are conveniently expressed in terms of the moments of the SDE solution.Secondly, there is a simple estimate for the error of the probability density approximation by Legendre polynomials if it is continuously differentiable r times (in this case, r defines the rate of convergence).
Let p(z) be the probability density of the random variable z = x(t n ), a ≤ z ≤ b, which is a component of the exact solution to SDE at time t n ; i.e., x is an arbitrary component of the random vector X.Next, we find the projection estimate of the marginal probability density p(z) of the SDE solution at time t n .
Consider the projection estimate (with expansion length m) of the function p ∈ L 2 [a, b] by the following expansion using orthonormal Legendre polynomials ψ i (z) [22] defined on a given interval [a,b]: where and the expantion coefficients a i are calculated as inner products, i.e., and they satisfy the recurrence relation For orthonormal Legendre polynomials on the interval [a,b], the following formula can be used: and the recurrence relation is satisfied for them as follows: This relation is obtained by substituting y = (2z − (b + a))/(b − a) into ( 7): The last equation together with (7) leads to (9).The randomization of the projection estimate ( 5) is obtained by calculating the linear functional by the Monte Carlo method using the numerical solution to SDE.
The randomized estimate of the marginal probability density of the SDE solution at time t n is constructed as follows: where x n is the numerical solution to SDE at time t n for the jth sample path, and N is the sample size.
For kth order unconditional moments of the SDE solution at time we have their statistical estimates: The limitation of the proposed method is only in the condition However, as a rule, the marginal probability density of the SDE solution is a continuous function Further, we propose two algorithms based on the above description for the randomized estimate of the marginal probability density of the SDE solution.The first algorithm is based on the Formula (8), which gives the explicit representation for Legendre polynomials.Therefore, this algorithm uses statistical estimates of moments of the numerical solution to calculate the coefficients in the Expansion (5).
Algorithm 1 for the numerical-statistical projection estimation of the probability density of the SDE solution on the interval [a,b] at time t n is based on m + 1 Legendre polynomials: (1) Simulate N paths of the SDE solution using a numerical method on the interval [0, t n ] and find the numerical solution at time t n : n , j = 1, . . ., N; (2) Calculate statistical estimates of moments of the numerical solution: (3) Calculate randomized estimates of coefficients a i in the Expansion (5): (4) Calculate the randomized projection estimate of the probability density: at nodes z s , where s = 0, 1, . . ., n z and n z is the number of nodes of the interval [a,b].The recurrence Relation (9) for Legendre polynomials avoids computer arithmetic errors for coefficients in the explicit Formula (8) for large m.Using the recurrence Relation ( 9) is recommended under condition m > 20 since in this case, computer arithmetic errors arise for the explicit Formula (8).However, in this case, it is more convenient to estimate the mean of Legendre polynomials to calculate the coefficients in the Expansion (5).Therefore, the second algorithm is proposed.Algorithm 2 for the numerical-statistical projection estimation of the probability density of the SDE solution on the interval [a,b] at time t n is based on m + 1 Legendre polynomials: (1) Simulate N paths of the SDE solution using a numerical method on the interval [0, t n ] and find the numerical solution at time t n : n , j = 1, . . ., N; (2) Calculate randomized estimates of coefficients a i in the Expansion (5) using the mean of Legendre polynomials and the recurrence Relation (9) for them: (3) Calculate the randomized projection estimate of the probability density: Otherwise, the additional estimation error arises, which increases with increasing the projection expansion length m.In particular, Step 2 in Algorithm 1 can be rewritten as follows: Step 2 in Algorithm 2 is modified in a similar way: In fact, the interval [a,b] should be chosen so that the probability of violating the condition

The Error Analysis
Let us estimate the deviation of the numerical-statistical projection estimate p m (z) from the marginal probability density p(z) of a random variable First, we use Jensen's inequality: Second, we apply Parseval's equality: where Finally, we use the variance definition and the following identity: where C 1 , C 2 , C 3 are positive constants independent of N, h, and m; r is the order of continuously differentiability of the marginal probability density p(z).In the Formula (10), we take into account both the variance of coefficient estimates a i [23] and the error of the estimates for functionals of the SDE solution obtained by numerical method that has the qth order of weak convergence.We also use the corollary of Theorem 4.10 from [22] about the error of the function approximation by Legendre polynomials: if a function p(z) is continuously differentiable r times on the interval [a,b], then Further, we consider the problem of the optimal, or consistent, choice of the following parameters of the statistical algorithm for calculating the projection estimate: the sample size N, the step size h in a numerical method, and the projection expansion length m.Theorem 2. Let p m (z) be the numerical-statistical projection estimate of the marginal density of the SDE solution at time t n ∈ [0, T] obtained by Algorithms 1 or 2.Then, the minimum computational cost of calculating the projection estimate is achieved for parameters where γ is the required calculation accuracy in the norm of space L 2 [a, b], i.e., Proof.Since the estimation requires functionals of the SDE solution, we should take into account Theorem 1. Setting h ≃ N −1/2q in the Formula (10), we have where C 4 is a positive constant independent of N, h, and m.To choose optimal parameters N opt and m opt , it is sufficient to write the equality for the resulting errors and to obtain the required orders from the relation We find that the following relations are conditionally optimal: Then, we have optimal orders for N, h, and m: If we choose parameters in this way, then where C is a positive constant independent of N, h, and m.Theorem 2 has been proved.□ Note that the one-dimensional probability density p(t,X) of the solution to the Itô SDE (1) satisfies the forward Kolmogorov equation, which is also called the Fokker-Planck-Kolmogorov equation: where b jl (t, X) are elements of the diffusion matrix B(t, X) = σ(t, X) σ T (t, X).
For the existence of a classical solution to the forward Kolmogorov equation, at least the condition p(t, • ) ∈ C 2 [a, b] should be satisfied for every t.Therefore, for a consistent choice of parameters in the Equation (10), the condition r = 2 is required.Then, Choosing parameters in this way, we have B(p, p m ) ≤ CN −3/8 , where C is a positive constant independent of N, h, and m.
The first two terms in the Formula (10) are quite small.The main term of the estimation error is given by the Formula (11).It decreases with increasing the order of continuous differentiability r of the marginal probability density p(z); consequently, fewer Legendre polynomials are required to achieve a given estimation error.In fact, the third term in Equation ( 10) is virtually zero for infinitely continuously differentiable functions, and the remaining terms linearly depend on m.Therefore, the projection estimate error can increase with increasing the projection expansion length m.Now we can compare the numerical-statistical projection estimate of the marginal probability density of the SDE solution constructed above with the histogram that is usually considered for estimating the probability densities.
The deviation of the histogram π * (z) from the marginal probability density p(z) of a random variable x(t n ) in the norm of space L 2 [a, b] is estimated in [2].The monograph [2] also considers the problem of optimizing the computational cost of constructing the histogram, i.e., the optimal, or consistent, choice of the following parameters: the sample size N; the step size h; and the histogram step h g , where h g = (b − a)/n g , and n g is the number of histogram nodes.The following theorem has been proved in [2]: Theorem 3. Let π * (z) be the histogram (constructed from a sample of size N), which is a statistical analogue for the marginal probability density p(z) of the SDE solution at timet n ∈ [0, T] obtained by the numerical method for solving SDE with the step size h.Then, the minimum computational cost of constructing the histogram is achieved for the parameters where γ is the required calculation accuracy in the norm of spaceL 2 [a, b], i.e., Choosing parameters in this way, we have B(p, π * ) ≤ CN −1/3 , where C is a positive constant independent of the sample size N, the step size h, and the histogram step h g .This is greater than the error of the numerical-statistical projection estimate proposed above.
Remark 2. In addition to the histogram, the kernel density estimate is also used as an estimate of the probability density of the SDE solution.In this context, the problem of the optimal, or consistent, choice is considered in [8] for the following parameters of the statistical algorithm: the sample size; the step size in a numerical method for solving SDEs; and the smoothing parameter.

Numerical Experiments
Real random impacts on dynamical systems are often differentiable processes, and there is a need to constructively define them.The problem of constructing the stationary Gaussian or non-Gaussian random processes when solving stochastic problems by statistical modeling is also of independent interest [24].Constructing a class of such processes-Markov diffusion processes-is possible by applying SDEs.
The monograph [2] describes a method for representing the random process with both a given one-dimensional probability density p(x) and an exponential correlation function as the solution to SDE where λ is a positive constant, and m 1 is arbitrary.If x(t) is the stationary solution to SDE (12), where then x(t) has the one-dimensional probability density p(x) and the correlation function If the probability density of the initial value x 0 is equal to p(x) in SDE (12), then its solution is stationary over the entire interval [0,T].If x 0 is given arbitrarily, then x(t) will be stationary after a certain transition time.By varying the parameter λ, we can change the correlation time for the process x(t): τ c ≈ 1/λ.Note that the stationary random process constructed in this way is a Markov diffusion process with continuous paths.
A numerical study of the proposed method for estimating marginal probability densities was carried out using SDE (12).This choice is due to the fact that the probability density is known over the entire interval [0,T].This allows one to estimate the approximation error for it and to confirm the theoretical result given by Theorem 2.
In the following examples, the probability densities of stationary random processes are infinitely continuously differentiable functions.As a numerical method for solving Stratonovich SDEs, the asymptotically unbiased method (4) is used.
For the numerical implementation of the proposed method, independent pairs of standard normally distributed random variables ξ 1 , ξ 2 were drawn by the Box-Muller approach: where α 1 , α 2 are independent uniformly distributed random variables on the interval (0,1).For modeling uniformly distributed random variables the pseudo-random number generator RND128 (with modulus 2 128 and multiplier 5 100109 ) [25] is used.
For all the examples, we have analytical expressions for both the stationary probability density and the moments of the solution.Therefore, the error of the numerical-statistical projection estimate can be calculated in the norm of space L 2 [a, b] as follows: All the tables for different sample sizes N show the projection estimate errors ε m , depending on the projection expansion length m and the histogram errors ε g for a consistent choice of the sample size N and the number of histogram nodes n g (see Theorem 3).

Example 1. The stationary solution to the Stratonovich
is the stationary random process with uniform distribution: , ν = 1, 2, . . .; The following parameters are used: α = 1; a = 0; b = 1.The initial realizations of the initial value x(0) were drawn by the pseudo-random number generator RND128 mentioned above.Table 1 gives the numerical results for t = 0 for different sample sizes N and projection expansion lengths m to compare root-mean-square errors of both the projection estimate and the histogram when the sample is exactly simulated.These results show that the projection estimate error decreases with increasing the sample size as N -1/2 , and it increases with increasing the projection expansion length m.The histogram error decreases with increasing sample size as N -1/3 in the case of a consistent choice of parameters, i.e., for n g,opt = N 1/3 .For all given projection expansion lengths, the projection estimate error is less than the histogram error.Table 2 shows the numerical results for t = T = 1 in order to compare the errors of both the projection estimate and the histogram, when the distribution is simulated by the numerical method (4) with step sizes h = 0.01 and h = 0.001.In this case, the projection estimate error also increases with increasing the projection expansion length This is similar to the case t = 0; therefore, Table 2 contains the projection estimation error only for m = 4 and m = 15.Theorem 2 implies that the statistical error dominates for the projection estimate with parameters N = 10 6 and h = 0.01; hence, reducing the step size h is justified, and the consistent parameters N = 10 6 and h = 0.001 lead to decreasing the error.
According to Theorem 3, parameters N = 10 6 and h = 0.01 for the histogram are consistent; consequently, reducing the step size h with the same sample size N does not lead to significant changes in the histogram error.
Table 2 certainly illustrates that the projection estimate error is smaller than the histogram error.
The graphs of the uniform probability density of the stationary solution to SDE, the projection estimate of the probability density, and the histogram corresponding to the numerical solution to SDE are shown in Figure 1 is the stationary random process with beta distribution: , where Γ denotes the gamma function.
The following parameters are used:  Example 2. The stationary solution to the Stratonovich SDE is the stationary random process with beta distribution: where Γ denotes the gamma function.
It is well known that the probability density of beta distribution for integer parameters γ and δ is the probability density of a γth order statistic for γ + δ -1 independent random values with uniform distribution on the interval (0,1) [26].Therefore, the simulation of the initial value x(0) in this example for γ = δ = 2 was carried out as the second order statistics for three independent random values with uniform distribution.Table 3 contains the numerical results for t = T = 1.It presents the root-mean-square errors of both the numerical-statistical projection estimate that depends on the projection It is well known that the probability density of beta distribution for integer parameters γ and δ is the probability density of a γth order statistic for γ + δ -1 independent random values with uniform distribution on the interval (0,1) [26].Therefore, the simulation of the initial value x(0) in this example for γ = δ = 2 was carried out as the second order statistics for three independent random values with uniform distribution.Table 3 contains the numerical results for t = T = 1.It presents the root-mean-square errors of both the numerical-statistical projection estimate that depends on the projection expansion length m and the histogram that depends on the number of histogram nodes n g .Table 3 is similar to Table 2 for uniform distribution.The results given in Table 3 demonstrate that with a consistent choice of parameters, the projection estimate error is smaller than the histogram error.For m = 4, the projection estimate very accurately approximates the probability density.
Figure 2 shows graphs of the probability density (beta distribution) of the stationary solution to SDE, the projection estimate of the probability density, and the histogram corresponding to the numerical solution to SDE (N = 10 4 ; h = 0.01; m = 4; n g = 20).Similar to Figure 1, Figure 2 demonstrates that the projection estimate more accurately approximates the exact probability density for the same sample size N.The results given in Table 3 demonstrate that with a consistent choice of parameters, the projection estimate error is smaller than the histogram error.For m = 4, the projection estimate very accurately approximates the probability density.
Figure 2 shows graphs of the probability density (beta distribution) of the stationary solution to SDE, the projection estimate of the probability density, and the histogram corresponding to the numerical solution to SDE (N = 10 4 ; h = 0.01; m = 4; ng = 20).Similar to Figure 1, Figure 2 demonstrates that the projection estimate more accurately approximates the exact probability density for the same sample size N. is the stationary random process with Gaussian distribution: )!!, 1, 2,... is the stationary random process with Gaussian distribution: The following parameters are used: α = 1; m 1 = 0; σ = 1.
The initial value x(0) was simulated exactly using the Box-Muller approach.The probability density of the solution was estimated on the interval [-4,4].Table 4 shows the numerical results.
In this example, the exact solution is the Gaussian process, and it takes values from (−∞, +∞).Therefore, we should use proposed algorithms, taking into account Remark 1, since Legendre polynomials are defined on the interval [ −4,4].So, statistical estimates for both the moments of the solution and the mean of Legendre polynomials should be calculated for solutions that belong to the interval [−4,4] only.
Table 4 shows the numerical results.In this example, the accuracy of the projection estimate depends on whether m is even or odd.In fact, the transition from even m to odd m + 1 does not increase the accuracy.This is apparently a specific property of Gaussian distribution, since corresponding moments with odd orders are zero.Due to the above property, the accuracy of the projection estimate increases with transition from m to m + 2. Figure 3 gives graphs of the Gaussian probability density of the stationary solution to SDE and the projection estimate of the probability density; it also contains the histogram of the numerical solution to SDE (N = 10 4 ; h = 0.01; m = 9; n g = 20).Figure 3 demonstrates that the projection estimate more accurately approximates the probability density of the solution for the same sample size N. Table 4 shows the numerical results.In this example, the accuracy of the projection estimate depends on whether m is even or odd.In fact, the transition from even m to odd m + 1 does not increase the accuracy.This is apparently a specific property of Gaussian distribution, since corresponding moments with odd orders are zero.Due to the above property, the accuracy of the projection estimate increases with transition from m to m + 2.

Conclusions
The main result of this article is two algorithms for calculating the numerical-statistical projection estimate for distributions of solutions to SDEs using Legendre polynomials.The root-mean-square error of this estimate is studied as a function of the projection expansion length, while the step of the numerical method and the sample size for expansion coefficients are fixed.The proposed technique is successfully verified on three one-dimensional SDEs that have stationary solutions with given one-dimensional probability densities and exponential correlation functions.A comparison of proposed numerical-statistical projection estimates and histograms demonstrates that the projection estimate more accurately approximates the probability density of the SDE solution.In addition, the projection estimate is presented in analytical form, it is a polynomial, i.e., smooth function.
Higher approximation accuracy is important for analyzing, filtering, predicting, or smoothing random processes in various dynamical systems described by SDEs.The smoothness of the projection estimate is more preferable when finding the solution to the boundary value problem in the form of a functional of the solution to the corresponding SDE.

Remark 1 .
at nodes z s , where s = 0, 1, . . ., n z , and n z is the number of nodes of the interval [a,b].The computational costs of both algorithms are approximately the same, but Algorithm 2 is more computationally reliable for large m.If the SDE solution is defined on the interval larger than the interval [a,b], on which Legendre polynomials are defined, then for Step 2 in Algorithms 1 and 2, the estimates of moments of the numerical solution and the estimates of the mean of Legendre polynomials should be calculated only for solutions that belong to the interval [a,b].

Mathematics 2024 ,
12, x FOR PEER REVIEW 13 of 17εg 0.03900 (ng = 20) 0.02268 (ng = 40) 0.00938 (ng = 100) The graphs of the uniform probability density of the stationary solution to SDE, the projection estimate of the probability density, and the histogram corresponding to the numerical solution to SDE are shown in Figure 1 (N = 10 4 ; h = 0.01; m = 4; ng = 20).

Table 4 .Figure 3
Figure3gives graphs of the Gaussian probability density of the stationary solution to SDE and the projection estimate of the probability density; it also contains the histogram of the numerical solution to SDE (N = 10 4 ; h = 0.01; m = 9; ng = 20).Figure3demonstrates that the projection estimate more accurately approximates the probability density of the solution for the same sample size N.
Figure3gives graphs of the Gaussian probability density of the stationary solution to SDE and the projection estimate of the probability density; it also contains the histogram of the numerical solution to SDE (N = 10 4 ; h = 0.01; m = 9; ng = 20).Figure3demonstrates that the projection estimate more accurately approximates the probability density of the solution for the same sample size N.

Table 1 .
Errors of estimating the probability density at time t = 0 (Example 1).

Table 2 .
Errors of estimating the probability density at time t = T = 1 (Example 1).

Table 3 .
Errors of estimating the probability density at time t = T = 1 (Example 2).

Table 4 .
Errors of estimating the probability density at time t = T = 1 (Example 3).