Spectral Representations of Iterated Stochastic Integrals and Their Application for Modeling Nonlinear Stochastic Dynamics

: Spectral representations of iterated Itô and Stratonovich stochastic integrals of arbitrary multiplicity, including integrals from Taylor–Itô and Taylor–Stratonovich expansions, are obtained by the spectral method. They are required for the implementation of numerical methods for solving Itô and Stratonovich stochastic differential equations with high orders of mean-square and strong convergence. The purpose of such numerical methods is the modeling of nonlinear stochastic dynamics in many ﬁelds. This paper contains necessary theoretical results, as well as the results of numerical experiments.


Introduction
Differential equations are used to describe dynamical systems, whose states change in continuous time.The classical theory of differential equations assumes that their solutions are sufficiently smooth deterministic functions.However, to take into account random actions or disturbances, it is necessary to consider differential equations with random functions (random processes).If these random functions are sufficiently regular, then the classical theory can be applied to study solutions.For irregular random functions -such as white noise -the classical theory is not applicable, and, in its development, the theory of stochastic differential equations appeared [1][2][3].Stochastic differential equations are used to describe stochastic dynamical systems in many fields [4][5][6][7].
In terms of applications, the methods for finding solutions to specific types of equations are important, but this is a very difficult problem.As a rule, when discussing the analytical solutions of Itô and Stratonovich stochastic differential equations, linear equations (or non-linear equations that are reduced to linear ones by a change in variables) are considered.Most of the examples of analytical solutions for stochastic differential equations are constructed via the inverse problem, i.e., according to the analytical expression of some random process.It is through this that the corresponding stochastic differential equation is derived.Since analytical solutions of stochastic differential equations can only be found in some particular cases, there is a need for methods that allow one to approximate solutions.
There are numerical methods based on the expansions of the solution of stochastic differential equations: Taylor-Itô and Taylor-Stratonovich expansions [17,21], as well as the unified Taylor-Itô and Taylor-Stratonovich expansions [20].These expansions are analogs of the Taylor series expansion of the solution of the ordinary differential equation.They contain iterated Itô and Stratonovich stochastic integrals over Wiener processes.The difference schemes for corresponding numerical methods contain the first terms of Taylor-Itô or Taylor-Stratonovich expansions, and the maximum multiplicity of iterated stochastic integrals in these terms determines the order of convergence for a specific numerical method.Although these expansions contain iterated stochastic integrals with the simplest weights (or, after transformations, with monomial weights), these integrals are generally not expressed in an elementary way in terms of random variables, for which there exist effective simulation algorithms.
For the approximate simulations of iterated Itô and Stratonovich stochastic integrals, we can use numerical integration; however, it is more preferable to apply the orthogonal expansion.This approach is described in [22], where the author approximates the iterated stochastic integrals of a second multiplicity using trigonometric functions.The numerical method using this approximation is called the Milstein method [13,17,20].A similar approximation for the iterated stochastic integrals of a third multiplicity appeared in [17].The approximation for integrals of a second multiplicity was discussed in [23,24], and their distribution was studied in [25].In general, without taking into account iterated stochastic integrals in a difference scheme, it is impossible to obtain a high order of mean-square or strong convergence of the numerical method [26].
Significant progress in the approximation of iterated Itô and Stratonovich stochastic integrals of arbitrary multiplicity was achieved in [20,27,28].The proposed approach involves the transition to corresponding multiple stochastic integrals with their subsequent representation by multiple series with random coefficients.In this case, the expansion of the deterministic functions of several variables into generalized multiple Fourier series is used.Most of the attention is paid to Legendre polynomials and trigonometric functions as the basis to represent the iterated stochastic integrals of arbitrary multiplicity [29], but Walsh and Haar functions are also used to represent the iterated stochastic integrals of a second multiplicity [20].
This paper proposes exact and approximate representations of iterated Itô and Stratonovich stochastic integrals using the spectral method [30].These representations are based on matrix expressions that use several spectral characteristics of functions, operators, and functionals.For iterated Itô stochastic integrals, it avoids relations with the indicators [20] or the Wick product [31] of random variables that correspond to Wiener processes.In this context, iterated Itô and Stratonovich stochastic integrals are considered earlier with the simplest weight functions [32][33][34].Here, we propose a general case of weight functions and consider mixed-type iterated stochastic integrals of arbitrary multiplicity.In addition to theoretical results, this paper contains the results of numerical experiments.

Itô Stochastic Differential Equations and Taylor-Itô Expansion
We consider the nonlinear Itô stochastic differential equation as follows [35]: where X(•) is the n-dimensional random process, t ∈ T = [t 0 , T], W(•) is the s-dimensional Wiener process with independent components, X 0 is the n-dimensional random vector, f (•) : T × R n → R n is the given vector function, and σ(•) : T × R n → R n×s is the given matrix function, R = (−∞, +∞).The random vector X 0 and the Wiener process W(•) are independent.
As a rule [35], it suffices to assume that f (•) and σ(•) are measurable functions, and they satisfy both the linear growth condition and the Lipschitz condition with respect to the second argument.However, we additionally assumed that f (•) and σ(•) are sufficiently smooth functions.Thus, all of the derivatives of these functions appearing in the Taylor-Itô expansion exist [17]: where ϑ t 0 , h > 0, and ϑ + h T. Expansion (2) uses the following notations (φ(•) : where σ * j (•) is the jth column of the matrix function σ(•).
The notation used for the iterated Itô stochastic integral (3) is based on the fact that it can be represented as a particular case of the multiple Itô stochastic integral [20].In fact, I J W(j 1 ...j k ) [ϑ,ϑ+h] is the linear operator that associates the function with the multiple Itô stochastic integral for this function.The number of function arguments should be the same as the integral multiplicity.In this case, the linear operator I J W(j 1 ...j k ) [ϑ,ϑ+h] acts on the function where and k(t) ≡ 1 for k = 1.However, below we consider a more general iterated Itô stochastic integral where ψ l (•) : T → R is the bounded weight function, l = 1, . . ., k.This notation means that the linear operator I J W(j 1 ...j k ) [ϑ,ϑ+h] acts on the function where In particular, if at least one of the values j 1 , . . ., j k is equal to zero (the corresponding integral is called the mixed-type iterated stochastic integral), then the integral (3) can be reduced to a linear combination of integrals (5) with a smaller multiplicity and the additional condition [20] The function k ψ (•) under condition ( 7) is denoted by k n 1 ...n k (•).For example, where the number of function arguments coincides with the integral multiplicity, and it suffices to use the notation J By substituting Relations (8) between iterated Itô stochastic integrals into Formula (2), the unified Taylor-Itô expansion is obtained [20].The difference is that Expansion (2) contains iterated stochastic integrals for functions (4) only, but some of them are mixedtype iterated stochastic integrals.In the unified expansion, there are no mixed-type integrals due to the transition to functions (6) under condition (7).These expansions are equivalent, but the latter expansion contains fewer types of iterated stochastic integrals.
Based on the Taylor-Itô expansion, a great deal of numerical methods for solving stochastic differential equations have been proposed [13,17,18,20].To describe some of them, we should introduce a discretization of the interval [t 0 , T] with a given constant step size h (a variable step size can also be considered).A discrete-time approximation for the random process X(•) is defined by a numerical method such as the following: The type and the order of convergence are very important for numerical methods.The numerical method has the order of mean-square convergence p if max i∈{1,...,N} Ch p , and the numerical method has the order of strong convergence p if max i∈{1,...,N} where E means the expectation, C > 0 is a constant independent of step size h, and h → 0. Note also that the second inequality follows from the first one.This can be proved using the Jensen inequality [36].
To construct a numerical method with the order of mean-square or strong convergence p, it suffices to leave only the first terms on the right-hand side of Expansion (2).For example, we can write the explicit difference schemes for the stochastic differential equation (1) as follows: p = 0.5 (Euler-Maruyama method) : p = 1.0 (Milstein-type method) : p = 1.5 (Strong 1.5 order method) : p = 2.0 (Strong 2.0 order method) :

Stratonovich Stochastic Differential Equations and Taylor-Stratonovich Expansion
We also consider the nonlinear Stratonovich stochastic differential equation [35]: for which the vector function a(•) : T × R n → R n is added to the notations for Equation (1), where the symbol • is to distinguish Itô and Stratonovich differential equations and stochastic integrals.If the condition holds, then Equations ( 1) and ( 13) are equivalent, i.e., they define the same random process X(•).
We assume that the function a(•) is sufficiently smooth, i.e., all of the derivatives of this function that are required in the Taylor-Stratonovich expansion exist [17]: L 0 a ϑ, X(ϑ) where ϑ t 0 , h > 0, and ϑ + h T.
The additional notation for Expansion ( 14) is where S J is the linear operator that associates the function with the multiple Stratonovich stochastic integral [20] for this function, and k(•) is defined by Formula (4).
For numerical methods, it suffices to restrict ourselves to functions k ψ (•) under condition (7), i.e., functions k n 1 ...n k (•), since the integral (15) can be reduced to a linear combination of integrals ( 16) if at least one of the values j 1 , . . ., j k is equal to zero.In particular, where the number of function arguments coincides with the integral multiplicity, and we can use the notation J (17) between iterated Stratonovich stochastic integrals are substituted into Expansion (14), then the unified Taylor-Stratonovich expansion is obtained [20].The difference is that Expansion ( 14) contains iterated stochastic integrals for functions (4) only, but some of them are mixed-type iterated stochastic integrals.In the unified expansion, there are no mixed-type integrals due to the transition to functions (6) under condition (7).These expansions are equivalent, but the latter expansion contains fewer types of iterated stochastic integrals.
Numerical methods based on the Taylor-Stratonovich expansion for solving stochastic differential equations can be proposed using the approach mentioned in the previous section [13,17,18,20].
To construct a numerical method with the order of mean-square or strong convergence p, it suffices to leave only the first terms on the right-hand side of Expansion (14).In particular, for the stochastic differential equation ( 13), we have the following explicit difference schemes: p = 1.0 (Milstein-type method) : p = 1.5 (Strong 1.5 order method) : p = 2.0 (Strong 2.0 order method) : where the leading term, which does not contain stochastic integrals, should be taken from Expansion (2) for the odd r = 2p.

Fundamentals of the Spectral Method
To obtain the main result, it is proposed to apply the spectral method [33,34].Here, it is described rather briefly, where only the notations and properties required below are given.
It is well known that all separable Hilbert spaces are isomorphic, i.e., it is possible to establish a one-to-one correspondence between the elements from separable Hilbert spaces.In this sense, the main space is the square-summable sequence space 2 [37].It is convenient to represent an element from 2 as an infinite column matrix in order to apply the matrix algebra techniques to it.
For example, let L 2 (H) be the space of square-integrable functions x(•) : with the standard inner product (•, •) L 2 (H) , and let {q(i, •)} ∞ i=0 be the orthonormal basis of L 2 (H).Then, any function x(•) ∈ L 2 (H) can be associated with the following infinite column matrix: , where The infinite column matrix X is called the spectral characteristic of the function x(•) with respect to the basis {q(i, •)} ∞ i=0 .Further, it is assumed that q(i, •) ∈ L ∞ (H), i = 0, 1, 2, . . ., where L ∞ (H) denotes the space of measurable and bounded functions Linear transformations of the element from the original Hilbert space are reduced to linear transformations of the corresponding element from 2 .If the original Hilbert space is L 2 (H), then a linear transformation of x(•) ∈ L 2 (H) corresponds to a linear transformation of its spectral characteristic X. Obviously, such transformations should be given by infinite matrices.Multilinear transformations of an ordered set of spectral characteristics are more complex, and they should be given by multidimensional infinite matrices with a dimension greater than 2. For example, a bilinear transformation is given by a three-dimensional infinite matrix.The elements of these infinite matrices generally depend on the basis {q(i, •)} ∞ i=0 .Let R be a multilinear operator such that then elements of the infinite (k + 1)-dimensional matrix R are defined as follows: i 1 , i 2 , . . . ,i k+1 = 0, 1, 2, . . ., where y i 2 ...i k+1 (•) = R(q(i 2 , •), . . ., q(i k+1 , •)).The infinite matrix R is called the spectral characteristic of the multilinear operator R with respect to the basis {q(i, •)} ∞ i=0 .Linear functionals defined on the Hilbert space correspond to linear functionals that are defined on 2 ; hence, it is also convenient to represent them by infinite column matrices.For multilinear functionals, it is proposed to use multidimensional infinite matrices with a dimension greater than 1.
Let F be a multilinear functional such that then elements of the infinite k-dimensional matrix F are given by the following equation and the infinite matrix or the column matrix F is called the spectral characteristic of the multilinear functional F with respect to the basis {q(i, •)} ∞ i=0 .Linear and multilinear operators and functionals need not be defined for all functions from L 2 (H).It suffices to consider them on some suitable linear subspaces.
For brevity, we introduce the notation S for the spectral transform to indicate the correspondence between functions, operators, functionals, and their spectral characteristics.Thus, By appropriately defining the multiplication of multidimensional matrices, we can reduce the linear and multilinear transformations of functions to operations with their spectral characteristics when using the matrix algebra techniques.An important point to consider is that algebraic operations should be carried out with infinite matrices, consequently, the problem of series convergence arises when we need to multiply matrices.However, it is solved quite simply if we consider only compact or bounded operators, as well as bounded functionals.
Next, we give some examples of the spectral characteristics of functions, linear and bilinear operators, and linear functionals, which are used below.
Consider the function k(t) ≡ 1, and indicate the elements of its spectral characteristic V 0 according to Equation ( 21): where ϑ is a parameter (which is omitted for notational simplicity).
The infinite column matrix V 0 is also the spectral characteristic of the linear functional J H , which associates the function x(•) ∈ L 2 (H) with its integral (bounded functional [37]): and this follows from Equation (23) As examples of linear operators, consider the integration operator D −1 (compact operator [37]) and the multiplication operator A ψ with bounded multiplier ψ(•) (bounded operator [38]): where elements of their spectral characteristics satisfy Equation (22) under condition k = 1: . .}, then it is natural to use the notation A n for the spectral characteristic of the multiplication operator A k n .Here, ϑ is also a parameter.
Another example is related to the bilinear multiplication operator M: and we apply Equation ( 22) under condition k = 2 to find the elements of its spectral characteristic: The right-hand side of the latter equality contains the product of the infinite threedimensional matrix and the infinite column matrix.The result of this operation is the infinite matrix: and other matrix operations, such as the multiplication of the infinite matrix and the infinite column matrix, are standard.
Additionally, we should highlight that if Ψ = S ψ(•) and for example, E = VV 0 and A = VF, where E is the infinite identity matrix (E is the spectral characteristic of the identity operator I).
The spectral method is convenient because operations with functions are reduced to operations with their spectral characteristics.For example, and Similar relations can also be written for random processes, as well as for random linear and multilinear operators and functionals.In this case, it is proposed to use the same basis {q(i, •)} ∞ i=0 .This approach leads to random spectral characteristics, and, here, it is natural to consider random processes, which have finite second moments and realizations that belong to L 2 (H) with probability 1.We use the notation L 2 (H) for the space of such random processes.Then, Equation ( 21) can be used to determine the elements of spectral characteristics of random processes from L 2 (H).For random linear and multilinear operators and functionals, we use Equations ( 22) and (23).
Consider the example of the random linear functional J W H that associates the function x(•) ∈ L 2 (H) with its stochastic integral over the Wiener process W(•): To determine the elements of its spectral characteristic V, we apply Equation ( 23) under condition k = 1, i.e., According to stochastic integral properties, V i are independent random variables that have a standard normal distribution (the normal distribution law is a consequence of the integration over the Gaussian random process, and the independence follows from the Itô isometry) [35,36].Then for example, As an example of the random linear operator, consider the stochastic integration operator D −1 W over the Wiener process W(•): where Equation ( 22) is applied under condition k = 1 to obtain the elements of its spectral characteristic: , where V corresponds to the spectral characteristic of the random linear functional J W H .It is known that the Wiener process W(•) is not differentiable in the ordinary sense.However, in the generalized sense, it is admissible to consider its derivative V(•), which is called Gaussian white noise [35].More strictly, Gaussian white noise is the generalized random process, i.e., the random linear functional, and this is precisely the functional J W H . Therefore, the spectral characteristic V is also called the spectral characteristic of Gaussian white noise V(•), i.e., S[J In [33,34], the spectral characteristic W of the Wiener process W(•) is proposed in the form W = P −1 V.
The random linear operator D −1 W and the random linear functional J W H can also act on the elements from L 2 (H), but under the condition that stochastic integrals are understood as Stratonovich stochastic integrals.This is relevant as it is this that only provides the ordinary integration rules and the possibility of applying the spectral method.
Corresponding random linear operators and functionals are denoted by S D −1 W and S J W H .In fact, S D −1 W = D −1 W and S J W H = J W H .For Itô stochastic integrals, a preliminary transformation to Stratonovich stochastic integrals is required [17], after which we can apply the spectral method.Corresponding random linear operators and functionals are denoted by I D −1 W and I J W H .In general, they differ from S D −1 W and S J W H , where, for example

Spectral Representation of Iterated Stochastic Integrals
In this section, we formulate the main result related to the spectral representations of iterated stochastic integrals.Since it is more natural to apply the spectral method to Stratonovich stochastic integrals, we start with their representation.Theorem 1.Let V 0 be the spectral characteristic of the function k(t) ≡ 1, and let V 1 , . . ., V s be the spectral characteristics of Gaussian white noises V 1 (•), . . ., V s (•), which correspond to the Wiener processes W 1 (•), . . ., W s (•).Let P −1 be the spectral characteristic of the integration operator, and let V be the spectral characteristic of the bilinear multiplication operator.Moreover, let Ψ 1 , . . ., Ψ k be the spectral characteristics of multiplication operators with multipliers ψ 1 (•), . . ., ψ k (•).The spectral characteristics V 0 , P −1 , V, Ψ 1 , . . ., Ψ k are defined with respect to the basis {q(i, •)} ∞ i=0 .Then, the iterated Stratonovich stochastic integral (16) is represented by relations or in the explicit form Proof of Theorem 1.The introduction of the notation for the spectral characteristics of functions and random processes: Moreover, and the representation P −1,V j = P −1 (VV j ) holds for the spectral characteristics P −1,V j .It is proved in [39], but it can also be applied for j = 0 since , where E is the infinite identity matrix.
The iterated Stratonovich stochastic integral S J W(j 1 ...j k ) H k ψ (•) may be represented as the composition of simple operations including deterministic (j l = 0) or stochastic (j l = 0) integrations, as well as the multiplication by weight functions ψ l (•), l = 1, . . ., k: where If the spectral transform is applied to both sides of Equation ( 27), then This means that since Ψ l and VV j l are symmetric matrices (the product of symmetric matrices is commutative), and the latter expression can be simplified under condition l = 1: The last step in the proof is to represent the random linear functional S J Therefore, we have the spectral analog of Equation ( 26): i.e., Relations (24) are proved.Finally, excluding X 1 , . . ., X k−1 from them, we prove Relation (25).
If we use the unified Taylor-Stratonovich expansion from [20], then we need the following iterated Stratonovich stochastic integrals and their spectral representation: When comparing Relations ( 17), (28), and (29), we derive the following useful relations: Moreover, we have Note that Relations ( 30) and ( 31) can be proved using only the properties of the spectral characteristics V 0 , F, P −1 , A, and V.However, in this case, it suffices to prove them indirectly by comparing iterated Stratonovich stochastic integrals and their spectral representation.
Thus, the final form for the spectral representation of iterated stochastic integrals from Expansion ( 14) is i.e., all the necessary iterated stochastic integrals to be simulated for the implementation of numerical methods ( 18)-( 20) for solving the stochastic differential equation ( 13) are expressed using algebraic operations with column matrices V 0 , V 1 , . . ., V s , F, matrices P −1 and A, and the three-dimensional matrix V.The spectral characteristics V 1 , . . ., V s are random, but other spectral characteristics are deterministic.An important point is that the spectral characteristics V 0 , F, P −1 , A, and V depend on the values ϑ and h.However, the property holds for iterated Stratonovich stochastic integrals, where is the function ( 6) under condition (7) with t 0 = 0, and δ j l 0 is the Kronecker delta, l = 1, . . ., k.Therefore, it suffices to find the spectral characteristics V 0 , F, P −1 , A, and V with respect to the basis {q( ).Further, we consider the spectral representation of iterated Itô stochastic integrals (5).Theorem 2. Let the conditions of Theorem 1 hold.Then, the Itô iterated stochastic integral (5) is represented by relations where δ * j l−1 j l = δ j l−1 j l (1 − δ j l 0 ), δ j l−1 j l and δ j l 0 are Kronecker deltas, l = 2, . . ., k.
Proof of Theorem 2. Using the approach from the proof of Theorem 1, we may represent the iterated stochastic integral I J W(j 1 ...j k ) H k ψ (•) as the composition of simple operations, including deterministic (j l = 0) or stochastic (j l = 0) integrations, as well as the multiplication by weight functions ψ l (•), l = 1, . . ., k: The next step is to express I J W(j 1 ...j k ) H k ψ (•) via Stratonovich stochastic integrals.For l = 1, it is necessary to change the notation only since the Itô and Stratonovich stochastic integrals of deterministic functions coincide, i.e., and for l > 1, it is required to use the relationship between the Itô and Stratonovich stochastic integrals [17]: and The use of the values δ * j l−1 j l instead of δ j l−1 j l is due to the fact that mixed-type iterated stochastic integrals are considered, i.e., not only the stochastic integration, but also the deterministic integration is allowed, in which the relationship between the Itô and Stratonovich stochastic integrals is not used.
Equation (35) is the same as for the iterated Stratonovich stochastic integral since if l = 1, then the integration is performed for the deterministic function.Therefore, If the spectral transform is applied to both sides of Equation ( 36) when using the notations from the proof of Theorem 1 and the linearity property (S is the linear transform), then The expression for the second term on the right-hand side of the latter equality is obtained in the proof of Theorem 1, so here we can restrict ourselves to the first term only: It remains to represent the random linear functional I J W(j k ) H , but its value includes S J W(j k ) H A ψ k X k−1 (•), which is obtained in the proof of Theorem 1.Therefore, it suffices to consider the first term on the right-hand side of Relation (37): Thus, Relations (34) are proved.
If we use the unified Taylor-Itô expansion from [20], then we need the following iterated Itô stochastic integrals and their spectral representation: On the right-hand sides of the above relations, all of the terms containing spectral characteristics V j 1 , . . ., V j 4 are iterated Stratonovich stochastic integrals.This follows from Theorem 1, and we can use Relations (30) and (31) for them.So, we obtain the final form Therefore, iterated stochastic integrals needed for the implementation of numerical methods ( 9)-( 12) for solving the stochastic differential equation ( 1) are expressed using algebraic operations with the column matrices V 0 , V 1 , . . ., V s , F, matrices P −1 and A, and the three-dimensional matrix V.
The same property holds for iterated Itô stochastic integrals as for iterated Stratonovich stochastic integrals, i.e., where all of the notations are introduced as above.Thus, if the spectral characteristics V 0 , F, P −1 , A and V, with respect to the basis {q(i, •)} ∞ i=0 of L 2 ([0, 1]), are known, then we have the base of simulation for all the necessary iterated Itô stochastic integrals.
Obviously, Properties (33) and ( 40) significantly reduce the computational costs for finding spectral characteristics.In addition, they allow one to determine the terms from Expansions (2) and ( 14) that should be included in the numerical scheme to guarantee the order of convergence p: Equations ( 32) and (39) give the exact representation, but, in practice, only the approximate simulation of iterated stochastic integrals is possible, except for some integrals.The approximate simulation requires one to use Equations ( 32) and ( 39), but all matrices should be truncated with some chosen order.
In this paper, the accuracy of the mean-square approximation of iterated stochastic integrals is not considered since this requires a great deal of additional theoretical results.However, it is possible to find the exact values or the estimates of approximation errors for all iterated stochastic integrals from Expansions (2) and (14).The results from [20,31] can be applied to this.For example, consider the iterated stochastic integrals of a second multiplicity I J W(j 1 j 2 ) [ϑ,ϑ+h] k(•) and S J W(j 1 j 2 ) [ϑ,ϑ+h] k(•) under conditions j 1 = j 2 and j 1 , j 2 = 0.The exact representation of them is written as Consequently, the approximate representation has the form where L is the truncation order and .
Then, we have where L 2 (H 2 ) is the space of square-integrable functions x(•) : . But this example is one of the simplest since the multiplicity k = 2 is not high and the iterated stochastic integrals I J In [31], Legendre polynomials, cosines, Walsh and Haar functions, and trigonometric functions are tested as the basis {q(i, •)} ∞ i=0 , and cosines (the Fourier series for an even function is based on this basis) provide the best approximation.The use of Walsh and Haar functions is similar to numerical integration [40].

Numerical Experiments
Here, we use the obtained spectral representations for the simulation of iterated Itô and Stratonovich stochastic integrals of the multiplicity k = 2, 3, 4 under condition H = [0, 1].
As the basis {q(i, •)} ∞ i=0 of L 2 ([0, 1]), we choose cosines and, for this basis, we have the following spectral characteristics of functions, operators, and functionals [30]: (1) The spectral characteristic V 0 of the function k(t) ≡ 1: (2) The spectral characteristic F of the function k 1 (t) = t: (3) The spectral characteristic P −1 of the integration operator D −1 : (4) The spectral characteristic A of the multiplication operator A k 1 : (5) The spectral characteristic V of the bilinear multiplication operator M: The above relations for spectral characteristics provide an approximate simulation of iterated Itô and Stratonovich stochastic integrals of arbitrary multiplicity from Expansions (2) and ( 14), or from the corresponding unified expansions.
For each truncation order L, we obtain the realizations of both iterated Itô stochastic integrals {I m } M m=1 and iterated Stratonovich stochastic integrals {S m } M m=1 , where M = 10 6 is the sample size.The expectation and the second moment may be estimated by formulas The estimates of the expectation and the second moment are given in Tables 1 and 2 (in these tables, the symbol * means "pairwise distinct"), respectively.The results for the iterated stochastic integrals of a fourth multiplicity are partially presented since the total number of different cases is too large, but the remaining integrals are similar in their properties.
For the iterated stochastic integrals, there are no analytical expressions to describe their distribution laws, except for some integrals.However, it is quite simple to find the expectation and the second moment, so we can compare the estimates of the expectation and the second moment with the known exact values.All the results are consistent with the exact values, which are given in the last column of Tables 1 and 2. In general, the estimation results contain errors of two types.The first error is the methodological error of the spectral method, and it is associated with the truncation of all matrices with the order L for any dimension.The second error is the Monte Carlo error, which depends on the sample size M only.
For example, the first error dominates for rows k = 2 (I, S) j 1 = j 2 , k = 3 (I, S) j 1 , j 2 , j 3 * , and k = 4 (I, S) j 1 , j 2 , j 3 , j 4 * in Table 2; however, the convergence to the exact values is obvious.There are many such rows in Tables 1 and 2. But for some rows, there exist only a second error, for example, this is true for rows k = 2 (I, S) j 1 = j 2 , k = 3 (I, S) j 1 , j 2 , j 3 * , and k = 4 (I, S) j 1 , j 2 , j 3 , j 4 * in Table 1.This is also true for some other rows of Tables 1 and 2; however, there, a more detailed analysis is required.

Discussion
This paper considers numerical methods for solving stochastic differential equations with the order of mean-square and strong convergence p = 0.5, 1.0, 1.5, 2.0; however, using the obtained spectral representations, methods with the order of convergence 2.5 and more can also be considered.There are no restrictions on both the multiplicity of the iterated stochastic integrals and the order of convergence of numerical methods.
For the Euler-Maruyama method, the spectral method is certainly not required (the strong 0.5 order method).For Milstein-type methods (the strong 1.0 order methods), the described approach gives the matrix notation and the possibility of using different bases only.The advantage of the spectral representations appears for methods with the order of convergence 1.5 and more.The approximate simulation of iterated Itô and Stratonovich stochastic integrals of arbitrary multiplicity uses matrix algebra techniques, and it therefore has a very simple implementation.In addition, there is no need to expand functions k(•) and k n 1 ...n k (•) with different k into generalized multiple Fourier series, but it is sufficient to use the spectral characteristics V 0 , F, P −1 , A, and V.
Through using the proposed approach, we can obtain a similar spectral representation of the iterated stochastic θ-integrals.From this more general representation, Theorems 1 and 2 are obtained as two particular cases: θ = 0 for iterated Itô stochastic integrals, and θ = 1/2 for iterated Stratonovich stochastic integrals.

Conclusions
The main result is the spectral representations of iterated Itô and Stratonovich stochastic integrals for the implementation of numerical methods for solving Itô and Stratonovich stochastic differential equations with high orders of mean-square and strong convergence.Their advantage is provided by simple final forms, namely each iterated stochastic integral is represented by matrix algebra techniques with known matrices.For the approximate simulations of iterated Itô and Stratonovich stochastic integrals, the standard mathematics software are sufficient, see -for example -the software based on the BLAS (Basic Linear Algebra Subprograms) specification or the software for mathematical calculations, such as Mathcad, MATLAB, etc.

Table 1 .
The estimates of the expectation of iterated Itô (I) and Stratonovich (S) stochastic integrals.

Table 2 .
The estimates of the second moment of iterated Itô (I) and Stratonovich (S) stochastic integrals.