Projected particle methods for solving McKean-Vlasov stochastic differential equations
Denis Belomestny John Schoenmakers




Projected particle methods for solving McKean-Vlasov stochastic differential equations
    Denis Belomestny^1 and John Schoenmakers ^2
    December 30, 2023
======================================================================================




We propose a novel projection-based particle method for solving McKean-Vlasov
stochastic differential equations. Our approach is based on a projection-type
estimation of the marginal density of the solution in each time step. The
projection-based particle method leads in many situation to a significant
reduction of numerical complexity compared to the widely used kernel density
estimation algorithms. We derive strong convergence rates and rates of density
estimation. The convergence analysis, particularly in the case of linearly
growing coefficients, turns out to be rather challenging and requires some new
type of averaging technique. This case is exemplified by explicit solutions to
a class of McKean-Vlasov equations with affine drift. The performance of the
proposed algorithm is illustrated by several numerical examples.



[1]Duisburg-Essen University, Thea-Leymann-Str. 9, D-45127
Essen, Germany,  [2]
Weierstrass Institute for Applied Analysis and Stochastics, Mohrenstr. 39,
10117 Berlin, Germany,  Keywords:
Mckean-Vlasov equations, particle systems, projection estimators, explicit
solutions.2010 MSC: 60H10, 60K35.



§ INTRODUCTION


Nonlinear Markov processes are stochastic processes whose transition functions
may depend not only on the current state of the process but also on the
current distribution of the process. These processes were introduced by
McKean <cit.> to model plasma dynamics. Later nonlinear Markov
processes were studied by a number of authors; we mention here the books of
Kolokoltsov <cit.> and
Sznitman <cit.>. These processes arise naturally in the
study of the limit behavior of a large number of weakly interacting Markov
processes and have a wide range of applications, including financial
mathematics, population dynamics, and neuroscience (see, e.g.,
<cit.> and the references therein).

Let [0,T] be a finite time interval and (Ω,ℱ,P) be
a complete probability space, where a standard m-dimensional Brownian motion
W is defined. We consider a class of McKean-Vlasov SDEs, i.e. stochastic
differential equation (SDE) whose drift and diffusion coefficients may depend
on the current distribution of the process of the form:

    {[                                                           X_t =ξ+∫_0^t∫_ℝ^da(X_s,y)μ_s(dy)ds+∫
    _0^t∫_ℝ^db(X_s,y)μ_s(dy)dW_s;                                                           μ_t                                          =Law(X_t),   t∈0,T], ].

where X_0=ξ is an ℱ_0-measurable random variable in
ℝ^d, a: ℝ^d×ℝ^d→ℝ^d and b: ℝ^d×ℝ^d→ℝ^d× m. If the functions a and b are smooth with
uniformly bounded derivatives and the random variable ξ has finite moments
of any order, then (see <cit.>) there is a unique strong
solution of (<ref>) such that for all p>1,

    E[  sup_s≤ T|X_s|^p]  ≤∞.

In the sequel we assume that there exists a unique strong solution of
(<ref>) such that (<ref>) holds and refer to
<cit.> for more general sufficient conditions for this.

Assume that d=1 and for any t≥0, the measure μ_t(du) possesses a
bounded density μ_t(u). Then the family of these densities satisfies a
nonlinear Fokker-Planck equation of the form

    ∂μ_t(x)/∂ t     =-∂/∂ x(
    (  ∫ a(x,y)μ_t(y) dy)   μ_t(x)) 
         +1/2∂^2/∂ x^2(  (  ∫
    b(x,y)μ_t(y) dy)  ^2 μ_t(x))  ,

which can be seen as an analogue of a well-known linear Fokker-Planck equation
in the case of linear stochastic differential equations. In
Section <ref> we will show that if the drift a is affine in x, and
the diffusion coefficient b is independent of x, then the system
(<ref>), and hence (<ref>), has an explicit solution. These
solutions, apart from being interesting in their own right, also provide
explicit cases of an explosive behavior.

The theory of the propagation of chaos developed in <cit.>,
states that (<ref>) is a limiting equation of the system of stochastic
interacting particles (samples) with the following dynamics

    X_t^i,N=ξ^i+∫_0^t∫_ℝ^da(X_s^i,N,y)μ
    _s^N(dy) ds+∫_0^t∫_ℝ^db(X_s^i,N,y)μ_s
    ^N(dy) dW_s^i

for i=1,…,N, where μ_t^N=1/N∑_i=1^Nδ
_X_t^i,N, ξ^i, i=1,…,N, are i.i.d copies of ξ,
distributed according the law μ_0, and W^i, i=1,...,N, are
independent copies of W. In fact it can be shown, under sufficient
regularity conditions on the coefficients, that convergence in law for
empirical measures on the path space holds, i.e., μ^N={μ_t^N
:t∈0,T]}→μ, N→∞, see
<cit.>.

Despite the numerous branches of research on stochastic particle systems,
results on numerical approximations of McKean-Vlasov-SDEs are very sparse. The
authors in <cit.> proposed to use the Euler scheme with
time-step h=T/L, that for l=0,…,L-1, yields

    X̅_t_l+1^i,N=X̅_t_l^i,N+1/N∑_j=1^N
    a(X̅_t_l^i,N,X̅_t_l^j,N) h+1/N∑_j=1
    ^Nb(X̅_t_l^i,N,X̅_t_l^j,N) Δ_l+1W^i

for i=1,…,N, t_l=hl, and Δ_l+1W^i=W_h(l+1)^i
-W_hl^i (see also <cit.> for more general MVSDEs
and <cit.> for Gauss-quadrature based approach).
Implementation of the above algorithm requires usually N^2× L
operations in every step of the Euler scheme. By using the algorithm presented
here one can significantly reduce the complexity of the particle simulation
especially if the coefficients of the corresponding McKean-Vlasov SDE are
smooth enough.

The contribution of this paper is twofold. On the one hand, we propose a new
approximation methodology based on a projection-type estimation of the
marginal densities of (<ref>). This methodology often leads to
numerically more efficient algorithms than the kernel-type approximation
algorithms, as they can profit from a global smoothness of coefficients
a,b and the corresponding marginal densities. On the other hand, we
present a comprehensive convergence analysis of the proposed algorithms in the
case of possibly linearly growing (in x) coefficients a and b. To the
best of our knowledge, no stability analysis of MVSDEs under this linear
growth assumption was done before. In fact such analysis is rather challenging
and requires a special type of averaging technique. And, last but not least,
we study a general class of MVSDEs with affine drift and derive their explicit
solutions, to the best of our knowledge, for the first time.



The paper is organized as follows. In Section <ref> we present the
idea of our projected particle method. Section <ref> is devoted to
the convergence analysis of the projected particle method. In particular, in
Section <ref> we derive the convergence rates for the corresponding
projected density estimate. Section <ref> presents a thorough study
of affine MVSDEs. Numerical examples for affine and convolution-type MVSDEs
are presented in Section <ref>. All proofs are collected in
Section <ref>.



§ PROJECTED PARTICLE METHOD


 Let w:ℝ^d→ℝ_+ be some weight
function with w>0, such that

    a(x,·),b(x,·)∈ L_2(ℝ^d,w)  for any 
    x∈ℝ^d.

Let further (φ_k, k=0,1,2,..) be a total orthonormal system in
L_2(ℝ^d,w). The corresponding (generalized) Fourier
coefficients of the functions a(x,·) and b(x,·) are given by

    α_k(x)       :=∫ a(x,u)φ_k(u) w(u) du∈ℝ
    ^d,
    β_k(x)       :=∫ b(x,u)φ_k(u) w(u) du∈ℝ^d×
    m

and the following series representation holds

    a(x,·)=∑_k=0^∞α_k(x)φ_k(·) and
      b(x,·)=∑_k=0^∞β_k(x)φ_k(·),  x∈ℝ^d,

in L_2(ℝ^d,w). Further it is assumed that each function
φ_k is bounded so that the functions

    γ_k(s):=E[  φ_k(X_s)]

are well defined. Now let us fix some natural number K>0 and consider a
projected particle approximation for (<ref>)

    X_t^i,K,N=ξ^i+∫_0^t∑_k=0^Kγ_k^N(s)α
    _k(X_s^i,K,N) ds+∫_0^t∑_k=0^Kγ_k^N(s)β
    _k(X_s^i,K,N) dW_s^i

for i=1,…,N, where

    γ_k^N(s):=1/N∑_j=1^Nφ_k(X_s^j,K,N)

can be regarded as an approximation to (<ref>). The projected system
(<ref>), with (<ref>), is heuristically motivated by assuming
that for any s≥0 the measure μ_s(du) possesses a density μ
_s(u) that, in view of

    γ_k(s)=∫μ_s(u)φ_k(u)du

(cf. (<ref>)), formally satisfies

    μ_s(u)=∑_k=0^∞γ_k(s)φ_k(u)w(u).

Then, we (formally) have the expansion

    ∫ a(x,u)μ_s(u)du=∑_k=0^∞α_k(x)γ_k(s),

and this motivates the drift term in (<ref>). For the diffusion
term in (<ref>) an analogue motivation applies. In order to solve
(<ref>) we may consider, for any fixed L>0, an Euler-type
approximation,

    X̅_t^i,K,N=X̅_η(t)^i,K,N+∑_k=0^Kγ_k
    ^N(  η(t))   α_k(X̅_η(t)^i,K,N
    ) (t-η(t))
    
    +∑_k=0^Kγ_k^N(  η(t))   β_k(X̅_η(t)^i,K,N) (W_t^i-W_η(t))

for i=1,…,N, and h=T/L, where η(t):=lh for t∈
lh,(l+1)h), l=1,…,L. Note that in order to generate a discretized
particle system (X̅_hl^i,K,N), i=1,…,N, l=1,…,L, we
need to perform (up to a constant depending on the dimension) NLK
operations. This should be compared to N^2L operations in
(<ref>). Thus if K is much smaller than N, we get a significant
cost reduction. Of course, this complexity analysis implicitly assumes that
the generalized Fourier coefficients α_k(x) and β_k(x) are
known in closed form or can be cheaply computed. For more details in this
respect see Remark <ref> below. 


 Many well known McKean-Vlasov type models used in physics and
engineering are constructed and formulated via certain Fourier type
expansions of the respective drift and/or diffusion coefficients. For
example, in the famous Kuramoto-Shinomoto-Sakaguchi model (see e.g.
<cit.>, eq. (5.214)) or in the coupled Brownian phase
oscillators (see <cit.>) the mean field potential is
given by its Fourier series, which entails a similar expansion for the
coefficient a(x,u) (b(x,u) is constant).  Let us also mention
a classical work of <cit.>, where a known power series
expansion for the coefficients of a nonlinear Fokker-Planck equation is
assumed. From another point of view, since the basis ( φ
_k) with the corresponding weight w can, in principle, be chosen
freely, it is natural to assume  they can be chosen such that the coefficients α_k(x) and
β_k(x) can be computed in closed form. In this respect, let us give
some further examples. If for any x, a(x,·) is a linear combination of
functions of the form:

    q_1(u_1)··· q_d(u_d)

where each q_i:ℝ→ ℝ is a polynomial with coefficients
possibly depending on x, then (  φ_k) may be taken to
be Hermite functions in ℝ^d, i.e.

    φ_α(  u)  =H_α_1(u_1)···
    H_α_d(u_d)e^-| u| ^2/2,  α
    =(α_1,…,α_d).

The latter situation appears for instance in the popular interaction case with
a(x,u)=A(x-u), where the function A has a given representation

    A(z)=∑_αc_αz_1^α_1·…· z_d
    ^α_d,   z∈ℝ^d,  α∈ℕ_0
    ^d.

As another example, note that the Fourier coefficients of any function of the
form

    u→ u_1^α_1·…· u_d^α_d
    e^-| u-c| ^2/σ,   c∈ℝ
    ^d,  α∈ℕ_0^d,σ>0,

with respect to the Hermite basis above can be expressed in closed form. One
so could also consider a(x,u), b(x,u) of the form

    ∑_r=1^Rq_r(x)u^α_re^-| u-c_r(x)|
    ^2/σ_r(x),

with free to choose q_r(x)∈ℝ, c_r(x)∈ℝ^d,
σ_r(x)∈ℝ_+, α_r∈ℕ_0^d,
R∈ℕ.




§ CONVERGENCE ANALYSIS


 In this section we first study the convergence of the
approximated particle system (<ref>) to the solution of the
original system (<ref>). As a first obvious but important observation,
we note that the distribution of the triple (  X_s^j,K,N,X_s
^K,N,X_s^j) with X_s^K,N:=(  X_s^1,K,N
,…,X_s^N,K,N) does not depend on j, and therefore we can
write

    (  X^j,K,N,X^K,N,X^j)  distr.=(
    X^·,K,N,X^K,N,X^·)   for j=1,...,N.

For ease of notation, henceforth we denote with |·|
:=|·| _ for a generic dimension  the
standard Euclidian norm in ℝ^. Let us make the following assumptions.



  (AF) The basis functions (φ_k) fulfil

    |φ_k(z)-φ_k(z^')|≤
    L_k,φ| z-z^'| ,  |φ
    _k(z)|≤ D_k,φ,   k=0,1,…

for all z,z^'∈ℝ^d and some constants L_k,φ,D_k,φ>0.


  (AC) The functions α_k(x), β_k(x), k=0,1,2,…
satisfy

    |α_k(x)|   ≤ A_k,α(1+|
    x| )  with   (  A_k,α)  _k=0,1,…∈ l_2,
    ∑_k=0^∞D_k,φA_k,α   ≤ D_φA_α,  and   ∑_k=0^∞L_k,φA_k,α≤
    L_φA_α,
    |β_k(x)|   ≤ A_k,β(1+|
    x| )  with   (  A_k,β)  _k=0,1,…∈ l_2,
     ∑_k=0^∞D_k,φA_k,β   ≤ D_φA_β,  and   ∑_k=0^∞L_k,φA_k,β≤
    L_φA_β,

for some constants A_α, A_β, D_φ, and L_φ>0, and further

    sup_x,x^'∈ℝ^d,x≠ x^'|α
    _k(x)-α_k(x^')|/|x-x^'|   ≤ B_k,α  with   ∑_k=0^∞D_k,φB_k,α≤ D_φB_α,
    sup_x,x^'∈ℝ^d,x≠ x^'|β
    _k(x)-β_k(x^')|/|x-x^'|   ≤ B_k,β   with   ∑_k=0^∞D_k,φB_k,β≤ D_φB_β,

for some B_α, B_β >0.

  (_p) For some p>0 the initial distribution μ_0 possesses a finite absolute moment of order p.



In the sequel, for any random variable ξ∈ℝ^ on
(Ω,ℱ,P) we shall use ‖ξ‖_p for the
norm of |ξ| in L_p(Ω). The following bound
on the strong error can be proved.


 For p≥2, it holds under assumptions (AC), (AF) and (AM_p)
that

    ‖sup_0≤ r≤ T| X_r^·,K,N-X_r^·|‖ _p   ≲ N^-1/2+∑_k=K+1^∞A_k,α ‖γ_k‖_L_p[0,T]
         +∑_k=K+1^∞A_k,β ‖γ_k‖_L_p[0,T],

where ≲ stands for an inequality with some (hidden) positive finite
constant depending only on A_α,A_β,B_α,B_β
,D_φ,L_φ, p, and T.



For 1≤ p^'≤2, we simply have

    ‖sup_0≤ r≤ T| X_r^·,K,N-X_r^·|‖ _p^'≤‖sup_0≤ r≤
    T| X_r^·,K,N-X_r^·|‖ _p

for any p≥2.


The next theorem, on the convergence of the Euler approximation
(<ref>) to the projected system (<ref>), can be
proved along the same lines as the proof of Theorem <ref>.


For p≥2, it holds under assumptions (AC), (AF) and (AM_p) that for any
natural K,N

    ‖sup_0≤ r≤ T|X̅_r^·,K,N-X_r
    ^·,K,N|‖ _p≲√(h),

where ≲ stands for an inequality with some (hidden) positive finite
constant depending only on A_α,A_β,B_α,B_β,D_φ,L_φ, p and T.




  
Discussion

The bound (<ref>) is proved under rather general assumptions on
the coefficients a(x,y) and b(x,y). In particular, we allow for linear
growth of these coefficients in x. This makes the proof of the bound in
Theorem <ref> rather challenging, since we need to avoid an explosion. In
order to overcome this problem, we employ a kind of averaging technique which,
being combined with the symmetry of the particle distribution and the
existence of moments (see Section <ref>), gives the desired bound.
Note that for this we have to assume existence and uniqueness of a strong
solution of the original MVSDE (<ref>). Funaki <cit.>
proved existence and uniqueness under (essentially) global Lipschitz
condition. However, one should be able to extend his results by exploiting a
kind of one sided Lipschitz condition like in <cit.> or
<cit.>.

The bound (<ref>) consists of stochastic and approximation
errors. While the first error is of order 1/√(N), the second one depends
on K and the properties of the coefficients a(x,y) and b(x,y). If these
coefficients are smooth in the sense that their generalized Fourier
coefficients (α_k) and (β_k) decay fast, then the
approximation error can be made small even for medium values of K.


The (normalized) Hermite polynomial of order j is given, for j≥0, by

    H_j(x)=c_j(-1)^je^x^2d^j/dx^j(e^-x^2
    ),   c_j=(  2^jj!√(π))  ^-1/2.

These polynomials satisfy: ∫_ℝH_j(x)H_ℓ(x)e^-x^2dx=δ_j,ℓ and, as a consequence,

    φ_k(u)=H_k(u)e^-u^2/2,   k=1,2,…,

is a total orthonormal system in L_2(  ℝ^d) (i.e.
here w=1). Moreover, (φ_k)_k≥0 fulfil the assumption (AF)
with D_k,φ and L_k,φ being uniformly bounded in k, see,
e.g. <cit.>, p. 242. Now let us suppose that
a(x,·),b(x,·)∈ L_2(ℝ^d) for any x∈ℝ, and
discuss the assumptions (AC).



Suppose that for any x∈ℝ, the functions (in u)

    a(x,u):=a(x,u)/√(1+x^2),   b(x,u):=b(x,u)/√(1+x^2)

admit derivatives in u up to order s>2 such that the functions (in u)

    u^ℓ∂_u^ma(x,u),   u^ℓ∂_u
    ^mb(x,u),   0≤ l+m≤ s

are bounded and belong to L_1(ℝ) (uniformly in x) together
with their first derivatives in x. Then the assumption (AC) is satisfied
and

    ‖sup_0≤ r≤ T| X_r^·,K,N-X_r^·|‖ _p≲ K^1-s/2+N^-1/2,

as K,N→∞.



We have

    α_k(x)       =√(1+x^2)∫a(x,u)H
    _k(u) e^-u^2/2 du =√(1+x^2) α_k(x) ,
    β_k(x)       =√(1+x^2)∫b(x,u)H
    _k(u) e^-u^2/2 du= √(1+x^2) β_k(x).

The identity

    (2k+2)^1/2H_k(x)=H^'_k+1(z)

and the integration-by-parts formula imply

    α_k(x)       =.  a(x,u)e^-u^2
    /2H_k+1(u)/(2k+2)^1/2|  _-∞^∞
         -1/(2k+2)^1/2∫_-∞^∞[  ∂a(x,u)/∂ u-ua(x,u)]  H
    _k+1(u)e^-u^2/2 du.

Note that |H_k(u)| e^-u^2/2≤1 uniformly in u and k
(see, e.g. <cit.>, p. 242.) Hence if a(x,u) is bounded and

    ∫|  ∂a(x,u)/∂ u-ua(x,u)|   du

is bounded uniformly in x, then α_k(x)=O(
k^-1/2) uniformly in x. The second integration-by-parts shows
that α_k(x)=O(  k^-1)  , provided the
functions

    u^2·a(x,u),∂^2a(x,u)/∂
    u^2,u·∂a(x,u)/∂ u

are integrable on ℝ with their L_1(ℝ) norms uniformly
bounded in x. Integrating by parts further, we derive the desired statement.



As a rule, one chooses N and K such that the errors in (<ref>) are
balanced, that is N^1/(s-2)∼ K, yielding a proportional reduction of
computational cost of order N· K/N^2 ∼ N^-(s-3)/(s-2).
Alternatively we can compare the complexity, that is the computational cost
for achieving a prescribed accuracy ε, for of the Euler schemes
(<ref>) and (<ref>). It is not difficult to see that, after incorporating the
path-wise time discretization error, the standard Euler scheme (<ref>) has complexity of
order ε^-6, while the projected one (<ref>) has complexity of
order ε^-(4s-6)/(s-2) which is significantly smaller when s>3.
Moreover, in <cit.> conditions are formulated, guaranteeing
that all measures μ_t, t≥0, possess infinitely smooth exponentially
decaying densities. In this case we can additionally profit from the decay of
the generalized Fourier coefficients (γ_k) such that the convergence
rates in (<ref>) give rise to a proportional reduction of
computational cost approaching N^-1, corresponding to a complexity of
order ε^-4 (modulo some logarithmic term) for the method (<ref>).




 §.§ Density estimation


 Let us now discuss the estimation of the densities μ_t,
t≥0. Let us assume that the formal relationship (<ref>) holds in
the sense that

    μ_s/w=∑_k=0^∞γ_k(s)φ_k

in L_2(  ℝ^d,w)  , i.e. μ_s^2/w∈
L_1(  ℝ^d)  . Fix some t>0, K_test
∈ℕ and set

    μ_t^K_test,K,N(x):=∑_k=1^K_test
    γ_k^N(t)φ_k(x)w(x)

with γ_k^N(t):=1/N∑_i=1^Nφ_k(X_t^i,K,N),
k=1,…,K_test. We obviously have

    E∫|μ_t^K_test,K,N(x)-μ_t
    (x)|^2w^-1(x) dx       =∑_k=1^K_testE[
    |γ_k^N(t)-γ_k(t)|^2] 
         +∑_k=K_test+1^∞|γ_k(t)|^2,

where (due to (AF))

    E[  |γ_k^N(t)-γ_k(t)|^2]      
    =E[  |1/N∑_j=1^Nφ_k
    (X_t^j,K,N)-E[  φ_k(X_t^·)]
    | ^2] 
       ≤2E[  |1/N∑_j=1^N(
    φ_k(X_t^j,K,N)-φ_k(X_t^j))  |
    ^2] 
         +2E[  |1/N∑_j=1^N(  φ
    _k(X_t^j)-E[  φ_k(X_t^j)]  )
    | ^2] 
       ≤2L_k,φ^2E[  | X_t^·,K,N
    -X_t^·| ^2]  +2/NVar[
    φ_k(X_t)]  ,

since the X^j are independent. Theorem <ref> now implies

    (  E∫|μ_t^K_test
    ,K,N(x)-μ_t(x)|^2 w^-1(x) dx)  ^1/2≲(
    1/N∑_k=1^K_test(L_k,φ^2+D_k,φ
    ^2))  ^1/2
    
    +(  ∑_k=1^K_testL_k,φ^2)  ^1/2(
    ∑_k=K+1^∞(A_k,α+A_k,β) ‖γ_k‖
    _L_p[0,T]) 
    
    +(  ∑_k=K_test+1^∞|γ_k(t)|^2)
    ^1/2.

The last term always converges to zero as K_test→∞,
since μ_t/w∈ L_2(ℝ^d,w). The first term can be controlled
for any fixed K_test by taking N large enough. Finally, for any
fixed K_test, the second term can be made small by taking K large
enough and using the condition (AC).



§ SPECIFIC MODELS






 §.§ Generalized Shimizu-Yamada Models




Inspired by the work of Shimizu and Yamada <cit.>,
<cit.> and <cit.>, we consider one-dimensional
McKean-Vlasov equations of the form (<ref>) with

    a(x,u):=a^0(u)+a^1(u)x,   b(x,u):=b(u).

This class of models allows for a linear dependence of drift on the
distribution of X through E[  a^0(X_t)] and
E[  a^1(X_t)]  . Let us define for polynomially
bounded and measurable functions a^j and b the generalized Gauss
transforms,

    H_a^j(p,q)       :=1/√(2π q)∫ a^j(u)e^-(p-u)^2
    /2qdu,   j=0,1,
    
    H_b(p,q)       :=1/√(2π q)∫ b(u)e^-(p-u)^2/2q
    du,  p∈ℝ, q>0.

Let moreover a^j and b be such that the partial derivatives,

    ∂_pH_a^j(p,q), ∂_qH_a^j(p,q)  
    j=0,1,  and   ∂_pH_b(p,q), ∂_qH_b
    (p,q), 
    extend continuously to any (p,q)∈ℝ× ℝ_≥0.

It is not difficult to see that (<ref>) holds if a^j and b are
entire functions for which the coefficients of their power series around u=0
decay fast enough to zero (which is trivially satisfied for any polynomial). A
complete characterization of a^j and b such that (<ref>) holds, is
connected with analytic vectors for semigroups related to the heat kernel and
considered beyond the scope of this paper however.


 Let a^j and b satisfy (<ref>). (i) Then
the following system of ODEs

    G_t^'     =H_b^2(  A_t,G_t)  +2H_a^1(
    A_t,G_t)  G_t
    
    A_t^'     =H_a^0(  A_t,G_t)  +H_a^1(
    A_t,G_t)  A_t,  (A_0,G_0)=(  x_0,0)
    ,

has for 0≤ t<t_∞≤∞, i.e. up to some possibly finite
exploding time t_∞, a unique solution (A_t,G_t)∈ℝ× ℝ_≥0. (ii) The McKean-Vlasov SDE

    dX_t=(E[  a^0(X_t)]  +X_t E[
    a^1(X_t)]  ) dt+E[  b(X_t)]  dW_t,  
    X_0=x_0

is then equivalent to

    dX_t=(  H_a^0(  A_t,G_t)  +H_a^1(
    A_t,G_t)  X_t)  dt+H_b(  A_t,G_t)
    dW_t,  X_0=x_0,

and has explicit solution,

    X_t     =x_0e^∫_0^tH_a^1(  A_s,G_s)  ds
    +∫_0^tH_a^0(  A_s,G_s)  e^∫_s^tH_a^1
    (  A_r,G_r)  drds
         +∫_0^tH_b(  A_s,G_s)  e^∫_s^tH_a^1
    (  A_r,G_r)  drdW_s,  0≤ t<t_∞
    ≤∞.

Note: the Wiener integral in (<ref>) can be interpreted by an ordinary
integral after partial integration, due to the smoothness of the
(deterministic) integrand.




 §.§ Affine structures


 Let us consider affine functions

    a^0(u)       =a_0^0+a_1^0u,
    
    a^1(u)       =a_0^1+a_1^1u,
    
    b(u)       =b_0+b_1u.

Then for c≡ a^0, c≡ a^1, and c≡ b, respectively, we
have

    H_c(p,q)       =1/√(2π q)∫ c(u)e^-(p-u)^2/2qdu
         =1/√(2π q)∫ c_0e^-(p-u)^2/2qdu+1/√(2π q)∫ c_1ue^-(p-u)^2/2qdu
         =c_0+c_1p

with c(u)=c_0+c_1u. In particular, the H_c(p,q) do not depend on
q, and so (<ref>) simplifies to

    A_t^'=a_0^0+(  a_1^0+a_0^1)  A_t+a_1
    ^1A_t^2,  A_0=x_0.

We first consider the case a_1^1=0, then (<ref>) reads
A_t^'=a_0^0+(  a_1^0+a_0^1)  A_t with
solution

    A_t     =(  x_0+a_0^0/a_1^0+a_0^1)
    e^(  a_1^0+a_0^1)  t-a_0^0/a_1^0+a_0
    ^1  if   a_1^0+a_0^1≠0,  and
    
    A_t     =x_0+a_0^0t  if   a_1^0+a_0^1=0.

For the case a_1^1≠0 the solution (checked by Mathematica) is as
follows. If D:=(  a_1^0+a_0^1)  ^2-4a_0^0a_1
^1<0, a_1^1≠0,

    A_t=-(  a_1^0+a_0^1)  /2a_1^1+√(-D)/2a_1^1tan[   1/2√(-D)t+arctan[
    a_1^0+a_0^1+2a_1^1x_0/√(-D)]  ]  .

If D>0, a_1^1≠0,

    A_t=1/2(  √(D)-a_1^0-a_0^1)  /a_1^1
    +x_0-1/2(  √(D)-a_1^0-a_0^1)
    /a_1^1/1+1/2(  √(D)+a_1^0+a_0^1+2a_1^1
    x_0)  (  e^ -√(D)t-1)  /√(D).

If D=0, a_1^1≠0,

    A_t=-a_1^0+a_0^1/2a_1^1+1/a_1^1a_1^0+a_0^1+2a_1^1x_0/2-(  a_1^0+a_0^1+2a_1
    ^1x_0)  t.

As a result, the McKean-Vlasov SDE

    dX_t=(a_0^0+a_1^0A_t+(  a_0^1+a_1^1A_t)
    X_t)dt+(  b_0+b_1A_t)  dW_t

has the following (unique) solution

    X_t     =x_0e^∫_0^t(  a_0^1+a_1^1A_s)
    ds+∫_0^t(  a_0^0+a_1^0A_s)  e^∫_s
    ^t(  a_0^1+a_1^1A_r)  drds
         +∫_0^t(  b_0+b_1A_s)  e^∫_s^t(
    a_0^1+a_1^1A_r)  drdW_s,

where A_t is given by (<ref>), (<ref>), (<ref>), or (<ref>).


By taking in Section <ref>

    a(x,u)=a_1^0u+a_0^1x,  b(x,u)=b_0,  a_1
    ^0+a_0^1<0,

we get essentially the Shimizu-Yamada model. From (<ref>) we then have

    A_t=x_0e^ (  a_1^0+a_0^1)  t,

and from (<ref>) we then get the explicit solution

    X_t=x_0e^(  a_1^0+a_0^1)  t+∫_0^tb_0
    e^a_0^1(t-s)dW_s

which is Gaussian with mean x_0e^(  a_1^0+a_0^1)  t
and variance b_0^2e^2a_0^1t-1/2a_0^1, and which is
consistent with the terminology in (<cit.>, Section 3.10),
where a_1^0+a_0^1=-γ and a_0^1=-γ-κ.



 By taking in Section <ref>

    a(x,u)=(  a_0^1+a_1^1u)  x,  b(x,u)=b_0,

we straightforwardly get from (<ref>),

    A_t=x_0e^ a_0^1t/1-a_1^1/a_0^1x_0(
    e^ a_0^1t-1)  ,

and

    X_t=x_0e^∫_0^t(  a_0^1+a_1^1A_s)  ds
    +∫_0^tb_0e^∫_s^t(  a_0^1+a_1^1A_r)
    drdW_s,

respectively. Plugging (<ref>) into (<ref>) then yields

    X_t=x_0e^ a_0^1t/1-a_1^1/a_0^1x_0(
    e^ a_0^1t-1)  +b_0e^ a_0^1t/1-a_1^1
    /a_0^1x_0(  e^ a_0^1t-1)  Γ_t

with Gaussian Γ_t=∫_0^t(  1-a_1^1/a_0^1
x_0(  e^ a_0^1s-1)  )  e^ -a_0^1sdW_s. In
particular, if a_0^1=0 we get

    A_t=x_0/1-a_1^1x_0t,

and

    X_t=x_0/1-a_1^1x_0t+b_0∫_0^t1-a_1^1
    x_0s/1-a_1^1x_0tdW_s.





 From Example <ref> it is clear that if a_1^1
≠0, the affine McKean-Vlasov solution may explode in finite time. This is
not surprising since in this case the derivative ∂_ua(x,u) is
unbounded and so the main results in <cit.> do not apply. On
the other hand, it is easy to check that for the case a_1^1≠0, the
affine solutions in Section <ref> are non-exploding whenever,

    D≥0 and √(D)≥ a_1^0+a_0^1+2a_1^1x_0.

That is, in the case D≥0, a_1^1≠0, it is always possible to
choose x_0 such that the solution does or does not explode.




 §.§ Kuramoto-Shinomoto-Sakaguchi type models


In the Kuramoto-Shinomoto-Sakaguchi model  the nonlinear one-dimensional
Fokker-Planck equation (<ref>) is considered in the domain (t,x)∈(
0,∞)  ×(  0,2π), where b=1, a(x,y)  =a(x-y)=-d/dxU_MF(x-y) with

    U_MF(z)  =-∑_n=1^∞c_ncos(nz)

and the process starts in u at time zero, for some fixed u∈(
0,2π), see for details <cit.> (Sect. 5.3.2).
Thus a is a 2π-periodic function related to a 2π-periodic potential.
Let us  consider the corresponding McKean-Vlasov SDE 

    {[                              X_t =u+∫_0^t∫_ℝa(X_s-y)μ_s(dy)ds+W_t;                              μ_t             =Law(X_t),   t∈0,T], ].

and define the integer valued function k(x):=max{  j∈ℤ:2π j≤
x}. Obviously, the process

    Y_t:=X_t-2π k(X_t)

has state space [  0,2π). Let ρ_t(x;u)=ρ_t(x) be
the density of Y_t, which is concentrated on (  0,2π).
Note that for any 2π-periodic function f we have by (<ref>) that

    ∫_0^2πf(x)ρ_t(x)dx=E[  f(Y_t)]  =E[
    f(X_t)]  =∫_-∞^∞f(x)μ_t(x)dx,

and for any test function g with support in (0,2π) it holds that,

    ∫_0^2πg(x)ρ_t(x)dx      =E[  g(Y_t)]  =E[
    g(X_t-2π k(X_t))]  
         =∑_j∈ℤ∫_2π j^2π(  j+1)  g(x-2π
    j)μ_t(x)dx
         =∫_0^2πg(z)∑_j∈ℤμ_t(z+2π j)dz,

that is

    ρ_t(z)=∑_j∈ℤμ_t(z+2π j)

for z∈ (0,2π).
Thus, in particular,

    ∫_ℝa(x-y)μ_t(y)dy=∫_0^2πa(x-y)ρ_t
    (y)dy.

an  (<ref>) is equivalent to

    X_t     =u+∫_0^t∫_0^2πa(X_s-y)ρ_s(y)dyds+W_t
    ρ_t     =Law(Y_t),   t∈0,T],

(see (<ref>)). Note that by using (<ref>) and (<ref>) it
straightforwardly follows that ρ_t(x)=ρ_t(x;u) satisfies
(<ref>) in the above context. Instead of taking the scalar product in
L_2(  ℝ^d,w), we now consider the scalar product in
L_2(  [  0,2π)  ), i.e. w≡1, and take for
(  φ_k) the standard (total) orthonormal trigonometric
basis consisting of the 2π-periodic functions (  2π)
^-1/2, π^-1/2cos(  my) and π^-1/2sin(
my)  , m=1,2,… suitably ordered. Thus, by defining

    γ_k(t)   =∫_0^2πρ_t(y)φ_k(  y)
    dy,   α_k(  x)     =∫_0^2πa(x-y)φ_k(
    y)  dy,

one has

    ρ_t(y)=∑_k=0^∞γ_k(t)φ_k(
    y)     and ∫_ℝa(x-y)μ_t(y)dy=∑
    _k=0^∞α_k(  x)  γ_k(t),

due to (<ref>). That is (<ref>) reads,

    X_t^i,K,N=u+∫_0^t∑_k=0^Kγ_k^N(s)α_k
    (X_s^i,K,N) ds+ W_t^i

with γ_k^N as in (<ref>). Next we may follow
(<ref>) for the corresponding Euler scheme. Finally, the
estimator for the density ρ_t reads

    ρ_t^K_test,K,N(y):=∑_k=1^K_test
    γ_k^N(t)φ_k(y)

(cf. the estimator for μ_t in Section 3.1).



§ NUMERICAL TEST CASES






 §.§ Affine MVSDE models


Let us now test the numerical performance of the projected particle approach
for the processes discussed in Section <ref>. Consider the situation
where

    a^0(u)=(1+u^M)exp(-u^2/2),   a^1(u)=ρexp(-u^2/2),  
    b(x,u)≡σ

for some M>0, ρ≥0 and σ>0. Then, by using the Hermite
functions (<ref>) with w≡1 and the well-known identity

    u^M=1/2^M∑_m=0^⌊M/2⌋
    M!/m!(M-2m)!H_M-2m(u),

we derive straightforwardly,

    α_k(x)       =∫(1+u^M)exp(-u^2)H_k(u) du+ρ
    x·∫exp(-u^2)H_k(u) du
         =
    
    0    if   k>M  or  k  is uneven 
    π^1/4/2^M-k/2M!/(  M-k/2)  !√(k!)   if   0≤ k≤ M, s.t.  M-k is even
    
         +(  1+ρ x)  ·
    π^1/4   if   k=0,
    
    0    if   k>0.

On the other hand, by some algebra we get

    H_a^0(p,q)       =1/√(2π q)∫(1+u^M)e^-u^2/2
    e^-(p-u)^2/2qdu=1/√(1+q)e^-p^2/2(1+q)
         +1/√(2π(1+q))e^-p^2/2(1+q)∫(  √(q/1+q)y+p/1+q)  ^Me^-y^2/2dy,

and

    H_a^1(p,q)=ρ/√(1+q)e^-p^2/2(1+q).

The explicit solution of the MVSDE

    dX_t=(  𝖤[  (1+X_t^M)exp(-X_t^2/2)]
    +ρ X_t𝖤[  exp(-X_t^2/2)]  )   dt+σ
    dW_t

is given by (<ref>) and (<ref>). Hence the density of X_t is
normal with mean

    x_0e^∫_0^tH_a^1(  A_s,G_s)  ds+∫_0
    ^tH_a^0(  A_s,G_s)  e^∫_s^tH_a^1(
    A_r,G_r)  dr ds

and variance

    σ^2∫_0^te^2∫_s^tH_a^1(  A_r,G_r)
    dr ds.

In our numerical example we take M=2, ρ=-1, σ=1 and x_0=0.
Our aim is to approximate the normal density of X_1 by using our projected
particle method based on Hermite basis. To this end, we first simulate N
paths of the process X̅^i,K,N, defined in (<ref>)
with a time step h=0.02. Since M=2, the case K=2 corresponds to a
perfect approximation of the integral ∫_ℝ^da(x,y)μ
_t(y) dy=∑_k=0^2α_k(x)γ_k(t). Next using the obtained
sample X̅_1^1,K,N,…,X̅_1^N,K,N, we construct
projection estimates for the density of X_1 by using Hermite basis
functions of order K_test∈{1,2,…,10}. The mean (0.727)
and the variance (0.487) of the true normal density are approximated by
solving the ODE system  (<ref>) using Euler method with time step
0.0001. The Figure <ref> shows the box plots of L_2-distance between μ_1 and μ_t^K_test,K,N for
K∈{1,2} based on 50 different replications of the process X̅^i,K,N. As can be seen, the choice of K_test is crucial
and depends on K and N. It also should be stressed that the truncation
error dominates the statistical one already for medium sample sizes. An
optimal balance between K, N and K_test can be found by
analyzing the right hand side of (<ref>) under various
assumptions on the coefficients (A_k,α), (A_k,β) and
(γ_k). 



 §.§ Convolution-type MVSDE models


Consider the MVSDE of the form:

    dX_t=_X^'[  Q(X_t-X_t^')]   dt+σ dW_t,   t∈0,1],   X_0∼𝒩(0,1),

i.e. of the form (<ref>) with a(x,y)=Q(x-y), b(x,y)=σ and
μ_0(x)=(1/√(2π))e^-x^2/2. Let us again use the Hermite basis
to approximate the density of X_t for any t∈0,1]. In the case
Q(x)=e^-x^2/2, we explicitly derive via repeated integration by parts

    ∫_ℝe^-(x-y)^2/2-x^2/2H_n(x) dx       =e^-y^2/4
    /2∫ e^-(z-y)^2/4H_n(z/2) dz
         =√(π) e^-y^2/4/2(  1/2)
    ^n-1(2y)^n.

As a result

    α_n(y)=∫ e^-(x-y)^2/2-x^2/2H_n(x) dx=π
    ^1/4(  1/2)  ^n/2y^n/√(n!)e^-y^2/4,

where H_n stands for the normalized Hermite polynomial of order
n. We take σ=0.1. Using the Euler scheme (<ref>) with time
step h=1/L=0.01, we first simulate N=500 paths of the time discretized
process X̅^·,N. Next, by means of the closed form expressions for
α_n, we generate N paths of the projected approximating process
X̅^·,K,N, K∈{1,…,20} (see (<ref>
)), using the same Wiener increments as for X̅^·,N, so that the
approximations X̅^·,N and X̅^·,K,N are coupled.
Finally, we compute the strong approximation error

    E_N,K=√(1/N∑_i=1^N(X̅_1^i,K,N-X̅
    _1^i,N)^2)

of the projective system relative to the system (<ref>) and record
times needed to compute approximations X̅_1^·,N and X̅_1^·,K,N, respectively. Figure <ref> shows the
(natural) logarithm of E_N,K versus the logarithm of the corresponding
(relative) computational time gain defined as (comp. time due to
(<ref>) - comp. time due to (<ref>))/comp. time
due to (<ref>), for values K∈{1,…,20}. As can be seen,
the relation between logarithmic strong error and logarithmic computational
time gain can be well approximated by a linear function. On the right-hand
side of Figure <ref> we depict the projection estimate for the
density of X_1 corresponding to K=20. Note that we compare two particle
systems (projected and non projected ones) for a fixed N and are mainly
interested in the dependence of their strong distance on K. In fact, the
choice of N doesn't have much influence on E_N,K, provided N is large
enough.





§ PROOFS






 §.§ Proof of Theorem <ref>


Let us introduce

    𝔞_K,N(x,y):=1/N∑_j=1^N∑_k=1^Kα
    _k(x)φ_k(y^j)=1/N∑_j=1^N∑_k=1^Kφ
    _k(y^j)∫ a(x,u)φ_k(u)w(u)du,


    𝔟_K,N(x,y):=1/N∑_j=1^N∑_k=1^Kβ
    _k(x)φ_k(y^j)=1/N∑_j=1^N∑_k=1^Kφ
    _k(y^j)∫ b(x,u)φ_k(u)w(u)du,

and

    𝔞_s(x)       :=∫_ℝ^da(x,u)μ_s(du)ds,
    𝔟_s(x)       :=∫_ℝ^db(x,u)μ_s(du)ds

for any x∈ℝ^d,y∈ℝ^d× N. We so have
that

    Δ_t^i:=X_t^i,K,N-X_t^i     =∫_0^t(
    𝔞_K,N(X_s^i,K,N,X_s^K,N) -𝔞_s(X_s
    ^i))  ds
         +∫_0^t(  𝔟_K,N(X_s^i,K,N,X_s
    ^K,N) -𝔟_s(X_s^i))  dW_s^i,

where W^i, i=1,...,N, are i.i.d. copies of the m-dimensional Wiener
process W. Hence,

    |Δ_t^i| ^p   ≤2^p-1t^p-1∫_0
    ^t|𝔞_K,N(X_s^i,K,N,X_s^K,N
    ) -𝔞_s(X_s^i)| ^pds
         +2^p-1d^p-1∑_q=1^d|∫_0^t(  𝔟
    _K,N^q(X_s^i,K,N,X_s^K,N) -𝔟_s^q(X_s
    ^i))  dW_s^i| ^p,

where

    𝔟_K,N^q:=(  𝔟_K,N^q,1,...,𝔟
    _K,N^q,m)  ,  q=1,...,d,

denote the rows of the ℝ^d× m valued 𝔟_K,N,
and so we have with

    Δ_t^p:=1/N∑_i=1^Nsup_s∈
    0,t]|Δ_s^i| ^p

the bound

    Δ_t^p   ≤2^p-1t^p-11/N∑_i=1^N
    ∫_0^t|𝔞_K,N(X_s^i,K,N,X_s
    ^K,N) -𝔞_s(X_s^i)| ^pds
         +2^p-1d^p-1∑_q=1^d1/N∑_i=1^Nsup_s∈0,t]|∫_0^s(  𝔟_K,N^q(X
    _s^i,K,N,X_s^K,N) -𝔟_s^q(X_s^i))  dW_s
    ^i| ^p
         =:2^p-1t^p-1 Term_1+2^p-1d^p-1 Term_2.

Assumption (AC) implies

    |𝔞_K,N(x,y)-𝔞_K,N(x^',y^')|     =|1/N∑_j=1^N∑_k=1^K(
    α_k(x)φ_k(y_j)-α_k(x^')φ_k
    (y_j^'))  |
       ≤1/N∑_j=1^N∑_k=1^K|α_k
    (x)-α_k(x^')||φ_k(y_j^')|
         +1/N∑_j=1^N∑_k=1^K|α_k(x)||φ_k(y_j)-φ_k(y_j^')|
       ≤| x-x^'| D_φB_α+L_φA_α/N(1+| x| )∑_j=1^N|
    y_j-y_j^'| .

Hence

    |𝔞_K,N(x,y)-𝔞_K,N(x^',y^')| ^p   ≤2^p-1| x-x^'|
    ^pD_φ^pB_α^p
         +2^p-1L_φ^pA_α^p(1+| x| )^p
    1/N∑_j=1^N| y_j-y_j^'| ^p.

So it holds that

    |𝔞_K,N(X_s^i,X_s)-𝔞_K,N(X
    _s^i,K,N,X_s^K,N) | ^p   ≤2^p-1D_φ
    ^pB_α^p|Δ_s^i|^p
         +2^p-1L_φ^pA_α^p(1+| X_s^i|
    )^p1/N∑_j=1^N|Δ_s^j|^p,

and then it follows that, with regard to Term_1,

    E[  Term_1]      ≤2^2p-2D_φ
    ^pB_α^p∫_0^tE[  Δ_s^p
    ]  ds
         +2^2p-2L_φ^pA_α^p∫_0^tE[
    Δ_s^p·1/N∑_i=1^N(1+| X_s
    ^i| )^p]  ds
         +2^p-11/N∑_i=1^N∫_0^tE[  |𝔞_K,N(X_s^i,X_s) -𝔞_s(X_s
    ^i)| ^p]  ds.

Let us now consider the middle term. Set

    ζ_s,N:=1/N∑_i=1^N(1+| X_s^i|
    )^p-1/N∑_i=1^NE[  (1+| X_s
    ^i| )^p]

so that

    E[  Δ_s^p·1/N∑_i=1
    ^N(1+| X_s^i| )^p]        =1/N
    ∑_i=1^NE[  (1+| X_s^i|
    )^p]  ·E[  Δ_s^p]
    
         +E[  ζ_s,N·Δ_s^p]  .

For arbitrary but fixed θ>0, it holds that

    E[  ζ_s,N·Δ_s^p]
    =E[  ζ_s,N·Δ_s^p 1_{ζ_s,N≤θ}]  +E[  ζ_s,N
    ·Δ_s^p 1_{ζ_s,N>θ}]  ,

where on the one hand

    E[  ζ_s,N·Δ_s^p 1_{ζ_s,N≤θ}]  ≤θE[
    Δ_s^p]

and on the other

    E[  ζ_s,N·Δ_s^p 1_{ζ_s,N>θ}]  ≤√(E[  ζ
    _s,N^21_{ζ_s,N>θ}]  )√(E[  (  Δ_s^p)  ^2]  ).

Due to (<ref>) we have that for any η>0, there exists
C_θ,η>0 such that

    E[  ζ_s,N^21_{ζ_s,N>θ}]  =1/NE[  (  √(N)ζ_s,N)
    ^21_{√(N)ζ_s,N>θ√(N)}]
    ≤C_θ,η^2/N^η+1,   0≤ s≤ T,

for N large enough and

    E[  (  Δ_s^p)  ^2]      ≤E[  1/N∑_j=1^Nsup_r∈
    0,T]|Δ_r^j| ^2p]  =E[
    sup_r∈0,T]|Δ_r^·| ^2p]
    
         =E[  sup_r∈0,T]| X_r^·,K,N
    -X_r^·| ^2p] 
       ≤2^2p-1E[  sup_r∈0,T]| X_r
    ^·,K,N| ^2p]  +2^2p-1E[  sup
    _r∈0,T]| X_r^·| ^2p] 
       ≤ D_1+D_2=:D^2,

where due to Theorem <ref>, Appendix <ref>,

    2^2p-1E[  sup_r∈0,T]| X_r^·
    ,K,N| ^2p]  ≤ D_1<∞  uniform in N
    and K,

and

    D_2:=2^2p-1E[  sup_r∈0,T]| X_r
    ^·| ^2p]  <∞

due to (<ref>). Thus, by combining (<ref>)–(<ref>), one has

    E[  Δ_s^p·1/N∑_i=1
    ^N(1+| X_s^i| )^p]  ≤ F_1^p
    ·E[  Δ_s^p]  +F_2
    /N^p/2+1/2

with F_1:=θ^1/p+sup_0≤ s≤ T‖ 1+|
X_s|‖ _p and F_2:=C_θ,pD, where we have
taken η=p. Set now

    H(s):=E[  Δ_s^p]  ,

then the estimate (<ref>) (cf. (<ref>)) reads

    1/N∑_i=1^N∫_0^tE[  |𝔞_K,N(X_s^i,K,N,X_s^K,N) -𝔞_s(X_s
    ^i)| ^p]  ds
    ≤(  2^2p-2D_φ^pB_α^p+2^2p-2L_φ
    ^pA_α^pF_1^p)  ∫_0^tH(s)ds+2^2p-2L_φ
    ^pA_α^pF_2/N^p/2+1/2t
    
    +2^p-11/N∑_i=1^N∫_0^tE[  |𝔞_K,N(X_s^i,X_s) -𝔞_s(X_s
    ^i)| ^p]  ds.

Regarding the term Term_2 we call upon the Burkholder-Davis-Gundy's
inequality which states that for any p≥1,

    ‖sup_s∈0,t]|∫_0^s(
    𝔟_K,N^q(X_s^i,K,N,X_s^K,N) -𝔟_s
    ^q(X_s^i))  dW_s^i|‖ _p
       ≤ C_p(  E[  (  ∫_0^t|(
    𝔟_K,N^q(X_s^i,K,N,X_s^K,N) -𝔟_s
    ^q(X_s^i))  | ^2ds)  ^p/2]  )
    ^1/p.

This implies that for p≥2,

    Esup_s∈0,t]|∫_0^s(
    𝔟_K,N^q(X_s^i,K,N,X_s^K,N) -𝔟_s
    ^q(X_s^i))  dW_s^i| ^p
       ≤ C_p^pE[  (  ∫_0^t|(
    𝔟_K,N^q(X_s^i,K,N,X_s^K,N) -𝔟_s
    ^q(X_s^i))  | ^2ds)  ^p/2] 
       ≤ C_p^pt^p/2-1E[  ∫_0^t|(
    𝔟_K,N^q(X_s^i,K,N,X_s^K,N) -𝔟_s
    ^q(X_s^i))  | ^pds] 
       ≤ C_p^pt^p/2-1E[  ∫_0^t|(
    𝔟_K,N(X_s^i,K,N,X_s^K,N) -𝔟_s(X_s
    ^i))  | ^pds]  .

Now, completely analogue to the derivation of (<ref>), we get

    1/N∑_i=1^N∫_0^tE[  |𝔟_K,N(X_s^i,K,N,X_s^K,N) -𝔟_s(X_s
    ^i)| ^p]  ds
    ≤(  2^2p-2D_φ^pB_β^p+2^2p-2L_φ
    ^pA_β^pF_1)  ∫_0^tH(s)ds+2^2p-2L_φ
    ^pA_β^pF_2/N^p/2+1/2t
    
    +2^p-11/N∑_i=1^N∫_0^tE[  |𝔟_K,N(X_s^i,X_s) -𝔟_s(X_s
    ^i)| ^p]  ds.

Now by taking expectations on both sides of (<ref>) and gathering all
together, we arrive at

    H(t)≤(  D_φ^pB_α^pT^p-1+L_φ^pA_α^pF_1^pT^p-1. 
    .  +C_p^pD_φ^pB_β^pd^pT^p/2-1+C_p
    ^pL_φ^pA_β^pd^pT^p/2-1F_1^p)  2^3p-3
    ∫_0^tH(s)ds
    
    +2^3p-3(  L_φ^pA_α^pT^p+d^pC_p^pL_φ^pA_β^pT^p/2)  F_2/N^p/2+1/2
    
    +2^2p-2T^p-11/N∑_i=1^N∫_0^tE[
    |𝔞_K,N(X_s^i,X_s) -𝔞_s
    (X_s^i)| ^p]  ds
    
    +2^2p-2d^pC_p^pT^p/2-11/N∑_i=1^N∫_0
    ^tE[  |𝔟_K,N(X_s^i
    ,X_s) -𝔟_s(X_s^i)| ^p]  ds.

We next proceed with explicit estimates for the last two terms above. Let us
write

    𝔞_K,N(X_s^i,X_s)-𝔞_s(X_s^i)=∑_k=1
    ^Kα_k(X_s^i)∑_j=1^N1/N(  φ_k
    (X_s^j)-γ_k(s))  -∑_k=K+1^∞α_k(X_s
    ^i)γ_k(s),

then we have by the Minkowski inequality,

    ‖𝔞_K,N(X_s^i,X_s)-𝔞_s(X_s
    ^i)‖ _p   ≤∑_k=1^K‖α_k(X_s
    ^i)1/N∑_j=1^Nξ_k^j‖ _p
         +∑_k=K+1^∞‖α_k(X_s^i)γ_k
    (s)‖ _p,

where ξ_k^j:=φ_k(X_s^j)-γ_k(s), j=1,…,N,
have mean zero. Let us now observe that

    E[  .  |∑_j=1^Nξ_k^j|
    ^p| X^i]        =E[  .  |ξ_k^i+∑_j≠ i^Nξ_k^j| ^p|
    X^i] 
       ≤2^p-1E[  .  |ξ_k^i|
    ^p+|∑_j≠ i^Nξ_k^j| ^p|
    X^i] 
       ≤2^2p-1D_k,φ^p+2^p-1E[  |∑_j≠ i^Nξ_k^j| ^p]

using (<ref>). For p≥2, it follows from the Rosenthal's inequality
that,

    E[  |∑_j≠ i^Nξ_k^j|
    ^p]  ≤ C_p^(1)(  (  ∑_j≠ i^NE
    |ξ_k^j| ^2)  ^p/2+∑_j≠ i
    ^NE|ξ_k^j| ^p)

for a constant C_p^(1) only depending on p, and, in fact, for p=2 we
have simply,

    E[  |∑_j≠ i^Nξ_k^j|
    ^p]  =∑_j≠ i^NE|ξ_k^j|
    ^2.

Thus, for p≥2,

    E[  .  |1/N∑_j=1^Nξ_k
    ^j| ^p| X_s^i]      ≤2^2p-1
    D_k,φ^p/N^p+2^p-1C_p^(1)/N^p(  (
    ∑_j≠ i^NE|ξ_k^j| ^2)
    ^p/2+∑_j≠ i^NE|ξ_k^j|
    ^p) 
       ≤2^2p-1D_k,φ^p/N^p+2^2p-1C_p
    ^(1)D_k,φ^p/N^p/2+2^2p-1C_p^(1)D_k,φ^p
    /N^p-1
       ≤(  C_p^(2))  ^pD_k,φ^p/N^p/2  for N>N_p and some constants  C_p^(2)>0, N_p>0.

So for any p≥2,

    ‖α_k(X_s^i)1/N∑_j=1^Nξ_k
    ^j‖ _p^p   ≤ A_k,α^pE[  (
    1+| X_s^i|)  ^pE[  .
    |1/N∑_j=1^Nξ_k^j| ^p|
    X_s^i]  ] 
       ≤ A_k,α^pD_k,φ^p(  C_p^(2))
    ^p/N^p/2E[  (  1+| X_s|)
    ^p]  ,

hence

    ‖α_k(X_s^i)1/N∑_j=1^Nξ_k
    ^j‖ _p≤ C_p^(2)A_k,αD_k,φF_3
    N^-1/2  with   F_3:=sup_0≤ s≤ T‖
    1+| X_s|‖ _p,

and further

    ∑_k=K+1^∞‖α_k(X_s^i)γ_k(s)‖
    _p≤ F_3∑_k=K+1^∞A_k,α|γ
    _k(s)| .

We thus obtain,

    ‖𝔞_K,N(X_s^i,X_s)-𝔞_s(X_s
    ^i)‖ _p≤ C_p^(2)A_αD_φF_3N^-1/2
    +F_3∑_k=K+1^∞A_k,α|γ_k(s)| ,

that is,

    E[  |𝔞_K,N(X_s^i,X_s)-𝔞
    _s(X_s^i)| ^p]      ≤2^p-1(  C_p
    ^(2))  ^pA_α^pD_φ^pF_3^pN^-p/2
         +2^p-1F_3^p(  ∑_k=K+1^∞A_k,α|γ_k(s)|)  ^p.

Analogously we get

    E[  |𝔟_K,N(X_s^i,X_s)-𝔟
    _s(X_s^i)| ^p]      ≤2^p-1(  C_p
    ^(2))  ^pA_β^pD_φ^pF_3^pN^-p/2
         +2^p-1F_3^p(  ∑_k=K+1^∞A_k,β|γ_k(s)|)  ^p.

Now, combining the estimates (<ref>) and (<ref>) with (<ref>),
yields for 0≤ t≤ T,


    H(t)     ≤(  C_p,φ,XT^p-1+D_p,φ,Xd^p
    T^p/2-1)  ∫_0^tH(s)ds
         +(  E_p,φ,XT^p+F_p,φ,Xd^pT^p/2+O(N^-1/2
    ))  N^-p/2
         +G_p,φ,XT^p-1∫_0^T(  ∑_k=K+1^∞
    A_k,α|γ_k(s)|)  ^pds
         +H_p,φ,Xd^pT^p/2-1∫_0^T(  ∑_k=K+1^∞A_k,β|γ_k(s)|)  ^pds

with abbreviations

    C_p,φ,X     =2^3p-3D_φ^pB_α^p+2^3p-3L_φ^pA_α^pF_1^p
    
    D_p,φ,X     =2^3p-3C_p^pD_φ^pB_β^p
    +2^3p-3C_p^pL_φ^pA_β^pF_1^p
    
    E_p,φ,X     =2^3p-3(  C_p^(2))  ^pA_α
    ^pD_φ^pF_3^p
    
    F_p,φ,X     =2^3p-3C_p^p(  C_p^(2))  ^p
    A_β^pD_φ^pF_3^p
    
    G_p,φ,X     =2^3p-3F_3^p
    
    H_p,φ,X     =2^3p-3C_p^pF_3^p.

Finally, the statement of the theorem follows from Gronwall's lemma by raising
the resulting inequality to the power 1/p, then using that (
∑_i=1^q|a_i|^p)  ^1/p≤∑_i=1^q|a_i| for
arbitrary a_i∈ℝ, p,q∈ℕ, a Minkowski type
inequality, and the observation that

    E[  Δ_T^p]  =1/N∑_i=1
    ^NE[  sup_s∈0,T]|Δ_s
    ^i| ^p]  =E[  sup_s∈
    0,T]|Δ_s^·| ^p]  .





 §.§ Proof of Theorem <ref>


(i): Under the assumption (<ref>) the functions H_a^j and H_b
are locally Lipschitz in ℝ× ℝ_≥0 and, obviously, their
extensions (p,q)→ H_a^j(p,|q|), H_b(p,|q|) are locally
Lipschitz in ℝ× ℝ. Thus, by standard ODE theory, there exists
a unique solution to the system

    G_t^'     =H_b^2(  A_t,| G_t|)  +2H_a^1(  A_t,| G_t|)  G_t
    
    A_t^'     =H_a^0(  A_t,| G_t|)  +H_a^1(  A_t,| G_t|)
    A_t,  (A_0,G_0)=(  x_0,0)  ,  0≤
    t<t_∞≤∞,

for some possibly finite explosion time t_∞. Then it can be
straightforwardly checked that this (unique) solution can be represented as

    G_t     =∫_0^tH_b^2(  A_s,| G_s|)  e^2∫_s^tH_a^1(  A_r,| G_r|)  drds
    
    A_t     =e^∫_0^tH_a^1(  A_s,| G_s|)  dsx_0+∫_0^tH_a^0(  A_s,| G_s
    |)  e^∫_s^tH_a^1(  A_r,|
    G_r|)  drds, 0≤ t<t_∞,

whence in particular G_t≥0 for 0≤ t<t_∞. This proves (i).

(ii): By straightforward differentiating with respect to t, it follows that
(<ref>) is a solution to (<ref>). Let us abbreviate in (<ref>)

    𝔞_t^0≡ H_a^0(  A_t,G_t)  ,  𝔞_t^1≡ H_a^1(  A_t,G_t)  ,  𝔟_t≡ H_b(  A_t,G_t)  ,  0≤ t<t_∞.

The characteristic function of X_t in (<ref>) then takes the form

    φ_t(v)=exp[  𝔦v∫_0^t𝔞_s
    ^0e^∫_s^t𝔞_r^1drds-1/2v^2∫_0
    ^t𝔟_s^2e^2∫_s^t𝔞_r^1drds+𝔦
    ve^∫_0^t𝔞_s^1dsx_0]  .

Since

    e^-(p-u)^2/2q/√(2π q)=1/2π∫
    e^-𝔦vuexp[  𝔦vp-v^2q/2]  dv,

we have for j=0,1,

    H_a^j(p,q)=1/2π∫ a^j(u)du∫exp[  𝔦
    vp-v^2q/2]  e^-𝔦vudv.

It then follows that

    H_a^j(e^∫_0^t𝔞_s^1dsx_0+∫_0
    ^t𝔞_s^0e^∫_s^t𝔞_r^1drds,∫_0
    ^t(  𝔟_s^0)  ^2e^2∫_s^t𝔞
    _r^1drds)
         =1/2π∫ a^j(u)du∫φ_t(v)e^-𝔦
    vudv
         =∫ a^j(u)μ_t(u)du=E[  a^j(X_t)]  ,  j=0,1,

with μ_t being the density of X_t, and similarly,

    H_b(e^∫_0^t𝔞_s^1dsx_0+∫_0^t𝔞
    _s^0e^∫_s^t𝔞_r^1drds,∫_0^t(
    𝔟_s^0)  ^2e^2∫_s^t𝔞_r^1
    drds)=E[  b(X_t)]  .

On the other hand, in view of (<ref>) and the fact that G≥0, one has

    ∫_0^t(  𝔟_s^0)  ^2e^2∫_s
    ^t𝔞_r^1drds       =G_t
    
    e^∫_0^t𝔞_s^1dsx_0+∫_0^t𝔞_s
    ^0e^∫_s^t𝔞_r^1drds       =A_t,

that is, by (<ref>), (<ref>), and (<ref>), we obtain (<ref>) from
(<ref>).



§ APPENDIX






 §.§ Existence of moments





 Fix some p≥2 and suppose that E[|X_0
|^p]<∞. Then it holds under assumptions (AC) and (AF),

    ‖sup_s∈0,T]| X_s^·,K,N|‖ _p<∞,

uniformly in K and N.



Fix some i∈{1,…,N} and for every R>0 introduce the stopping time

    τ_i,R=inf{  t∈0,T] :| X_t^i,K,N-X_0
    ^i| >R}  .

We obviously have

    sup_t∈0,T]| X_t∧τ_i,R^i,K,N|≤
    R+| X_0^i|

so that the non-decreasing function f_R(t):=‖sup_s∈0,t]| X_s∧τ_i,R^i,K,N|‖
_p, t∈0,T], is bounded by R+‖ X_0^i‖
_p. On the other hand

    sup_s∈0,t]| X_s∧τ_i,R^i,K,N|   ≤| X_0^i| +∫_0^t∧τ_i,R|𝔞_K,N(X_r^i,K,N,X_r^K,N)| dr
         +sup_s∈0,t]|∫_0^s∧τ_i,R
    𝔟_K,N(X_r^i,K,N,X_r^K,N) dW_r^i|
       ≤| X_0^i| +∫_0^t∧τ_i,R|𝔞_K,N(X_r^i,K,N,X_r^K,N)| dr
         +∑_q=1^dsup_s∈0,t]|∫_0^s∧τ
    _i,R𝔟_K,N^q(X_r^i,K,N,X_r^K,N) dW_r
    ^i|

(cf. (<ref>)). It then follows from the Minkowski and BDG inequality
that

    f_R(t)     ≤‖ X_0‖ _p+∫_0^t‖
    1_{s≤τ_i,R}𝔞_K,N(X_s^i,K,N,X_s^K,N
    )‖ _p ds
         +dC_p^BDG‖√(∫_0^t∧τ_i,R|𝔟_K,N(X_s^i,K,N,X_s^K,N)| ^2 ds)‖
    _p
       ≤‖ X_0‖ _p+A_αD_φ∫_0
    ^t‖(  1+| X_s∧τ_i,R^i,K,N|)  ‖ _p ds
         +A_βD_φdC_p^BDG‖√(∫_0^t|(  1+| X_s∧τ_i,R^i,K,N|)
    | ^2ds)‖ _p
       ≤‖ X_0‖ _p+A_αD_φ∫_0
    ^t(  1+‖| X_s∧τ_i,R^i,K,N|‖ _p)   ds
         +A_βD_φdC_p^BDG(  √(t)+(  ∫_0
    ^t‖| X_s∧τ_i,R^i,K,N|
    ^2‖ _p/2 ds)  ^1/2)

again by the Minkowski inequality (p≥2). Consequently, the function
f_R satisfies

    f_R(t)≤‖ X_0‖ _p+A_αD_φ∫_0
    ^t(  1+f_R(s))   ds+A_βD_φdC_p^BDG(
    √(t)+(  ∫_0^tf_R^2(s) ds)  ^1/2)  ,

that is,

    f_R(t)     ≤‖ X_0‖ _p+A_αD_φt+A_βD_φdC_p^BDG√(t)
         +A_αD_φ∫_0^tf_R(s) ds+A_βD_φ
    dC_p^BDG(  ∫_0^tf_R^2(s) ds)  ^1/2.

By Lemma <ref> (see Appendix) it follows that

    ‖sup_s∈0,T]| X_s∧τ_i,R
    ^i,K,N|‖ _p   ≤2e^(  2A_α
    D_φ+A_β^2D_φ^2d^2(  C_p^BDG)
    ^2)  T×
       (  ‖ X_0‖ _p+A_αD_φT+A_βD_φdC_p^BDG√(T))  .



Now note that the stopping times τ_i,R are non-decreasing in R, and
thus converges non-decreasingly to τ_i,∞ say, with τ
_i,∞∈0,T]∪{∞}. Thus,

    R→sup_s∈0,T]| X_s∧τ_i,R
    ^i,K,N|

is nondecreasing with

    lim_R↑∞sup_s∈0,T]| X_s∧τ_i,R
    ^i,K,N| ={c]c
    sup_s∈0,T]| X_s^i,K,N|  on {τ_i,∞=∞}
    ∞  on {τ_i,∞≤ T}
      .  .

Indeed, on the set {τ_i,∞≤ T} we have for any R>0,
| X_τ_i,R^i,K,N-X_0^i|≥ R with
τ_i,R≤ T, so that

    sup_s∈0,T]| X_s∧τ_i,R^i,K,N|≥| X_τ_i,R^i,K,N|≥| X_τ_i,R
    ^i,K,N|≥ R-| X_0^i| .

The Fatou lemma (<ref>) implies (with 0:=∞·0),

    ‖lim_R↑∞1_{τ_i,∞≤ T}sup_s∈0,T]| X_s∧τ_i,R^i,K,N|‖ _p     =∞· P(  {τ_i,∞≤
    T}) 
       ≤lim inf_R‖ 1_{τ_i,∞≤ T}
    sup_s∈0,T]| X_s∧τ_i,R^i,K,N|‖ _p
       ≤lim inf_R‖sup_s∈0,T]| X_s∧τ_i,R^i,K,N|‖ _p<∞,

because of (<ref>). So P(  {τ_i,∞≤ T})  =0, i.e. τ_i,∞=∞ almost surely. Again by the Fatou
lemma, (<ref>) then implies

    ‖sup_s∈0,T]| X_s^i,K,N|‖ _p≤lim inf_R‖sup_s∈0,T]|
    X_s∧τ_i,R^i,K,N|‖ _p<∞,

uniformly in K and N, because of (<ref>) again.


The following lemma is consequence of Gronwall's theorem.


 Let f: [0,T]→ℝ_+ and ψ:
[0,T]→ℝ_+ be two non-negative non-decreasing functions
satisfying

    f(t)≤ A∫_0^tf(s) ds+B(  ∫_0^tf^2(s) ds)
    ^1/2+ψ(t),   t∈0,T],

where A,B are two positive real constants. Then

    f(t)≤2e^(  2A+B^2)  t ψ(t),   t∈0,T].





It follows from the elementary inequality √(xy)≤1/2(
x/B+By)  , x,y≥0,B>0, that

    (  ∫_0^tf^2(s) ds)  ^1/2≤(  f(t)∫_0
    ^tf(s) ds)  ^1/2≤f(t)/2B+B/2∫_0
    ^tf(s) ds.

Plugging this into (<ref>) yields

    f(t)≤(2A+B^2)∫_0^tf(s) ds+2ψ(t).

Now the standard Gronwall inequality yields the desired result.


plain


