# A Control Variate Method Driven by Diffusion Approximation

Josselin Garnier, Laurent Mertz

► To cite this version:

Josselin Garnier, Laurent Mertz. A Control Variate Method Driven by Diffusion Approximation. Communications on Pure and Applied Mathematics, 2022, 75 (3), pp.455-492. 10.1002/cpa.21976 . hal-03897282

HAL Id: hal-03897282

https://inria.hal.science/hal-03897282v1

Submitted on 31 Oct 2025

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A control variate method driven by diffusion approximation

Josselin Garnier and Laurent Mertz

August 10, 2020

## Abstract

In this paper we examine a control variate estimator for a quantity that can be expressed as the expectation of a function of a random process, that is itself the solution of a differential equation driven by fast mean-reverting ergodic forces. The control variate is the same function for the limit diffusion process that approximates the original process when the mean-reversion time goes to zero. To get an efficient control variate estimator, we propose a coupling method for the original process and the limit diffusion process. We show that the correlation between the two processes indeed goes to one when the mean reversion time goes to zero and we quantify the convergence rate, which makes it possible to characterize the variance reduction of the proposed control variate method. The efficiency of the method is illustrated on a few examples.

## 1 Introduction

In this paper we consider a system driven by external time-dependent random forces and we aim to compute a quantity of interest that is the expectation of a function of the system. The system state is the solution of an ordinary differential equation (or a system of ordinary differential equations) driven by external forces which are modeled as stationary random processes. The driving processes may have complicated spectra that have to be taken into account to compute the quantity of interest. This happens for instance in seismic probabilistic risk assessment studies or in the analysis of the structural performance of installations under seismic excitations [31] or under other loading sources such as wind or waves [17,33]. For instance, the reliability of complex systems such as fixed or floating offshore wind turbines depends on its resistance against fatigue damage. Fatigue damage can be assessed by time-domain simulations in which the structure is subjected to wind, wave and current loads [5]. The different loads can be described by (locally) stationary Gaussian processes with tabulated power spectral densities (such as the JONSWAP spectrum [15]). We may then wish to estimate the mean cumulative fatigue damage or a probability of failure which corresponds to the exceedance of a threshold value.

Monte Carlo simulations are standard to estimate the quantities of interest but they may be very time consuming. We look for an efficient variance reduction technique in this framework. It is known from the diffusion approximation theory [6, 8, 21] that the driving forces can often be approximated by white noises and the responses of the system can then be modeled by stochastic differential equations. This makes it possible to implement a partial differential equation approach to compute the quantity of interest. However, the bias due to the approximation of the original driving force by a white noise may be significant and difficult to assess. To compensate for this bias, one may think at a control variate method [12]. Such a strategy has already been implemented in a Markov chain Monte Carlo context, where the goal was to sample from a complex invariant
probability distribution of a Markov chain for which an approximate distribution has a known expression. The expectation of the approximate distribution then provides an initial guess, which can be corrected by simulating the two coupled processes to estimate the difference (in expected values) between the true distribution and the approximate distribution [13]. The implementation of a control variate method in our framework requires to be able to simulate the system driven by the original driving force with its complicated spectrum and the limit system driven by the white noise in such a way that both systems are strongly correlated. Unfortunately, most diffusion approximation results are established in a weak sense [6,8]. Some strong results have been obtained but only when the drift is a term of order one [11, 19, 25], not when it is a zero-mean large term as we deal with in this paper. In this paper we build an efficient coupling between the original and limit systems, we establish a strong convergence result by quantifying the mean square distance between the original and limit processes, and we characterize the variance reduction of the control variate method. We show by our theoretical results and numerical simulations that the variance reduction can be dramatic.

Our method is relevant when the quality of the approximation of the driving forces by a white noise is moderate. If it is very accurate, then the quantity of interest can be estimated (up to a very small and negligible bias) by resolution of a Kolmogorov equation based on the limit diffusion system (or by a brute force Monte Carlo method applied to the limit system), so there is no need to apply a control variate method. If it is very poor, then the limit diffusion system is not correlated to the original system and the control variate method is not efficient. If it is moderate, then the bias of the estimation method that consists in replacing the original system by the limit one is non-negligible, and the two systems are correlated, so that the control variate method turns out to be very efficient.

The paper is organized as follows. In Section 2 we introduce the random ordinary differential equations addressed in this paper and we state the main results of the paper. Note that, motivated by applications in engineering mechanics and physics such as the study of the risk analysis of failure for mechanical structures subjected to random vibrations [1, 2, 7, 23] or the modeling of the stochastic dynamics of fluid-structure interaction in turbulent thermal convection [16], we also consider the case of multivalued ordinary differential equations. Sections 3-4 consider random ordinary differential equations. In Section 3 we state the diffusion approximation theorem that gives the convergence in probability of the original process to the limit process. In Section 4 we apply the control variate method to a few examples. The results are extended to the multivalued case in Sections 5-7. In particular Section 7 report numerical results for relevant engineering mechanics problems. The concluding remark of Section 8 connects our findings to the multilevel Monte Carlo literature.

## 2 Main results

We consider the $\mathbb{R}^n$-valued process $X^\epsilon = (X_t^\epsilon)_{t \in [0,T]}$ solution of the ordinary differential equation (ODE)<fnref n="1" />

$$ \frac{dX^\epsilon}{dt} = b(X^\epsilon) + \frac{1}{\epsilon}\sigma(X^\epsilon)\eta^\epsilon, \quad X_0^\epsilon = x_0, \qquad (1) $$

<footnote n="1">¹Throughout the paper, symbols of scalar quantities are printed in italic type, symbols of vectors are printed in bold italic type, and symbols of matrices are printed in bold type.</footnote>
where $b(\boldsymbol{x})$ is a Lipschitz function from $\mathbb{R}^n$ to $\mathbb{R}^n$, $\boldsymbol{\sigma}(\boldsymbol{x})$ is a function of class $C^2$ with bounded derivatives from $\mathbb{R}^n$ to $M_{n,d}(\mathbb{R})$, and $\boldsymbol{\eta}^\varepsilon$ is a $\mathbb{R}^d$-valued rapidly varying mean-reverting process, with a mean equal to zero, a unique invariant distribution, and a mean reversion time of the order of $\varepsilon^2$. More exactly, in this paper we address the case when $\boldsymbol{\eta}^\varepsilon$ is a multivariate $d$-dimensional Ornstein-Uhlenbeck process

$$d\boldsymbol{\eta}^{\varepsilon} = \frac{\mathbf{K}}{\varepsilon}d\mathbf{W}_t - \frac{\mathbf{A}}{\varepsilon^2}\boldsymbol{\eta}^{\varepsilon}dt, \qquad (2)$$

where $\mathbf{A}$ is a $d \times d$ matrix, whose eigenvalues have positive real parts, $\mathbf{K}$ is a $d \times d'$ matrix, and $\mathbf{W}$ is a $d'$-dimensional Brownian motion. This model is classical. It can be encountered in earthquake engineering [24] and also in finance [35]. It can model stationary Gaussian processes with very general spectra (see Section 3).

Our main motivation is to estimate a quantity of the form

$$I^{\varepsilon} \triangleq \mathbb{E}[F(\mathbf{X}^{\varepsilon})] \qquad (3)$$

for a fixed, small or moderate, parameter $\varepsilon$, for a smooth real-valued function $F$ defined on the space of continuous functions over $[0, T]$. We may think at $F(\mathbf{X}) = f(\mathbf{X}_T)$ where $f$ is smooth with polynomial growth, or $F(\mathbf{X}) = \int_0^T h(\mathbf{X}_s)ds + f(\mathbf{X}_T)$. By the Feynman-Kac formula it is possible to get the value of $I^\varepsilon$ for the model (1-2) by solving a parabolic equation, but this equation is formulated in a $d+n$-dimensional space and it possesses large terms (of order $\varepsilon^{-2}$) that give rapid fluctuations. These rapid fluctuations need to be resolved by the numerical scheme, which imposes to take a time step smaller than $\varepsilon^2$. The numerical resolution (with a finite difference method) is, therefore, challenging, if not impossible, and we look for other resolution methods. It is also possible to estimate $I^\varepsilon$ by a brute force Monte Carlo method. The Monte Carlo method, however, requires many simulations to get an accurate estimation, and each simulation requires to resolve the rapid fluctuations at the scale $\varepsilon^2$, so we would like to propose an efficient variance reduction method. The main idea is to find a limiting process $\mathbf{X}^0$ that approximates $\mathbf{X}^\varepsilon$ in a strong sense when $\varepsilon \to 0$ and for which the value

$$I^0 = \mathbb{E}[F(\mathbf{X}^0)] \qquad (4)$$

is known or can be estimated efficiently. It is then possible to propose a control variate method to estimate $I^\varepsilon$ for a fixed $\varepsilon$.

We consider the limiting $\mathbb{R}^n$-valued process $\mathbf{X}^0$ solution of the stochastic differential equation (SDE)

$$d\mathbf{X}^0 = \tilde{\mathbf{b}}(\mathbf{X}^0)dt + \Gamma(\mathbf{X}^0)d\mathbf{W}_t, \qquad (5)$$

where $\mathbf{X}^0$ share the same driving Brownian motion as $\boldsymbol{\eta}$, with the functions $\tilde{\mathbf{b}}(\boldsymbol{x})$ from $\mathbb{R}^n$ to $\mathbb{R}^n$ and $\Gamma(\boldsymbol{x})$ from $\mathbb{R}^n$ to $M_{n,d'}(\mathbb{R})$ given by

$$\tilde{b}_j(\boldsymbol{x}) \triangleq b_j(\boldsymbol{x}) + \sum_{i=1}^{n} ((\partial_{x_i} \boldsymbol{\sigma}(\boldsymbol{x}))\mathbf{A}^{-1}\mathbf{C}\boldsymbol{\sigma}(\boldsymbol{x})^T)_{ji}, \qquad (6)$$

$$\Gamma(\boldsymbol{x}) \triangleq \boldsymbol{\sigma}(\boldsymbol{x})\mathbf{A}^{-1}\mathbf{K}, \qquad (7)$$

and $\mathbf{C}$ is the $d \times d$ matrix defined by

$$\mathbf{C} \triangleq \int_{0}^{\infty} e^{-\mathbf{A}s} \mathbf{K} \mathbf{K}^{T} e^{-\mathbf{A}^{T}s} ds. \qquad (8)$$
The matrix **C** is the covariance matrix of the stationary distribution of the process $\eta^\epsilon$. We show in Proposition 3.5 that the continuous process ($X^\epsilon - X^0$) converges in probability to zero as $\epsilon \to 0$. The fact that the continuous process $X^\epsilon$ converges in distribution to $X^0$ is well-known [8, Chapter 6], but here we get a stronger result with a particular coupling between the two processes $X^\epsilon$ and $X^0$, that is needed to implement the control variate method that we have in mind.

The form of the limiting equation (5) is not surprising. Indeed, by (2), we can anticipate that $\frac{1}{\epsilon}\eta^\epsilon dt \simeq \mathbf{A}^{-1}\mathbf{K}d\mathbf{W}_t+\text{corrections}$, which explains the form (7) of the diffusion $\mathbf{\Gamma}$. The form (6) of the drift $\tilde{\mathbf{b}}$ is a manifestation of the Itô-versus-Stratonovich problem [32]. This problem is whether one should interpret the stochastic integral in the limiting equation in Itô sense, Stratonovich sense, or another sense. The Wong-Zakai theory [36] claims that the limiting diffusion should be a Stratonovich equation when $d=1$. Indeed, Eqs. (6-8) then reduce to $\mathbf{\Gamma}(\mathbf{x}) = \frac{1}{A}\sigma(\mathbf{x})\mathbf{K}$, $\mathbf{C} = \frac{1}{2A}\mathbf{K}\mathbf{K}^T$,

$$ \tilde{b}_j(\mathbf{x}) - b_j(\mathbf{x}) = \frac{1}{2A^2} \sum_{i=1}^{n} (\partial_{x_i}(\sigma(\mathbf{x})\mathbf{K})(\sigma(\mathbf{x})\mathbf{K})^T)_{ji} = \frac{1}{2} \sum_{i=1}^{n} (\partial_{x_i}\mathbf{\Gamma}(\mathbf{x})\mathbf{\Gamma}(\mathbf{x})^T)_{ji}, $$

so that (5) can be written as

$$ d\mathbf{X}^0 = \mathbf{b}(\mathbf{X}^0)dt + \mathbf{\Gamma}(\mathbf{X}^0) \circ d\mathbf{W}_t, \qquad (9) $$

where $\circ$ stands for the Stratonovich integral, because

$$ \begin{aligned} (\mathbf{\Gamma}(\mathbf{X}^0) \circ d\mathbf{W}_t)_j &= (\mathbf{\Gamma}(\mathbf{X}^0)d\mathbf{W}_t)_j + \frac{1}{2} \sum_{i=1}^{n} \sum_{j'=1}^{d'} \partial_{x_i} \Gamma_{jj'}(\mathbf{X}^0) d\langle X_i^0, W_{j'} \rangle_t \\ &= (\mathbf{\Gamma}(\mathbf{X}^0)d\mathbf{W}_t)_j + \frac{1}{2} \sum_{i=1}^{n} (\partial_{x_i} \mathbf{\Gamma}(\mathbf{X}^0) \mathbf{\Gamma}(\mathbf{X}^0)^T)_{ji} dt. \end{aligned} $$

The form (9) is valid when $d=1$ and it looks simpler than (5), but we have chosen to write the stochastic integral in (5) in Itô's sense and to add the appropriate Itô-Stratonovich drift correction $\tilde{\mathbf{b}}-\mathbf{b}$, because it is a natural starting point for numerical schemes [20] and it is the appropriate form to express the martingale problems used in the proofs (see Appendix). When $d > 1$ the difference between $\tilde{\mathbf{b}}$ and $\mathbf{b}$ is an Itô-Stratonovich correction that is more complex and the limiting equation (5) cannot be reduced to (9).

We can now introduce the Monte Carlo method for the estimation of $I^\epsilon$. Let $\mathbf{W}^k$, $k=1, \dots, N$, be $N$ independent and identically distributed $d'$-dimensional Brownian motions. We consider three Monte Carlo-type estimators of $I^\epsilon$:

1) The brute force Monte Carlo estimator is

$$ \hat{I}_N^\epsilon \triangleq \frac{1}{N} \sum_{k=1}^{N} F(\mathbf{X}^\epsilon(\mathbf{W}^k)), \qquad (10) $$

where $\mathbf{X}^\epsilon(\mathbf{W}^k)$ is the solution of (1-2) with $\mathbf{W}^k$. The estimator $\hat{I}_N^\epsilon$ is unbiased and its variance is

$$ \mathrm{Var}(\hat{I}_N^\epsilon) = \frac{1}{N}\mathrm{Var}(F(\mathbf{X}^\epsilon)). \qquad (11) $$
It is asymptotically normal as $N \to +\infty$:

$$ \sqrt{N}(\hat{I}_{N}^{\varepsilon} - I^{\varepsilon}) \xrightarrow{dist.} \mathcal{N}(0, \sigma_{I^{\varepsilon}}^{2}), \qquad (12) $$

with the asymptotic variance

$$ \sigma_{I^{\varepsilon}}^{2} = \mathrm{Var}(F(\mathbf{X}^{\varepsilon})), \qquad (13) $$

which has the following behavior as $\varepsilon \to 0$ when $F$ is continuous and bounded (because $\mathbf{X}^{\varepsilon}$ weakly converges to $\mathbf{X}^{0}$):

$$ \sigma_{I^{\varepsilon}}^{2} = \mathrm{Var}(F(\mathbf{X}^{0})) + o(1). \qquad (14) $$

2) The control variate estimator [12] is

$$ \hat{J}_{N}^{\varepsilon} \triangleq I^{0} + \frac{1}{N} \sum_{k=1}^{N} F(\mathbf{X}^{\varepsilon}(\mathbf{W}^{k})) - F(\mathbf{X}^{0}(\mathbf{W}^{k})), \qquad (15) $$

where $I^0 = \mathbb{E}[F(\mathbf{X}^0)]$ is supposed to be known exactly (or with high accuracy). The value $I^0$ can be obtained by solving a Kolmogorov equation in a $n$-dimensional framework and without large term; if this is not possible (because $n$ is too large for instance), then the value $I^0$ can be obtained by a brute force Monte Carlo method which is easier than for $I^\varepsilon$ because there is no large term of order $\varepsilon^{-2}$, so that a standard Euler scheme for stochastic differential equations can be used [20]. The control variate estimator $\hat{J}_N^\varepsilon$ is unbiased and its variance is

$$ \mathrm{Var}(\hat{J}_{N}^{\varepsilon}) = \frac{1}{N}\mathrm{Var}(F(\mathbf{X}^{\varepsilon}) - F(\mathbf{X}^{0})). \qquad (16) $$

It is asymptotically normal as $N \to +\infty$:

$$ \sqrt{N}(\hat{J}_{N}^{\varepsilon} - I^{\varepsilon}) \xrightarrow{dist.} \mathcal{N}(0, \sigma_{J^{\varepsilon}}^{2}), \qquad (17) $$

with the asymptotic variance

$$ \sigma_{J^{\varepsilon}}^{2} = \mathrm{Var}(F(\mathbf{X}^{\varepsilon}) - F(\mathbf{X}^{0})). \qquad (18) $$

When $F$ is continuous and bounded, we have by Proposition 3.5 that $\sigma_{J^\varepsilon}^2$ goes to zero as $\varepsilon \to 0$. More quantitatively, if $F(\mathbf{X}) = f(\mathbf{X}_T)$ for a smooth $f$ with bounded derivatives, then the asymptotic variance has the following behavior as $\varepsilon \to 0$ (by Lemma 3.8):

$$ \sigma_{J^{\varepsilon}}^{2} \leq C\varepsilon^{2}. \qquad (19) $$

The order of magnitude $\varepsilon^2$ of the asymptotic variance of $\hat{J}_N^\varepsilon$ is confirmed by the numerical simulations that we report in Section 4.

3) The theoretical optimal control variate estimator is

$$ \hat{O}_{N}^{\varepsilon} \triangleq \rho^{\varepsilon} I^{0} + \frac{1}{N} \sum_{k=1}^{N} F(\mathbf{X}^{\varepsilon}(\mathbf{W}^{k})) - \rho^{\varepsilon} F(\mathbf{X}^{0}(\mathbf{W}^{k})), \qquad (20) $$

with

$$ \rho^{\varepsilon} = \mathrm{Cov}(F(\mathbf{X}^{\varepsilon}), F(\mathbf{X}^{0})) / \mathrm{Var}(F(\mathbf{X}^{0})). \qquad (21) $$
This estimator is unbiased and has the minimal variance

$$ \mathrm{Var}(\hat{O}_N^\epsilon) = \frac{1}{N} \mathrm{Var}(F(\mathbf{X}^\epsilon) - \rho^\epsilon F(\mathbf{X}^0)), \qquad (22) $$

amongst all control variate estimators of the form

$$ \rho I^0 + \frac{1}{N} \sum_{k=1}^{N} F(\mathbf{X}^{\epsilon}(\mathbf{W}^k)) - \rho F(\mathbf{X}^0(\mathbf{W}^k)). $$

Note that $\rho = 0$ corresponds to the brute force Monte Carlo estimator $\hat{I}_N^\epsilon$, $\rho = 1$ corresponds to the control variate estimator $\hat{J}_N^\epsilon$, and $\rho = \rho^\epsilon$ corresponds to the optimal control variate estimator $\hat{O}_N^\epsilon$. The estimator $\hat{O}_N^\epsilon$ is asymptotically normal as $N \to +\infty$:

$$ \sqrt{N}(\hat{O}_N^\epsilon - I^\epsilon) \xrightarrow{\mathrm{dist.}} \mathcal{N}(0, \sigma_{O^\epsilon}^2), \qquad (23) $$

with the asymptotic variance

$$ \sigma_{O^\epsilon}^2 = \mathrm{Var}(F(\mathbf{X}^\epsilon) - \rho^\epsilon F(\mathbf{X}^0)). \qquad (24) $$

The estimator $\hat{O}_N^\epsilon$ is, however, not practical as it depends on $\rho^\epsilon$ which is unknown. The practical optimal control variate estimator [12] is

$$ \hat{K}_N^\epsilon \triangleq \hat{\rho}_N^\epsilon I^0 + \frac{1}{N} \sum_{k=1}^{N} F(\mathbf{X}^\epsilon(\mathbf{W}^k)) - \hat{\rho}_N^\epsilon F(\mathbf{X}^0(\mathbf{W}^k)), \qquad (25) $$

where $\hat{\rho}_N^\epsilon$ is the empirical correlation

$$ \hat{\rho}_N^\epsilon = \frac{\sum_{k=1}^{N} (F(\mathbf{X}^\epsilon(\mathbf{W}^k)) - \hat{I}_N^\epsilon)(F(\mathbf{X}^0(\mathbf{W}^k)) - \hat{I}_N^0)}{\sum_{k=1}^{N} (F(\mathbf{X}^0(\mathbf{W}^k)) - \hat{I}_N^0)^2}, \qquad (26) $$

with $\hat{I}_N^\epsilon = \frac{1}{N} \sum_{k=1}^N F(\mathbf{X}^\epsilon(\mathbf{W}^k))$ as in (10) and $\hat{I}_N^0 = \frac{1}{N} \sum_{k=1}^N F(\mathbf{X}^0(\mathbf{W}^k))$. This estimator is a practical and approximate version of the theoretical optimal control variate estimator $\hat{O}_N^\epsilon$ in which the unknown correlation coefficient $\rho^\epsilon$ has been replaced by its empirical estimator $\hat{\rho}_N^\epsilon$. The estimator $\hat{K}_N^\epsilon$ may be slightly biased and may have a variance slightly larger than (22) because of the empirical estimation of $\rho^\epsilon$. $\hat{K}_N^\epsilon$ is, however, asymptotically normal with an asymptotic variance that is the same one as that of the optimal estimator $\hat{O}_N^\epsilon$, as shown by the following proposition.

**Proposition 2.1.** As $N \to +\infty$,

$$ \sqrt{N}(\hat{K}_N^\epsilon - I^\epsilon) \xrightarrow{\mathrm{dist.}} \mathcal{N}(0, \sigma_{K^\epsilon}^2), \qquad (27) $$

with

$$ \sigma_{K^\epsilon}^2 = \sigma_{O^\epsilon}^2 = \mathrm{Var}(F(\mathbf{X}^\epsilon) - \rho^\epsilon F(\mathbf{X}^0)). \qquad (28) $$

Furthermore, if $F$ is continuous and bounded, then $\sigma_{K^\epsilon}^2$ goes to zero as $\epsilon \to 0$. If $F(\mathbf{X}) = f(\mathbf{X}_T)$, with $f$ with bounded derivatives, then there exists $C > 0$ such that

$$ \sigma_{K^\epsilon}^2 \le C\epsilon^2, \quad 0 \le \sigma_{J^\epsilon}^2 - \sigma_{K^\epsilon}^2 \le C\epsilon^4. \qquad (29) $$
Proof. By the law of large numbers, $\hat{\rho}_N^\epsilon$ converges to $\rho^\epsilon$ as $N \to +\infty$. The convergence holds almost surely, hence in probability. We have

$$ \hat{K}_N^\epsilon - I^\epsilon = (\hat{O}_N^\epsilon - I^\epsilon) - (\hat{\rho}_N^\epsilon - \rho^\epsilon)(\hat{I}_N^0 - I^0), $$

so we get (27-28) from Slutsky's theorem.
Furthermore, we have

$$ \sigma_{J^\epsilon}^2 - \sigma_{K^\epsilon}^2 = \mathrm{Var}(F(\mathbf{X}^0))(1 - \rho^\epsilon)^2. $$

If $F$ is continuous and bounded, then $\rho^\epsilon$ goes to one and $\sigma_{J^\epsilon}^2$ goes to zero as $\epsilon \to 0$ by Proposition 3.5. If $F(\mathbf{X}) = f(\mathbf{X}_T)$, then, by Lemma 3.8, $1 - \rho^\epsilon$ and $\sigma_{J^\epsilon}^2$ are of order $O(\epsilon^2)$ for small $\epsilon$. This shows the desired result (29). $\square$

Proposition 2.1 shows that the asymptotic variances of the estimators $\hat{K}_N^\epsilon$ and $\hat{J}_N^\epsilon$ are equivalent for vanishingly small $\epsilon$ and of the order of $O(\epsilon^2)$, and that the asymptotic variance of the estimator $\hat{K}_N^\epsilon$ is slightly smaller than that of $\hat{J}_N^\epsilon$ for moderately small $\epsilon$. These statements are confirmed by the numerical simulations that we report in Section 4.

In addition, motivated by the examples that we address in Section 7, we consider the case where the $\mathbb{R}^n$-valued process $\mathbf{X}^\epsilon$ satisfies a multivalued ODE of the form

$$ \frac{d\mathbf{X}^\epsilon}{dt} + \partial\varphi(\mathbf{X}^\epsilon) \ni \mathbf{b}(\mathbf{X}^\epsilon) + \frac{1}{\epsilon}\sigma\eta^\epsilon, \quad \mathbf{X}_0^\epsilon = \mathbf{x}_0, \qquad (30) $$

and the case where $\mathbf{X}^\epsilon$ together with a $\mathbb{R}^m$-valued process $\mathbf{Z}^\epsilon$ satisfy the multivalued ODE

$$ \left\{ \begin{array}{l} \frac{d\mathbf{X}^\epsilon}{dt} + \partial\varphi(\mathbf{X}^\epsilon) \ni \mathbf{b}^X(\mathbf{X}^\epsilon, \mathbf{Z}^\epsilon) + \frac{1}{\epsilon}\sigma\eta^\epsilon, \quad \mathbf{X}_0^\epsilon = \mathbf{x}_0, \\ \\ \frac{d\mathbf{Z}^\epsilon}{dt} + \partial\psi(\mathbf{Z}^\epsilon) \ni \mathbf{b}^Z(\mathbf{X}^\epsilon, \mathbf{Z}^\epsilon), \quad \mathbf{Z}_0^\epsilon = \mathbf{z}_0. \end{array} \right. \qquad (31) $$

Here $\sigma \in \mathcal{M}_{n,d}(\mathbb{R})$ is constant, $\mathbf{b}(\mathbf{x})$ from $\mathbb{R}^n$ to $\mathbb{R}^n$, $\mathbf{b}^Z(\mathbf{x}, \mathbf{z})$ from $\mathbb{R}^{n+m}$ to $\mathbb{R}^m$ and $\mathbf{b}^X(\mathbf{x}, \mathbf{z})$ from $\mathbb{R}^{n+m}$ to $\mathbb{R}^n$ are Lipschitz functions. The operators $\partial\varphi$ and $\partial\psi$ are the subdifferentials of some lower semi continuous (l.s.c.) convex functions $\varphi$ from $\mathbb{R}^n$ to $[0, +\infty]$ and $\psi$ from $\mathbb{R}^m$ to $[0, +\infty]$. Stronger hypotheses will be assumed on $\varphi$ compared to $\psi$ as explained in Section 5 and important examples motivate the two situations as shown in Section 7. It is important to observe that the multivalued operators that appear in the differential inclusions above are subdifferential of convex functions, therefore existence and uniqueness are guaranteed [4, page 72]. For the reader's convenience, proofs of existence and uniqueness are given in Appendix C. It is worth mentioning that there is an alternative formulation using the language of variational inequalities, that is equivalent to differential inclusions. Eq. (30) is equivalent to

$$ \forall \xi \in \mathbb{R}^n, \forall t > 0, \quad \left( \mathbf{b}(\mathbf{X}^\epsilon) + \frac{1}{\epsilon}\sigma\eta^\epsilon - \frac{d\mathbf{X}^\epsilon}{dt} \right) \cdot (\xi - \mathbf{X}^\epsilon) + \varphi(\mathbf{X}^\epsilon) \le \varphi(\xi), \quad \mathbf{X}_0^\epsilon = \mathbf{x}_0, $$

and Eq. (31) is equivalent to

$$ \begin{array}{l} \forall \xi \in \mathbb{R}^n, \forall \zeta \in \mathbb{R}^m, \forall t > 0, \\ \left( \mathbf{b}^X(\mathbf{X}^\epsilon, \mathbf{Z}^\epsilon) + \frac{1}{\epsilon}\sigma\eta^\epsilon - \frac{d\mathbf{X}^\epsilon}{dt} \right) \cdot (\xi - \mathbf{X}^\epsilon) + \varphi(\mathbf{X}^\epsilon) \le \varphi(\xi), \quad \mathbf{X}_0^\epsilon = \mathbf{x}_0, \\ \left( \mathbf{b}^Z(\mathbf{X}^\epsilon, \mathbf{Z}^\epsilon) - \frac{d\mathbf{Z}^\epsilon}{dt} \right) \cdot (\zeta - \mathbf{Z}^\epsilon) + \psi(\mathbf{Z}^\epsilon) \le \psi(\zeta), \quad \mathbf{Z}_0^\epsilon = \mathbf{z}_0. \end{array} $$
Propositions 5.1 and 5.2 show that the multi-valued process $X^\epsilon$ strongly converges to a limiting process solution of a multivalued SDE. Eqs. (64) and (65) show that the control variate estimators have asymptotic variances of order $\epsilon^2$ for (30) and $\epsilon$ for (31).

To demonstrate the efficiency of our method on a practical problem, we consider a two-degree of freedom (TDOF) system as shown in Figure 1. It can describe a broad class of TDOF structures, including a two-storey building as presented in [34, Figure 4.6(a)].

<figure><img src="image_2.png" /><figcaption>Figure 1: A rheological model of a two-degree of freedom system. Two masses $m_1$ and $m_3$ are associated in series with elements which are themselves an association of dampers and springs. Elements ① and ③ represent a spring and a damper respectively, both possibly nonlinear or hysteretic. Here $c_1$ is a damping coefficient associated to the linear damper connecting the mass $m_1$ to the foundation and $k_3$ is a stiffness coefficient of the linear spring linking the masses $m_1$ and $m_3$. A random forcing $\epsilon^{-1}\eta^\epsilon$ is applied to the mass $m_1$ (e.g. wind forces on a two-storey building).</figcaption></figure>

When the external force $\eta^\epsilon$ is a colored noise such as an Ornstein-Uhlenbeck process, the equation of motion can be written in the form of Equation (1) with $n = 4$, where $(X_1^\epsilon, X_2^\epsilon)$, resp. $(X_3^\epsilon, X_4^\epsilon)$, represents the position and the velocity of the mass $m_1$, resp. $m_3$, shown in Figure 1. Many nonlinear behaviours enter into this framework, we have in mind a nonlinear spring of the linear-plus-quadratic cubic type and a nonlinear damper of the linear-plus-quadratic type (see Example 4.3 and Figure 2). Similarly Equations (30) and (31) arise in the description of nonlinear behaviours with hysteresis such as elasto-plasticity and friction, see [34, Chapter 8] and Section 7 (see Example 7.3 and Figure 3). In Figures 2-3 we compare the behaviors of the brute force Monte Carlo estimator $\hat{I}_N^\epsilon$ with the ones of the control variate estimators $\hat{J}_N^\epsilon$ and $\hat{K}_N^\epsilon$. We also plot the empirical estimators of the asymptotic variances of the estimators $\hat{I}_N^\epsilon, \hat{J}_N^\epsilon$, and $\hat{K}_N^\epsilon$ as described in Subsection 7.3. In addition, for each value of $\epsilon \in \{0.1, 0.5, 0.9\}$, error bars (95% confidence interval) are shown for each of the estimators (in $\hat{I}_N^\epsilon, \hat{J}_N^\epsilon, \hat{K}_N^\epsilon$ order from the left to the right). Here $I^0$ is obtained by a massive Monte Carlo estimation of the limit process, which is possible with a coarse grid step as there is no large term involved. We can observe that the control variate estimator $\hat{K}_N^\epsilon$ has always the minimal variance. When $\epsilon$ is small and the original system and the limit system are poorly correlated $\rho^\epsilon \approx 0$, it behaves as the standard Monte Carlo estimator $\hat{I}_N^\epsilon$. When $\epsilon$ is small and the original system and the limit system are strongly correlated $\rho^\epsilon \approx 1$, it behaves as the control variate estimator $\hat{J}_N^\epsilon$. We can also observe that the variance reduction is by a factor of order $\epsilon^2$ when the quantity to be estimated is the expectation of a smooth function, while it is of order $\epsilon$ when the quantity to be estimated is the expectation of an indicator function.
<figure><img src="image_3.png" /><figcaption>Figure 2: Example 4.3 of a TDOF modeling a two-storey building with nonlinear spring of the linear-plus-quadratic cubic type and nonlinear damper of the linear-plus-quadratic type, driven by an Ornstein-Uhlenbeck noise. In the top row the target is to estimate $I^\epsilon = \mathbb{E}[(X_{1,T}^\epsilon)^2 + (X_{3,T}^\epsilon)^2]$ for $T = 1$ and the expectation of the control variate $I^0 = \mathbb{E}[(X_{1,T}^0)^2 + (X_{3,T}^0)^2]$ is obtained by an intensive Monte Carlo computation with coarse time step. In the bottom row the target is to estimate $I^\epsilon = \mathbb{P}(|X_{1,T}^\epsilon| \le a, |X_{3,T}^\epsilon| \le b)$ for $T = 1, a = 0.1, b = 0.1$ and the expectation of the control variate $I^0 = \mathbb{P}(|X_{1,T}^0| \le a, |X_{3,T}^0| \le b)$ is also obtained by an intensive Monte Carlo computation with coarse time step. The numerical procedure is standard and is described in Section 4 (Euler-Maruyama time discretization with time step $\delta t = 10^{-5}$). Here the number of Monte Carlo samples is $N = 10^4$ and $m_1 = m_3 = c_1 = k_3 = 1$. The details of the elements ① and ③ with nonlinear behaviours can be found in Example 4.3.</figcaption></figure>

### 3 Diffusion approximation for a driving multivariate Ornstein-Uhlenbeck process

We consider the $\mathbb{R}^n$-valued process $X^\epsilon$ solution of the ODE (1) when $\eta^\epsilon$ is the multivariate $d$-dimensional Ornstein-Uhlenbeck process (2). We give several explicit examples.

**Example 3.1.** $\eta^\epsilon$ is a one-dimensional Ornstein-Uhlenbeck process, $d = d' = 1$, $A, K > 0$,

$$d\eta^\epsilon = -\frac{A}{\epsilon^2}\eta^\epsilon dt + \frac{K}{\epsilon}dW_t. \qquad (32)$$

**Example 3.2.** $\eta^\epsilon$ is a Langevin process

$$d\eta_1^\epsilon = \frac{1}{\epsilon^2}\eta_2^\epsilon dt, \qquad (33)$$

$$d\eta_2^\epsilon = -\frac{1}{\epsilon^2}[\mu\eta_1^\epsilon + \gamma\eta_2^\epsilon]dt + \frac{K}{\epsilon}dW_t, \qquad (34)$$
<figure><img src="image_4.png" /><figcaption>Figure 3: Example 7.3 of a TDOF modeling a two-storey building with an hysteretic spring of elasto-plastic type, driven by an Ornstein-Uhlenbeck noise. The quantities presented here are similar to those presented in Figure 2. The details of the elements ① and ③ with nonlinear and hysteretic behaviours can be found in Example 7.3.</figcaption></figure>

which corresponds to $d = 2, d' = 1$, $A = \begin{pmatrix} 0 & -1 \\ \mu & \gamma \end{pmatrix}$, and $K = \begin{pmatrix} 0 \\ K \end{pmatrix}$. The process $\eta_1^\epsilon$ is a white-noise driven linear oscillator with stiffness $\mu > 0$ and damping $\gamma > 0$. It can be encountered in earthquake engineering because it is considered to be a realistic type of random forcing to represent seismic excitation (it is the so-called Kanai-Tajimi model [24]).

**Example 3.3.** If $\eta^\epsilon$ is a real-valued zero-mean stationary Gaussian process with power spectral density $\text{PSD}^\epsilon(\omega) = \epsilon^2\text{PSD}(\epsilon^2\omega)$, $\text{PSD}(\omega) = \sum_{k=1}^q \frac{\sigma_k^2}{1+\omega^2/\Delta\Omega_k^2}$, then it has the same distribution as the process $\sum_{k=1}^q \sigma_k \eta_k^\epsilon$ where $\eta^\epsilon$ is solution of (2) with $d = d' = q$ and

$$A = K = \text{diag}(\Delta\Omega_k, k = 1, \dots, q).$$

This shows that any zero-mean stationary Gaussian process with power spectral density that can be decomposed as a sum of centered Lorentzians belongs to the model (2).

**Example 3.4.** If $\eta^\epsilon$ is a real-valued zero-mean stationary Gaussian process with power spectral density $\text{PSD}^\epsilon(\omega) = \epsilon^2\text{PSD}(\epsilon^2\omega)$, $\text{PSD}(\omega) = \frac{1}{2} \sum_{k=1}^q \frac{\sigma_k^2}{1+(\omega-\omega_k)^2/\Delta\Omega_k^2} + \frac{\sigma_k^2}{1+(\omega+\omega_k)^2/\Delta\Omega_k^2}$, then it has the same distribution as the process $\sum_{k=1}^q \sigma_k \eta_{2k-1}^\epsilon$ where $\eta^\epsilon$ is solution of (2) with $d = d' = 2q$ and

$$A = \bigoplus_{k=1}^{q} \begin{pmatrix} \Delta\Omega_k & -\omega_k \\ \omega_k & \Delta\Omega_k \end{pmatrix}, \quad K = \bigoplus_{k=1}^{q} \begin{pmatrix} \Delta\Omega_k & 0 \\ 0 & \Delta\Omega_k \end{pmatrix}.$$

This shows that any zero-mean stationary Gaussian process with power spectral density that can be decomposed as a sum of non-centered Lorentzian functions belongs to the model (2).
We also consider the limiting $\mathbb{R}^n$-valued process $\mathbf{X}^0$ solution of the SDE (5). The continuous process $(\mathbf{X}^\epsilon - \mathbf{X}^0)$ converges in probability to zero as $\epsilon \to 0$ as stated in the following proposition.

**Proposition 3.5.** If $\mathbf{X}_0^\epsilon = \mathbf{X}_0^0$, then the continuous process $(\mathbf{X}^\epsilon - \mathbf{X}^0)$ converges in probability to zero as $\epsilon \to 0$. The convergence holds in the space of continuous functions equipped with the topology associated to the uniform norm over compact intervals.

The proof of Proposition 3.5 is based on the perturbed test function method as described first in [22, Chapter 7] or in [8, Chapter 6]. It is given in Appendix A.

**Example 3.6.** We consider the process $\mathbf{X}^\epsilon$ solution of the ODE (1) where $\eta^\epsilon$ is the rapidly varying mean-reverting process (32). We also consider the limiting process

$$d\mathbf{X}^0 = \mathbf{b}(\mathbf{X}^0)dt + \frac{K}{A}\sigma(\mathbf{X}^0)dW_t + \frac{K^2}{2A^2}(\sigma(\mathbf{X}^0) \cdot \nabla_{\mathbf{x}^0})\sigma(\mathbf{X}^0)dt,$$

driven by the same Brownian motion. The continuous process $(\mathbf{X}^\epsilon - \mathbf{X}^0)$ converges in probability to zero as $\epsilon \to 0$.

**Example 3.7.** We consider the process $\mathbf{X}^\epsilon$ solution of the ODE (1) where $\eta^\epsilon$ is the rapidly varying mean-reverting process (33-34). We also consider the limiting process

$$d\mathbf{X}^0 = \mathbf{b}(\mathbf{X}^0)dt + \frac{K}{\mu}\sigma(\mathbf{X}^0)dW_t + \frac{K^2}{2\mu^2}(\sigma(\mathbf{X}^0) \cdot \nabla_{\mathbf{x}^0})\sigma(\mathbf{X}^0)dt,$$

driven by the same Brownian motion. The continuous process $(\mathbf{X}^\epsilon - \mathbf{X}^0)$ converges in probability to zero as $\epsilon \to 0$.

The proof that the optimal control variate estimator $\hat{K}_N^\epsilon$ and the control variate estimator $\hat{J}_N^\epsilon$ have asymptotic variances of the order of $\epsilon^2$ as stated in Proposition 2.1 follows from the following lemma.

**Lemma 3.8.** Let $f,g$ be smooth functions from $\mathbb{R}^n$ to $\mathbb{R}$ with bounded derivatives. Let $T > 0$. There exists $C > 0$ such that, for any $t \in [\epsilon, T]$,

$$|\mathbb{E}[g(\mathbf{X}_t^0)(f(\mathbf{X}_t^\epsilon) - f(\mathbf{X}_t^0))]| \le C\epsilon^2, \quad \mathbb{E}[(f(\mathbf{X}_t^\epsilon) - f(\mathbf{X}_t^0))^2] \le C\epsilon^2. \quad (35)$$

The important hypothesis is that $f$ should be smooth. We could certainly relax the hypothesis on the bounded derivatives by using uniform estimates of high-order moments of the process $\mathbf{X}^\epsilon$. Lemma 3.8 is proved in Appendix B.

# 4 Numerical simulations

In this section, we illustrate our control variate method and report the numerical results on different types of dynamical systems driven by colored noises. The two examples are smooth oscillators that can be described by Equation (1) (one being linear with time-dependent coefficients and the other being of Van der Pol type). Other examples with non-smooth dynamical systems will be addressed in Section 7.

We use the Euler-Maruyama approximation method to compute the approximate numerical solution of a SDE [20]. In Subsection 4.1, we recall the two types of colored noise that we consider and provide their time discretization. Then, in Subsection 4.2, some details and discretization of the dynamical systems under consideration are given. Finally, in Subsection 4.3, numerical experiments on the control variate estimators are provided and discussed in each case.
## 4.1 Colored noise models and their discretization

The two models of noise are shown in Eq. (32) (OU) and in the system of equations (33-34) (Langevin). The OU noise has two parameters $A, K_{ou} > 0$ whereas the Langevin has three parameters $\mu, \gamma, K_{lan} > 0$. Their discretization works as follows. Let $T > 0$ and $N_T \in \mathbb{N}$ be the number of time steps such that $T = N_T \delta t$. Let $N$ be the number of Monte Carlo samples. Consider a sequence of independent and identically distributed standard Gaussian variables

$$ \{\Delta W_n^k \sim \mathcal{N}(0, 1), 0 \le n \le N_T - 1, 1 \le k \le N\}. $$

Let $\varepsilon > 0$. For each $1 \le k \le N$, we overload the notation by denoting the discretized noise in both cases by $\{\hat{\eta}_n^{\varepsilon,k}, 0 \le n \le N_T\}$.

* **Ornstein-Uhlenbeck noise:** $\hat{\eta}_0^{\varepsilon,k} \sim \mathcal{N}\left(0, \frac{K_{ou}^2}{2A}\right)$ and for $0 \le n \le N_T - 1$,

$$ \hat{\eta}_{n+1}^{\varepsilon,k} = \hat{\eta}_{n}^{\varepsilon,k} \left(1 - \delta t \frac{A}{\varepsilon^2}\right) + \sqrt{\delta t} \frac{K_{ou}}{\varepsilon} \Delta W_n^k. $$

* **Langevin noise:** $\hat{\eta}_0^{\varepsilon,k}$ and $\hat{\eta}_{2,0}^{\varepsilon,k}$ are independent variables with

$$ \hat{\eta}_0^{\varepsilon,k} \sim \mathcal{N}\left(0, \frac{K_{lan}^2}{2\gamma}\right), \quad \hat{\eta}_{2,0}^{\varepsilon,k} \sim \mathcal{N}\left(0, \frac{K_{lan}^2}{2\gamma\mu}\right) $$

and for $0 \le n \le N_T - 1$,

$$ \hat{\eta}_{n+1}^{\varepsilon,k} = \hat{\eta}_{n}^{\varepsilon,k} + \frac{\delta t}{\varepsilon^2} \hat{\eta}_{2,n}^{\varepsilon,k}, \quad \hat{\eta}_{2,n+1}^{\varepsilon,k} = \hat{\eta}_{2,n}^{\varepsilon,k} - \frac{\delta t}{\varepsilon^2} [\mu \hat{\eta}_{n}^{\varepsilon,k} + \gamma \hat{\eta}_{2,n}^{\varepsilon,k}] + \sqrt{\delta t} \frac{K_{lan}}{\varepsilon} \Delta W_n^k. $$

## 4.2 Details and discretization of the illustrative dynamical systems

We consider systems of the form of (1). We first consider the case of smooth systems that can have time-dependent coefficients,

$$ \frac{dX_1^\varepsilon}{dt} = X_2^\varepsilon, \quad \frac{dX_2^\varepsilon}{dt} = -h(X_1^\varepsilon, X_2^\varepsilon, t) + \frac{1}{\varepsilon}\eta^\varepsilon. \qquad (36) $$

Here we are interested in $\mathbb{E}[\|\mathbf{X}_T^\varepsilon\|^2]$ and in $\mathbb{P}(|X_{1,T}^\varepsilon| \le 1)$ for $T = 1$. Note that the second case corresponds to an expectation $\mathbb{P}(|X_{1,T}^\varepsilon| \le 1) = \mathbb{E}[f(\mathbf{X}_T^\varepsilon)]$ with a non-smooth function $f(\mathbf{x}) = \mathbf{1}_{|x_1| \le 1}$. As $\varepsilon \to 0$, $\mathbf{X}^\varepsilon = (X_1^\varepsilon, X_2^\varepsilon)$ converges to $\mathbf{X}^0 = (X_1^0, X_2^0)$ where

$$ dX_1^0 = X_2^0 dt, \quad dX_2^0 = -h(X_1^0, X_2^0, t)dt + CdW, \qquad (37) $$

$C = K_{ou}A^{-1}$ for an OU noise and $C = K_{lan}\mu^{-1}$ for a Langevin noise. For the stochastic simulation of (36) and (37), we proceed as follows:

* $\hat{X}_{1,0}^{\varepsilon,k} = x_{1,0}$, $\hat{X}_{2,0}^{\varepsilon,k} = x_{2,0}$ and for $0 \le n \le N_T - 1$,

$$ \left\{ \begin{array}{l} \hat{X}_{1,n+1}^{\varepsilon,k} = \hat{X}_{1,n}^{\varepsilon,k} + \delta t \hat{X}_{2,n}^{\varepsilon,k}, \\ \hat{X}_{2,n+1}^{\varepsilon,k} = \hat{X}_{2,n}^{\varepsilon,k} - \delta t h(\hat{X}_{1,n}^{\varepsilon,k}, \hat{X}_{2,n}^{\varepsilon,k}, n\delta t) + \frac{\delta t}{\varepsilon} \hat{\eta}_n^{\varepsilon,k}. \end{array} \right. $$
• $\hat{X}_{1,0}^{0,k} = x_{1,0}$, $\hat{X}_{2,0}^{0,k} = x_{2,0}$ and for $0 \le n \le N_T - 1$,

$$ \begin{cases} \hat{X}_{1,n+1}^{0,k} = \hat{X}_{1,n}^{0,k} + \delta t \hat{X}_{2,n}^{0,k}, \\ \hat{X}_{2,n+1}^{0,k} = \hat{X}_{2,n}^{0,k} - \delta t \hat{X}_{1,n}^{0,k}, \end{cases} $$

$\hat{X}_n^{\varepsilon,k}$ and $\hat{X}_n^{0,k}$ are independent (in $k$) copies that are meant to approximate $X_{n\delta t}^{\varepsilon}$ and $X_{n\delta t}^{0}$.

**Example 4.1 (linear oscillator with time-dependent coefficients).** We take $h(x_1, x_2, t) \triangleq p(t)x_1 + q(t)x_2$ where $p(t) \triangleq 1 + \cos(t)$ and $q(t) \triangleq 1 + \sin(t)$ (the choice is purely arbitrary). Here, in both OU and Langevin cases, the limiting process $\mathbf{X}^0 = (X_1^0, X_2^0)$ is a Gaussian process provided that the initial condition is deterministic or Gaussian. This is useful to derive the expectation of the control variate. The distribution of $\mathbf{X}_t^0 = (X_{1,t}^0, X_{2,t}^0)$ is characterized by its first-order moment $\mathbf{m}(t) \triangleq \mathbb{E}[\mathbf{X}_t^0] \in \mathbb{R}^2$ and second-order moment $\mathbf{M}(t) \triangleq (\mathbb{E}[X_{i,t}^0 X_{j,t}^0])^2_{i,j=1} \in \mathcal{M}_{2,2}(\mathbb{R})$ which satisfy the following systems of differential equations:

• first-order moment

$$ - (m_1(0), m_2(0)) = (x_0, \dot{x}_0), $$

$$ - \dot{m}_1(t) = m_2(t), $$

$$ - \dot{m}_2(t) = -p(t)m_1(t) - q(t)m_2(t). $$

• second-order moment

$$ \left\{ \begin{aligned} & (M_{11}(0), M_{22}(0), M_{12}(0)) = (x_0^2, \dot{x}_0^2, x_0 \dot{x}_0), \\ & \dot{M}_{11}(t) = 2M_{12}(t), \\ & \dot{M}_{22}(t) = -2p(t)M_{12}(t) - 2q(t)M_{22}(t) + C^2, \\ & \dot{M}_{12}(t) = M_{22}(t) - p(t)M_{11}(t) - q(t)M_{12}(t). \end{aligned} \right. \qquad (38) $$

The expectation of the control variate $\mathbb{E}[\|\mathbf{X}_T^0\|^2]$ with $T=1$ is estimated by solving numerically, with an Euler method, the differential equations for the first- and second-order moments.

**Example 4.2 (Van der Pol oscillator).** We take $h(x_1, x_2) = x_1 - \nu(1 - x_1^2)x_2$ where $\nu > 0$. The expectation of the control variate can be represented by $\mathbb{E}[\|\mathbf{X}_T^0\|^2] = \mathbf{c}(\mathbf{x}_0, 0)$ with $T=1$, where $\mathbf{c}$ satisfies the following backward in time PDE

$$ \left\{ \begin{aligned} & \partial_t \mathbf{c} + \frac{C^2}{2} \partial_{x_2}^2 \mathbf{c} - h(x_1, x_2) \partial_{x_2} \mathbf{c} + x_2 \partial_{x_1} \mathbf{c} = 0, \quad \text{in } \mathbb{R}^2 \times [0, 1) \\ & \mathbf{c}(\mathbf{x}, 1) = \|\mathbf{x}\|^2 \quad \text{in } \mathbb{R}^2. \end{aligned} \right. \qquad (39) $$

The expectation $\mathbb{E}[\|\mathbf{X}_T^0\|^2]$ is estimated by solving the PDE (39) with a finite difference method.

**Example 4.3 (nonlinear TDOF).** A TDOF modeling a nonlinear spring of the linear-plus-quadratic cubic type and a nonlinear damper of the linear-plus-quadratic type, can be seen as a coupling between two systems of the form (36)

$$ \left\{ \begin{aligned} \frac{dX_1^\varepsilon}{dt} &= X_2^\varepsilon, & \frac{dX_2^\varepsilon}{dt} &= -g_2(X_1^\varepsilon, X_2^\varepsilon, X_3^\varepsilon, X_4^\varepsilon) + \frac{1}{\varepsilon}\eta^\varepsilon, \\ \frac{dX_3^\varepsilon}{dt} &= X_4^\varepsilon, & \frac{dX_4^\varepsilon}{dt} &= -g_4(X_1^\varepsilon, X_2^\varepsilon, X_3^\varepsilon, X_4^\varepsilon). \end{aligned} \right. \qquad (40) $$
In addition to $\mathbb{E}[(X_{1,T}^{\varepsilon})^2 + (X_{3,T}^{\varepsilon})^2]$ for $T = 1$, we are interested in $\mathbb{P}(|X_{1,T}^{\varepsilon}| \le a, |X_{3,T}^{\varepsilon}| \le b)$. As $\varepsilon \to 0$, $\mathbf{X}^{\varepsilon} \to \mathbf{X}^0$ where

$$ \left\{ \begin{array}{ll} \mathrm{d}X_1^0 = X_2^0 \mathrm{d}t, & \mathrm{d}X_2^0 = -g_2(X_1^0, X_2^0, X_3^0, X_4^0)\mathrm{d}t + C\mathrm{d}W, \\ \\ \mathrm{d}X_3^0 = X_4^0 \mathrm{d}t, & \mathrm{d}X_4^0 = -g_4(X_1^0, X_2^0, X_3^0, X_4^0)\mathrm{d}t. \end{array} \right. \qquad (41) $$

Here

$$ g_2(x_1, x_2, x_3, x_4) \triangleq k_1 x_1 (1 + \min(x_1^2, \tilde{L})) + c_1 x_2 - k_3 (x_3 - x_1) - c_3 (x_4 - x_2) (1 + \min(|x_4 - x_2|, \tilde{L})) $$

and

$$ g_4(x_1, x_2, x_3, x_4) \triangleq k_3(x_3 - x_1) + c_3(x_4 - x_2)(1 + \min(|x_4 - x_2|, \tilde{L})). $$

In the original model of Spanos $\tilde{L} = \infty$, see pages 189-190 in [34]. For any positive finite value of $\tilde{L}$, the system above enters into the scope of our results. The simulation of (40) and (41) is similar to what is done for (36) and (37). We take $\tilde{L} = 1000$, $c_1 = c_3 = k_1 = k_3 = 1$.

## 4.3 Numerical experiments

We report our numerical results for the two systems mentioned above. In each of the two figures below, there are four subfigures (a)-(b)-(c)-(d). For subfigures (a) and (b), the driving force is an Ornstein-Uhlenbeck noise (32) with $A = K = 1$. In subfigure (a), the dashed black, dotted blue, and solid red lines represent the standard MC estimator $\hat{I}_N^\varepsilon$ and the control variate estimators $\hat{K}_N^\varepsilon$, respectively. For each value of $\varepsilon \in \{0.1, 0.5, 0.9\}$, error bars (95% confidence interval) are shown for each of the estimators (in $\hat{I}_N^\varepsilon$, $\hat{J}_N^\varepsilon$, $\hat{K}_N^\varepsilon$ order from the left to the right). The black dotted line represents the expectation of the control variate. The objective of the subfigure (b) is to illustrate the bound (19) and to show that the $\varepsilon^2$-behavior is actually sharp. The same description applies to (c) and (d), except they correspond to the case of a Langevin noise (33-34) with $\mu = \gamma = K = 1$. In the figures the asymptotic variance of the standard MC estimator $\hat{I}_N^\varepsilon$ is estimated by

$$ \hat{\sigma}_{I^{\varepsilon}, N}^2 = \frac{1}{N} \sum_{k=1}^{N} \left( F(\mathbf{X}^{\varepsilon}(\mathbf{W}^k)) \right)^2 - (\hat{I}_N^{\varepsilon})^2, \qquad (42) $$

the asymptotic variance of the control variate estimator $\hat{J}_N^\varepsilon$ is estimated by

$$ \hat{\sigma}_{J^{\varepsilon}, N}^2 = \frac{1}{N} \sum_{k=1}^{N} \left( F(\mathbf{X}^{\varepsilon}(\mathbf{W}^k)) - F(\mathbf{X}^0(\mathbf{W}^k)) + I^0 \right)^2 - (\hat{J}_N^{\varepsilon})^2, \qquad (43) $$

and the asymptotic variance of the optimal control variate estimator $\hat{K}_N^\varepsilon$ is estimated by

$$ \hat{\sigma}_{K^{\varepsilon}, N}^2 = \frac{1}{N} \sum_{k=1}^{N} \left( F(\mathbf{X}^{\varepsilon}(\mathbf{W}^k)) - \hat{\rho}_N^{\varepsilon} F(\mathbf{X}^0(\mathbf{W}^k)) + \hat{\rho}_N^{\varepsilon} I^0 \right)^2 - (\hat{K}_N^{\varepsilon})^2, \qquad (44) $$

with $\hat{\rho}_N^\varepsilon$ defined by (26). $\hat{\sigma}_{I^{\varepsilon}, N}^2$, $\hat{\sigma}_{J^{\varepsilon}, N}^2$, and $\hat{\sigma}_{K^{\varepsilon}, N}^2$ are consistent estimators of $\sigma_{I^{\varepsilon}}^2$, $\sigma_{J^{\varepsilon}}^2$, and $\sigma_{K^{\varepsilon}}^2$, respectively.

We use $N = 10^4$ samples with a time step of $\delta t = 10^{-5}$ (note that $\delta t/\varepsilon^2 = 0.1$ for the smallest $\varepsilon = 10^{-2}$ used in the numerical results). We report the numerical results for the linear oscillator
with time-dependent coefficients in Figure 4 and for the Van der Pol oscillator in Figure 5. The numerical results concern the estimation of $I^\varepsilon = \mathbb{E}[\|\mathbf{X}_T^\varepsilon\|^2]$ or $I^\varepsilon = \mathbb{P}(|X_{1,T}^\varepsilon| \le 1)$ with $T=1$ where $\mathbf{X}^\varepsilon$ satisfies (36) and thus the expectation of the control variate is $I^0 = \mathbb{E}[\|\mathbf{X}_T^0\|^2]$ or $\mathbb{P}(|X_{1,T}^0| \le 1)$ where $\mathbf{X}^0$ satisfies (37).

<figure><img src="image_5.png" /><figcaption>Figure 4: Example 4.1 (linear oscillator with time-dependent coefficients) with $h(x_1, x_2, t) \triangleq p(t)x_1 + q(t)x_2$, $p(t) \triangleq 1 + \cos(t)$ and $q(t) \triangleq 1 + \sin(t)$. The target is to estimate $I^\varepsilon = \mathbb{E}[\|\mathbf{X}_T^\varepsilon\|^2]$ for $T=1$ and the expectation of the control variate $I^0 = \mathbb{E}[\|\mathbf{X}_T^0\|^2]$ is obtained by solving the set of differential equations (38).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5: Example 4.2 (Van der Pol oscillator) with $h(x_1, x_2) \triangleq x_1 - (1 - x_1^2)x_2$. In the top row the target is to estimate $I^\varepsilon = \mathbb{E}[\|\mathbf{X}_T^\varepsilon\|^2]$ for $T=1$ and the expectation of the control variate $I^0 = \mathbb{E}[\|\mathbf{X}_T^0\|^2]$ is obtained by solving the PDE (39). In the bottom row the target is to estimate $I^\varepsilon = \mathbb{P}(|X_{1,T}^\varepsilon| \le 1)$ for $T=1$ and the expectation of the control variate $I^0 = \mathbb{P}(|X_{1,T}^0| \le 1)$ is obtained by solving the PDE (39) with the suitable final condition.</figcaption></figure>

The theoretical predictions provided by Eq. (19) and Proposition 2.1 are based on the condition that $f$ has bounded derivatives. As we have discussed above, the assumption that $f$ is smooth is important but the hypothesis on the boundedness of the derivatives can certainly be relaxed. The numerical results shown in Figures 4-5 are actually in good agreement with the theoretical predictions: the asymptotic variances $\sigma_{J^\varepsilon}^2$ and $\sigma_{K^\varepsilon}^2$ behave as $O(\varepsilon^2)$. The only cases where the behavior is $O(\varepsilon)$, and not $O(\varepsilon^2)$, are when the quantity of interest is of the form $\mathbb{E}[f(\mathbf{X}_T^\varepsilon)]$ with a
function f that is not smooth, which is not surprising.

# 5 Diffusion approximation for the multivalued case

In this section we consider multivalued ODEs of the form (30) or (31).

## 5.1 Basic properties on differential inclusions

We recall that the subdifferential of a convex function $F: \mathbb{R}^r \to (-\infty, \infty]$ such that $\text{Dom}(F) = \{x \in \mathbb{R}^r, F(x) < \infty\}$ is not empty, is the map from $\mathbb{R}^r$ to $\mathcal{P}(\mathbb{R}^r)$ (the set of subsets of $\mathbb{R}^r$) defined by $\partial F(x) = \{\xi \in \mathbb{R}^r, \forall z \in \mathbb{R}^r, \xi \cdot (z - x) + F(x) \le F(z)\}$ for $x \in \text{Dom}(F)$ and $\partial F(x) = \emptyset$ for $x \notin \text{Dom}(F)$. To grasp quickly the idea when $r = 1$, $\partial F(x)$ can be seen as the set of sub-slopes of $F$ at the point $x$ and when $F$ is differentiable at the point $x$, $\partial F(x) = \{F'(x)\}$. See [4] for more details.

One way to construct a solution to a multivalued ODE of the form (30) or (31) is to proceed by penalization. The inclusion is replaced by an equality involving the Moreau-Yosida regularisation of $F: \mathbb{R}^r \to (-\infty, +\infty]$ (with $F = \varphi, r = n$ or $F = \psi, r = m$), that is

$$ \forall p \ge 1, \forall x \in \mathbb{R}^r, F_p(x) \triangleq \inf_{z \in \mathbb{R}^r} \left\{ F(z) + \frac{p}{2} \|x - z\|^2 \right\}. \quad (45) $$

We recall from Annex B in [30] some properties of $F_p$:

1. $F_p: \mathbb{R}^r \to \mathbb{R}$ is a convex differentiable function,

2. $\forall x \in \mathbb{R}^r$, $\partial F_p(x) = \{\nabla F_p(x)\}$ and $\nabla F_p(x) \in \partial F(J_p x)$ where $J_p x \triangleq x - \frac{1}{p} \nabla F_p(x)$,

3. $\exists C > 0, \forall x \in \mathbb{R}^r, \forall p, \|J_p x\| \le \|x\| + C$,

4. $\forall x, y \in \mathbb{R}^r, \|\nabla F_p(x) - \nabla F_p(y)\| \le p\|x - y\|$,

5. $\forall x, y \in \mathbb{R}^r, (\nabla F_p(x) - \nabla F_p(y)) \cdot (x - y) \ge 0$,

6. $\forall x \in \mathbb{R}^r$,

$$ x \cdot \nabla F_p(x) \ge 0, \quad (46) $$

7. $\forall x, y \in \mathbb{R}^r$,

$$ (\nabla F_p(x) - \nabla F_{p'}(y)) \cdot (x - y) \ge - \left( \frac{1}{p} + \frac{1}{p'} \right) \nabla F_p(x) \cdot \nabla F_{p'}(y), \quad (47) $$

8. as a consequence of properties 2 and 3 above, we also have

$$ \sup_{p \ge 1} \sup_{x \in \mathbb{R}^r} \frac{\|\nabla F_p(x)\|}{p(1 + \|x\|)} < \infty. \quad (48) $$

Thus, the penalized versions of (30) and (31) are

$$ \frac{d\mathbf{X}^{\varepsilon,p}}{dt} + \nabla\varphi_p(\mathbf{X}^{\varepsilon,p}) = \mathbf{b}(\mathbf{X}^{\varepsilon,p}) + \frac{1}{\varepsilon}\sigma\eta^{\varepsilon}, \quad \mathbf{X}_0^{\varepsilon,p} = \mathbf{x}_0, \quad (49) $$
and

$$ \left\{ \begin{aligned} \frac{d\mathbf{X}^{\varepsilon,p}}{dt} + \nabla\varphi_p(\mathbf{X}^{\varepsilon,p}) &= \mathbf{b}^X(\mathbf{X}^{\varepsilon,p}, \mathbf{Z}^{\varepsilon,p}) + \frac{1}{\varepsilon}\sigma\eta^{\varepsilon}, \quad \mathbf{X}_0^{\varepsilon,p} = \mathbf{x}_0, \\ \frac{d\mathbf{Z}^{\varepsilon,p}}{dt} + \nabla\psi_p(\mathbf{Z}^{\varepsilon,p}) &= \mathbf{b}^Z(\mathbf{X}^{\varepsilon,p}, \mathbf{Z}^{\varepsilon,p}), \quad \mathbf{Z}_0^{\varepsilon,p} = \mathbf{z}_0. \end{aligned} \right. \qquad (50) $$

It can be shown [4] that, if $\varphi$ satisfies the condition:

$$ \sup_{p \ge 1} \sup_{\mathbf{x} \in \mathbb{R}^n} \|\nabla \varphi_p(\mathbf{x})\| < \infty, \qquad (51) $$

where $\varphi_p$ is the Yosida approximation (45) of $\varphi$, then the sequence of solutions of (49) $\{\mathbf{X}^{\varepsilon,p}, p \ge 1\}$ is a Cauchy sequence in $C([0, T]; \mathbb{R}^n)$, the limit $\mathbf{X}^{\varepsilon}$ satisfies the differential inclusion (30) and its solution is unique.

A similar statement using the sequence of solutions of (50) $\{(\mathbf{X}^{\varepsilon,p}, \mathbf{Z}^{\varepsilon,p}), p \ge 1\}$ in $C([0, T]; \mathbb{R}^n \times \mathbb{R}^m)$, holds for the existence and uniqueness of a solution for (31) when $\varphi$ (but not necessarily $\psi$) satisfies the condition (51), while $\psi$ satisfies the assumption:

$$ \sup_{p \ge 1} \psi_p(\mathbf{z}_0) < \infty. \qquad (52) $$

For the convenience of the reader we give the proofs of these results in Appendix C.

## 5.2 Diffusion approximation for Equation (30)

We consider the $\mathbb{R}^n$-valued process $\mathbf{X}^{\varepsilon}$ solution of the multivalued ODE (30) when $\eta^{\varepsilon}$ is given by (2). We assume that $\mathbf{b}$ is Lipschitz and that $\varphi$ satisfies the condition (51). We also consider the limiting $\mathbb{R}^n$-valued process $\mathbf{X}^0$ solution of the multivalued SDE

$$ d\mathbf{X}^0 + \partial\varphi(\mathbf{X}^0)dt \ni \mathbf{b}(\mathbf{X}^0)dt + \Gamma d\mathbf{W}_t, \qquad (53) $$

driven by the same Brownian motion, with $\Gamma = \sigma\mathbf{A}^{-1}\mathbf{K}$. Existence and uniqueness of the solution of (53) is the same one as in Proposition C.1 and is discussed in Appendix D. The following proposition gives the convergence of the process $(\mathbf{X}^{\varepsilon} - \mathbf{X}^0)$ to zero. It is proved in Appendix E.

**Proposition 5.1.** 1. We have for all $p \ge 1$:

$$ \sup_{\varepsilon} \mathbb{E} \left[ \sup_{t \le T} \| \mathbf{X}_t^{\varepsilon,p} - \mathbf{X}_t^{\varepsilon} \|^2 \right] \le \frac{C_T}{p}, \qquad (54) $$

where $\mathbf{X}^{\varepsilon,p}$ is the approximation (49) of $\mathbf{X}^{\varepsilon}$.

2. The continuous process $(\mathbf{X}^{\varepsilon} - \mathbf{X}^0)$ converges in probability to zero as $\varepsilon \to 0$.

## 5.3 Diffusion approximation for Equation (31)

We consider the $\mathbb{R}^n \times \mathbb{R}^m$-valued process $(\mathbf{X}^{\varepsilon}, \mathbf{Z}^{\varepsilon})$ solution of the multivalued ODE (31) when $\eta^{\varepsilon}$ is given by (2). We assume that $\mathbf{b}^X$ and $\mathbf{b}^Z$ are Lipschitz, that $\varphi$ satisfies (51), and that $\psi$ satisfies (52). We also consider the limiting $\mathbb{R}^n \times \mathbb{R}^m$-valued process $(\mathbf{X}^0, \mathbf{Z}^0)$ solution of the multivalued SDE

$$ d\mathbf{X}^0 + \partial\varphi(\mathbf{X}^0)dt \ni \mathbf{b}^X(\mathbf{X}^0, \mathbf{Z}^0)dt + \Gamma d\mathbf{W}_t, \quad d\mathbf{Z}^0 + \partial\psi(\mathbf{Z}^0)dt \ni \mathbf{b}^Z(\mathbf{X}^0, \mathbf{Z}^0)dt, \qquad (55) $$

driven by the same Brownian motion, with $\Gamma = \sigma\mathbf{A}^{-1}\mathbf{K}$. Existence and uniqueness of the solution of (55) is the same one as in Proposition C.2. The following proposition gives the convergence of the continuous process $(\mathbf{X}^{\varepsilon} - \mathbf{X}^0, \mathbf{Z}^{\varepsilon} - \mathbf{Z}^0)$ to zero. It is proved in Appendix F.
**Proposition 5.2.** 1. For all $p \ge 1$, we have

$$ \sup_{\epsilon} \mathbb{E} \left[ \sup_{t \le T} \left\{ \| \mathbf{X}_t^{\epsilon, p} - \mathbf{X}_t^{\epsilon} \|^2 + \| \mathbf{Z}_t^{\epsilon, p} - \mathbf{Z}_t^{\epsilon} \|^2 \right\} \right] \le \frac{C_T}{p}, \qquad (56) $$

where $(\mathbf{X}^{\epsilon,p}, \mathbf{Z}^{\epsilon,p})$ is the approximation (50) of $(\mathbf{X}^{\epsilon}, \mathbf{Z}^{\epsilon})$.

2. The continuous process $(\mathbf{X}^{\epsilon} - \mathbf{X}^0, \mathbf{Z}^{\epsilon} - \mathbf{Z}^0)$ converges in probability to zero as $\epsilon \to 0$.

# 6 Control variate method in the multivalued case

We here consider the multivalued case. Let $\mathbf{X}^{\epsilon}$ satisfy (30) or $(\mathbf{X}^{\epsilon}, \mathbf{Z}^{\epsilon})$ satisfy (31). We want to estimate $I^{\epsilon}$ defined by (3) when $\mathbf{X}^{\epsilon}$ satisfies (30) or $I^{\epsilon} = \mathbb{E}[F(\mathbf{X}^{\epsilon}, \mathbf{Z}^{\epsilon})]$ when $(\mathbf{X}^{\epsilon}, \mathbf{Z}^{\epsilon})$ satisfies (31). The control variate method can be applied in this framework as in the ODE case addressed in Section 2. The control variate estimator $\hat{J}_N^\epsilon$ and the optimal control variate estimator $\hat{K}_N^\epsilon$ are defined by (15) and (25), respectively, for $\mathbf{X}^{\epsilon}$ satisfying (30), they are asymptotically normal and their asymptotic variances are (18) and (28), respectively. For $(\mathbf{X}^{\epsilon}, \mathbf{Z}^{\epsilon})$ satisfying (31) the control variate estimator

$$ \hat{J}_N^\epsilon \triangleq \frac{1}{N} \sum_{k=1}^{N} F(\mathbf{X}^{\epsilon}(\mathbf{W}^k), \mathbf{Z}^{\epsilon}(\mathbf{W}^k)) - F(\mathbf{X}^0(\mathbf{W}^k), \mathbf{Z}^0(\mathbf{W}^k)) + I^0, \qquad (57) $$

with $I^0 = \mathbb{E}[F(\mathbf{X}^0, \mathbf{Z}^0)]$, is asymptotically normal with an asymptotic variance given by

$$ \sigma_{J^{\epsilon}}^2 = \mathrm{Var}(F(\mathbf{X}^{\epsilon}, \mathbf{Z}^{\epsilon}) - F(\mathbf{X}^0, \mathbf{Z}^0)). \qquad (58) $$

The practical optimal control variate estimator is

$$ \hat{K}_N^\epsilon \triangleq \hat{\rho}_N^\epsilon I^0 + \frac{1}{N} \sum_{k=1}^{N} F(\mathbf{X}^{\epsilon}(\mathbf{W}^k), \mathbf{Z}^{\epsilon}(\mathbf{W}^k)) - \hat{\rho}_N^\epsilon F(\mathbf{X}^0(\mathbf{W}^k), \mathbf{Z}^0(\mathbf{W}^k)), \qquad (59) $$

where $\hat{\rho}_N^\epsilon$ is the empirical correlation

$$ \hat{\rho}_N^\epsilon = \frac{\sum_{k=1}^{N} (F(\mathbf{X}^{\epsilon}(\mathbf{W}^k), \mathbf{Z}^{\epsilon}(\mathbf{W}^k)) - \hat{I}_N^\epsilon)(F(\mathbf{X}^0(\mathbf{W}^k), \mathbf{Z}^0(\mathbf{W}^k)) - \hat{I}_N^0)}{\sum_{k=1}^{N} (F(\mathbf{X}^0(\mathbf{W}^k), \mathbf{Z}^0(\mathbf{W}^k)) - \hat{I}_N^0)^2}, \qquad (60) $$

$\hat{I}_N^\epsilon$ is the standard Monte Carlo estimator

$$ \hat{I}_N^\epsilon \triangleq \frac{1}{N} \sum_{k=1}^{N} F(\mathbf{X}^{\epsilon}(\mathbf{W}^k), \mathbf{Z}^{\epsilon}(\mathbf{W}^k)), \qquad (61) $$

and $\hat{I}_N^0 = \frac{1}{N} \sum_{k=1}^{N} F(\mathbf{X}^0(\mathbf{W}^k), \mathbf{Z}^0(\mathbf{W}^k))$. The estimator $\hat{K}_N^\epsilon$ is asymptotically normal with an asymptotic variance given by

$$ \sigma_{K^{\epsilon}}^2 = \mathrm{Var}(F(\mathbf{X}^{\epsilon}, \mathbf{Z}^{\epsilon}) - \rho^{\epsilon} F(\mathbf{X}^0, \mathbf{Z}^0)), \qquad (62) $$

with $\rho^{\epsilon} = \mathrm{Cov}(F(\mathbf{X}^{\epsilon}, \mathbf{Z}^{\epsilon}), F(\mathbf{X}^0, \mathbf{Z}^0))/\mathrm{Var}(F(\mathbf{X}^0, \mathbf{Z}^0))$. The asymptotic variance of the estimator $\hat{I}_N^\epsilon$ is

$$ \sigma_{I^{\epsilon}}^2 = \mathrm{Var}(F(\mathbf{X}^{\epsilon}, \mathbf{Z}^{\epsilon})). \qquad (63) $$
For small $\epsilon$ the asymptotic variance $\sigma_{J\epsilon}^2$ is approximately equal to $\text{Var}(F(\mathbf{X}^0, \mathbf{Z}^0))$ and the asymptotic variances $\sigma_{J\epsilon}^2$ and $\sigma_{K\epsilon}^2$ are small, by Propositions 5.1 and 5.2. More quantitatively, if $F(\mathbf{X}, \mathbf{Z}) = f(\mathbf{X}_T, \mathbf{Z}_T)$ and $f$ is a smooth function with bounded derivatives, then the asymptotic variances $\sigma_{J\epsilon}^2$ and $\sigma_{K\epsilon}^2$ are of order $\epsilon^2$ when $\mathbf{X}^\epsilon$ satisfies (30):

$$ \sigma_{K\epsilon}^2 \le C\epsilon^2, \quad 0 \le \sigma_{J\epsilon}^2 - \sigma_{K\epsilon}^2 \le C\epsilon^4, \qquad (64) $$

or of order $\epsilon$ when $(\mathbf{X}^\epsilon, \mathbf{Z}^\epsilon)$ satisfies (31):

$$ \sigma_{K\epsilon}^2 \le C\epsilon, \quad 0 \le \sigma_{J\epsilon}^2 - \sigma_{K\epsilon}^2 \le C\epsilon^2. \qquad (65) $$

Eqs. (64-65) are consequences of the following lemma proved in Appendix G.

**Lemma 6.1.** 1. Let $\mathbf{X}^\epsilon$ satisfy (30). If $f,g: \mathbb{R}^n \to \mathbb{R}$ are smooth functions with bounded derivatives and $T > 0$, then there exists $C > 0$ such that, for any $t \in [\epsilon, T]$,

$$ |\mathbb{E}[g(\mathbf{Z}_t^0)(f(\mathbf{X}_t^\epsilon) - f(\mathbf{X}_t^0))]| \le C\epsilon^2, \quad \mathbb{E}[(f(\mathbf{X}_t^\epsilon) - f(\mathbf{X}_t^0))^2] \le C\epsilon^2. \qquad (66) $$

2. Let $(\mathbf{X}^\epsilon, \mathbf{Z}^\epsilon)$ satisfy (31). If $f,g: \mathbb{R}^{n+m} \to \mathbb{R}$ are smooth functions with bounded derivatives and $T > 0$, then there exists $C > 0$ such that, for any $t \in [\epsilon, T]$,

$$ |\mathbb{E}[g(\mathbf{X}_t^0, \mathbf{Z}_t^0)(f(\mathbf{X}_t^\epsilon, \mathbf{Z}_t^\epsilon) - f(\mathbf{X}_t^0, \mathbf{Z}_t^0))]| \le C\epsilon, \quad \mathbb{E}[(f(\mathbf{X}_t^\epsilon, \mathbf{Z}_t^\epsilon) - f(\mathbf{X}_t^0, \mathbf{Z}_t^0))^2] \le C\epsilon. \qquad (67) $$

# 7 Numerical simulations in the multi-valued case

We present examples which are non-smooth dynamical systems that are prevalent in engineering mechanics. Examples 7.1-7.2 are oscillators involving friction or/and elasto-plastic behaviours, Exemple 7.3 is a nonlinear and nonsmooth two-degree of freedom (TDOF) system, they can be described by Eqs. (30) and (31). Examples 7.4-7.5 which do not fall within the scope of any aforementioned case correspond to an obstacle problem and to the reflection of the integral of a colored noise, respectively.

## 7.1 Non-smooth systems in the form of Equations (30) and (31)

**Example 7.1 (friction behaviour).** With Equation (30) in mind, we take $\forall x \in \mathbb{R}$, $\varphi(x) \triangleq c_f|x|$ where $c_f > 0$ is a friction coefficient. The $\mathbb{R}$-valued process $\mathbf{X}^\epsilon$ represents the velocity of a material point (stick-slip motion) subjected to friction and colored noise. See for instance [29] for an explanation of the physics behind and [3] for the use of SDEs with multivalued drift for modeling. As $\epsilon \to 0$, $\mathbf{X}^\epsilon \to \mathbf{X}^0$ where $\mathbf{X}^0$ satisfies Equation (53). For the stochastic simulation, we proceed as follows:

$$ \bullet \hat{X}_0^{\epsilon,k} = x_0 \text{ and for } 0 \le n \le N_T - 1, $$

$$ \hat{X}_{n+1}^{\epsilon,k} = \hat{X}_{n}^{\epsilon,k} + \frac{\delta t}{\epsilon} \hat{\eta}_{n}^{\epsilon,k} - \delta t \text{proj}_{[-c_f, c_f]} \left( \frac{\hat{X}_{n}^{\epsilon,k}}{\delta t} + \frac{1}{\epsilon} \hat{\eta}_{n}^{\epsilon,k} \right), $$

with $\hat{\eta}_{n}^{\epsilon,k}$ described in Subsection 4.1.
• $\hat{X}_0^{0,k} = x_0$ and for $0 \le n \le N_T - 1$,

$$ \hat{X}_{n+1}^{0,k} = \hat{X}_n^{0,k} + (\delta t)^{1/2} C \Delta W_n^k - \delta t \text{proj}_{[-c_f, c_f]} \left( \frac{\hat{X}_n^{0,k}}{\delta t} + \frac{C}{\sqrt{\delta t}} \Delta W_n^k \right). $$

We are interested in $\mathbb{E}[(X_T^\varepsilon)^2]$ for $T = 1$. The expectation of the control variate is $\mathbb{E}[(X_T^0)^2]$. The latter can be represented as $c(x_0, 0)$ where $c$ satisfies the following backward in time partial differential inclusion

$$ \left\{ \begin{array}{l} \partial_t c(x,t) + \frac{C^2}{2} \partial_x^2 c(x,t) \in \partial \varphi(x) \partial_x c(x,t) = 0, \quad \text{for } (x,t) \in \mathbb{R} \times [0,1), \\ c(x,1) = x^2, \quad \text{for } x \in \mathbb{R}. \end{array} \right. \qquad (68) $$

It can be estimated by solving this partial differential inclusion with a finite difference method. We proceed as follows. For every $t > 0$, the function $x \mapsto c(x,t)$ is smooth and even, provided that the initial condition is smooth and even. Indeed, this comes from the probabilistic representation and the fact that, for any starting point $x \in \mathbb{R}$, $\{X_t^x, t \ge 0\}$ and $\{X_t^{-x}, t \ge 0\}$ have the same distribution because $\varphi$ is even. Therefore we must have $\forall t > 0$, $\partial_x c(0,t) = 0$. The solution of (68) is thus estimated by applying a finite difference method to

$$ \left\{ \begin{array}{l} \partial_t c(x,t) + \frac{C^2}{2} \partial_x^2 c(x,t) - c_f \partial_x c(x,t) = 0, \quad \text{for } (x,t) \in (0, \infty) \times [0,1), \\ \partial_x c(0,t) = 0, \quad \text{for } t \in [0,1), \\ c(x,1) = x^2, \quad \text{for } x \in [0, +\infty). \end{array} \right. \qquad (69) $$

The whole function $x \mapsto c(x,t)$ can be recovered by using the symmetry property.

**Example 7.2 (elasto-plastic behaviour).** With Eq. (31), we consider $\varphi \triangleq 0$, $\psi \triangleq \chi_D$ the indicator function of $D \triangleq [-c_{ep}, c_{ep}]$ in the sense of convex analysis, that is $\chi_D(x) = 0$ if $x \in D$ and $+\infty$ otherwise. Here $c_{ep} > 0$ is an elasto-plastic coefficient. The real-valued process $X^\varepsilon$ represents the velocity of a material point subjected to an elasto-plastic restoring force and colored noise. The process $Z^\varepsilon$ taking values in $[-c_{ep}, c_{ep}]$ represents the restoring force. See for instance [27] for an explanation of the physics and the use of SDEs with multivalued drift for modeling. Here we are interested in $\mathbb{E}[(X_T^\varepsilon)^2 + (Z_T^\varepsilon)^2]$ and in $\mathbb{P}(|Z_T^\varepsilon| = c_{ep})$ for $T = 1$. As $\varepsilon \to 0$, $(X^\varepsilon, Z^\varepsilon) \to (X^0, Z^0)$ where $(X^0, Z^0)$ satisfies (55). For the stochastic simulation, we proceed as follows:

• $\hat{X}_0^{\varepsilon,k} = x_0$ and $\hat{Z}_0^{\varepsilon,k} = z_0$ and for $0 \le n \le N_T - 1$,

$$ \left\{ \begin{array}{l} \hat{Z}_{n+1}^{\varepsilon,k} = \text{proj}_{[-c_{ep}, c_{ep}]} (\hat{Z}_n^{\varepsilon,k} + \delta t \hat{X}_n^{\varepsilon,k}), \\ \hat{X}_{n+1}^{\varepsilon,k} = \hat{X}_n^{\varepsilon,k} - \delta t \hat{Z}_n^{\varepsilon,k} + \frac{\delta t}{\varepsilon} \hat{\eta}_n^{\varepsilon,k}. \end{array} \right. $$

• $\hat{X}_0^{0,k} = x_0$ and $\hat{Z}_0^{0,k} = z_0$ and for $0 \le n \le N_T - 1$,

$$ \left\{ \begin{array}{l} \hat{Z}_{n+1}^{0,k} = \text{proj}_{[-c_{ep}, c_{ep}]} (\hat{Z}_n^{0,k} + \delta t \hat{X}_n^{0,k}), \\ \hat{X}_{n+1}^{0,k} = \hat{X}_n^{0,k} - \delta t \hat{Z}_n^{0,k} + \sqrt{\delta t} C \Delta W_n^k. \end{array} \right. $$

The expectations of the control variates are $\mathbb{E}[(X_T^0)^2 + (Z_T^0)^2]$ and $\mathbb{P}(|Z_T^0| = c_{ep})$ for $T = 1$. They are estimated using the PDE method of [28].
**Example 7.3 (nonlinear and nonsmooth TDOF).** A TDOF with an elasto-plastic element can be represented as a system of the form (30) which becomes

$$ \left\{ \begin{array}{l} \frac{dX_1^\epsilon}{dt} + \partial\chi_D(X_1^\epsilon) = X_2^\epsilon, \quad \frac{dX_2^\epsilon}{dt} = -\tilde{g}_2(X_1^\epsilon, X_2^\epsilon, X_3^\epsilon, X_4^\epsilon) + \frac{1}{\epsilon}\eta^\epsilon, \\ \\ \frac{dX_3^\epsilon}{dt} = X_4^\epsilon, \quad \frac{dX_4^\epsilon}{dt} = -\tilde{g}_4(X_1^\epsilon, X_2^\epsilon, X_3^\epsilon, X_4^\epsilon). \end{array} \right. \qquad (70) $$

Here $D = [-c_p, c_p]$, $c_p = 0.25$. In addition to $\mathbb{E}[(X_{1,T}^\epsilon)^2 + (X_{3,T}^\epsilon)^2]$ for $T = 1$, we are interested in $\mathbb{P}(|X_{1,T}^\epsilon| \le a, |X_{3,T}^\epsilon| \le b)$. As $\epsilon \to 0$, $X^\epsilon \to X^0$ where

$$ \left\{ \begin{array}{l} dX_1^0 + \partial\chi_D(X_1^0)dt = X_2^0dt, \quad dX_2^0 = -\tilde{g}_2(X_1^0, X_2^0, X_3^0, X_4^0)dt + CdW, \\ \\ dX_3^0 = X_4^0dt, \quad dX_4^0 = -\tilde{g}_4(X_1^0, X_2^0, X_3^0, X_4^0)dt. \end{array} \right. \qquad (71) $$

Here, for simplicity the other elements are linear,

$$ \tilde{g}_2(x_1, x_2, x_3, x_4) \triangleq k_1x_1 + c_1x_2 - k_3(x_3 - x_1) - c_3(x_4 - x_2) $$

and

$$ \tilde{g}_4(x_1, x_2, x_3, x_4) \triangleq k_3(x_3 - x_1) + c_3(x_4 - x_2). $$

The simulation of (70) and (71) is similar to what is explained above. We take $c_1 = c_3 = k_1 = k_3 = 1$.

## 7.2 Non-smooth systems : beyond Equations (30) and (31)

The two models presented in this subsection do not fall in the scope of our theoretical results, though they are not too far off. The presentation of the impact problem remains formal. The behaviour of the control variate estimator is investigated via numerical experiments.

**Example 7.4 (impact problem).** The pair displacement-velocity $X^\epsilon = (X_1^\epsilon, X_2^\epsilon)$ (taking values in $\mathbb{R}^2$) of a colored noise driven oscillator constrained by an obstacle can be formulated in terms of an equation of the form (36) when $|X_{1,t}^\epsilon| < P_O$ with the condition (that expresses the switch of the velocity at collision): for all $t$, $|X_{1,t}^\epsilon| = P_O \implies X_{2,t+}^\epsilon = -\epsilon X_{2,t-}^\epsilon$ where $P_O$ is the location of the obstacle and $\epsilon \in [0, 1]$ is the coefficient of restitution of energy. The notations $X_{2,t\pm}^\epsilon$ stand for the velocity immediately before and after the collision. Here we are interested in $\mathbb{E}[(X_{2,T}^\epsilon)^2]$ for $T = 1$. Formally, as $\epsilon \to 0$, the $\mathbb{R}^2$-valued limit process $X^0 = (X_1^0, X_2^0)$ is a white noise driven oscillator constrained by an obstacle that can be formulated similarly to the former case, except that we replace (36) by (37). When $\epsilon = 1$ (resp. $0 \le \epsilon < 1$), we say that the collisions are elastic (resp. inelastic). It is important to stress that obstacle problems with inelastic collisions deserve more attention for practical purposes since in real world phenomena kinetic energy is dissipated through heat or plastic deformation. With elastic collisions, there is no loss of kinetic energy. For the stochastic simulation, we use the same numerical procedure as for (36) and (37), except that if we find out that the $(n+1)^{st}$ point does not satisfy the obstacle condition, i.e. $|\hat{X}_{1,n+1}^\epsilon| > P_O$, we adjust the time step length to $\theta_{n+1}\delta t$ with $\theta_{n+1} \triangleq \frac{\pm P_O - \hat{X}_{1,n}^\epsilon}{\hat{X}_{1,n+1}^\epsilon - \hat{X}_{1,n}^\epsilon}$ and set $t_{n+1} \triangleq t_n + \theta_{n+1}\delta t$, $\hat{X}_{1,n+1}^\epsilon \triangleq P_O$, $\hat{X}_{2,n+1}^\epsilon \triangleq -\epsilon(\hat{X}_{2,n}^\epsilon - \theta_{n+1}\delta t f(\hat{X}_{1,n}^\epsilon, \hat{X}_{2,n}^\epsilon) + \theta_{n+1}\frac{\delta t}{\epsilon}\hat{\eta}_n^\epsilon)$ and $\hat{\eta}_{2,n+1}^\epsilon \triangleq \hat{\eta}_n^\epsilon(1 - \theta_{n+1}\delta t \frac{A}{\epsilon^2}) + P_O$.
$\sqrt{\theta_{n+1}\delta t \frac{K_{\text{ou}}}{\varepsilon}} \Delta W_n^k$. A similar adjustment is done in the other cases with Langevin and white noises. The expectation of the control variate is $\mathbb{E}[(X_{2,T}^0)^2]$ for $T=1$ which is estimated using the PDE method of [28].

**Example 7.5 (reflection of an integrated colored noise).** Define $E \triangleq [0, \infty)$ and consider the indicator function of $E$, that is $\chi_E(x) = 0$ if $x \in E$ and $+\infty$ otherwise. The reflection of an integrated colored noise corresponds to the case where $X^\varepsilon$ satisfies

$$ \frac{dX^\varepsilon}{dt} + \partial\chi_E(X^\varepsilon) \ni \frac{1}{\varepsilon}\eta^\varepsilon, \qquad (72) $$

and $X^0$, the limit process as $\varepsilon \to 0$, is a reflected Brownian motion

$$ dX^0 + \partial\chi_E(X^0)dt \ni C dW. \qquad (73) $$

We are interested in $\mathbb{E}[X_T^\varepsilon]$ for $T=1$. For the stochastic simulation of (72) and (73), we use the following scheme: $\hat{X}_0^{\varepsilon,k} = x_0$, $\hat{X}^{0,k} = x_0$ and for $0 \le n \le N_T - 1$,

$$ \bullet \hat{X}_{n+1}^{\varepsilon,k} = \text{proj}_E \left( \hat{X}_n^{\varepsilon,k} + \frac{\delta t}{\varepsilon} \hat{\eta}_n^{\varepsilon,k} \right), $$

$$ \bullet \hat{X}_{n+1}^{0,k} = \text{proj}_E \left( \hat{X}_n^{0,k} + \sqrt{\delta t} C \Delta W_n^k \right). $$

The expectation of the control variate is given by an explicit formula $\mathbb{E}[X_1^0] = \sqrt{2/\pi}$. Indeed, the backward Kolmogorov equation for the reflected Brownian motion in (73) is

$$ \partial_t w = C \partial_x^2 w, \quad x > 0, \ t > 0, \quad w(x, t=0) = x, \quad x > 0, \quad \partial_x w(0, t) = 0, \ t > 0. $$

It has an explicit solution

$$ w(x,t) = \frac{1}{\sqrt{4C\pi t}} \int_0^\infty y \left( \exp\left(-\frac{(x-y)^2}{4Ct}\right) + \exp\left(-\frac{(x+y)^2}{4Ct}\right) \right) dy, $$

which gives $\mathbb{E}[X_T^0] = w(0,T) = \sqrt{2/\pi}$ for $T=1$. In this case, we can provide an ad hoc proof to get an estimate similar to (35) (see Appendix H):

$$ \mathbb{E}[(X_T^\varepsilon - X_T^0)^2] \le C\varepsilon^2 |\log \varepsilon|. \qquad (74) $$

The $\log \varepsilon$ correction comes from a maximal inequality for the OU process [14] and a standard result on the maxima of Gaussian processes [26].

## 7.3 Numerical experiments

We report our numerical results for the four systems mentioned above. The convention is as in Subsection 4.3. In each of the four figures below, there are four subfigures (a)-(b)-(c)-(d). For subfigures (a) and (b), the driving force is an Ornstein-Uhlenbeck noise (32) with $A = K = 1$. In subfigure (a), the dashed black lines, the dotted blue lines, and the solid red lines represent the standard MC estimator $\hat{I}_N^\varepsilon$ and the control variate estimators $\hat{J}_N^\varepsilon$ and $\hat{K}_N^\varepsilon$, respectively. The dotted black line represents the expectation of the control variate $I^0$. The objective of the subfigure (b)
is to illustrate the bounds (64) and (65). The same description applies to (c) and (d), except they correspond to the case of a Langevin noise (33-34) with $\mu = \gamma = K = 1$. For the examples in which $\mathbf{X}^\epsilon$ satisfies (30), the asymptotic variance the standard MC estimator $\hat{I}_N^\epsilon$ is estimated by (42), the asymptotic variance of the control variate estimator $\hat{J}_N^\epsilon$ is estimated by (43), and the asymptotic variance of the optimal control variate estimator $\hat{K}_N^\epsilon$ is estimated by (44). For the examples in which $(\mathbf{X}^\epsilon, \mathbf{Z}^\epsilon)$ satisfies (31), the asymptotic variance of the standard MC estimator $\hat{I}_N^\epsilon$ is estimated by

$$ \hat{\sigma}_{I^\epsilon, N}^2 = \frac{1}{N} \sum_{k=1}^{N} \left( F(\mathbf{X}^\epsilon(\mathbf{W}^k), \mathbf{Z}^\epsilon(\mathbf{W}^k)) \right)^2 - (\hat{I}_N^\epsilon)^2, \quad (75) $$

the asymptotic variance of the control variate estimator $\hat{J}_N^\epsilon$ is estimated by

$$ \hat{\sigma}_{J^\epsilon, N}^2 = \frac{1}{N} \sum_{k=1}^{N} \left( F(\mathbf{X}^\epsilon(\mathbf{W}^k), \mathbf{Z}^\epsilon(\mathbf{W}^k)) - F(\mathbf{X}^0(\mathbf{W}^k), \mathbf{Z}^0(\mathbf{W}^k)) + I^0 \right)^2 - (\hat{J}_N^\epsilon)^2, \quad (76) $$

and the asymptotic variance of the optimal control variate estimator $\hat{K}_N^\epsilon$ is estimated by

$$ \hat{\sigma}_{K^\epsilon, N}^2 = \frac{1}{N} \sum_{k=1}^{N} \left( F(\mathbf{X}^\epsilon(\mathbf{W}^k), \mathbf{Z}^\epsilon(\mathbf{W}^k)) - \hat{\rho}_N^\epsilon F(\mathbf{X}^0(\mathbf{W}^k), \mathbf{Z}^0(\mathbf{W}^k)) + \hat{\rho}_N^\epsilon I^0 \right)^2 - (\hat{K}_N^\epsilon)^2, \quad (77) $$

with $\hat{\rho}_N^\epsilon$ defined by (60). $\hat{\sigma}_{I^\epsilon, N}^2$, $\hat{\sigma}_{J^\epsilon, N}^2$, and $\hat{\sigma}_{K^\epsilon, N}^2$ are consistent estimators of $\sigma_{I^\epsilon}^2$, $\sigma_{J^\epsilon}^2$, and $\sigma_{K^\epsilon}^2$, respectively.

Similarly to what was presented in Section 4, we use $N = 10^4$ samples with a time step of $\delta t = 10^{-5}$. In Figures 6 and 7, we report the numerical results for the friction and elasto-plastic problems, which are of the form (30) and (31), respectively. In Figures 8 and 9, we report the numerical results for the obstacle problem and for the reflection of the integral of a colored noise. The numerical results include errors bars on the estimators for each value of $\epsilon \in \{0.1, 0.5, 0.9\}$ in the $\hat{I}_N^\epsilon$, $\hat{J}_N^\epsilon$, $\hat{K}_N^\epsilon$ order.

The theoretical predictions provided by (64) and (65) are based on the condition that $f$ has bounded derivatives. The assumption that $f$ is smooth is important but the hypothesis on the boundedness of the derivatives can certainly be relaxed. The numerical results shown in Figures 6 and 7 are in good agreement with the theoretical predictions: the asymptotic variances $\sigma_{J^\epsilon}^2$ and $\sigma_{K^\epsilon}^2$ behave as $O(\epsilon^2)$. The only cases where the behavior is $O(\epsilon)$, and not $O(\epsilon^2)$, are when the quantity of interest is of the form $\mathbb{E}[f(\mathbf{X}_T^\epsilon)]$ or $\mathbb{E}[f(\mathbf{X}_T^\epsilon, \mathbf{Z}_T^\epsilon)]$ with a function $f$ that is not smooth, which is not surprising. In Figure 7, we also observe that $\sigma_{J^\epsilon}^2$ and $\sigma_{K^\epsilon}^2$ behave as $O(\epsilon^2)$ which is better than the behaviour $O(\epsilon)$ expected from (65) (which is an upper bound). In Figures 8 and 9, the numerical results concern two problems which do not fall within the scope of our theoretical predictions. The first one (Figure 8) is the impact problem that cannot be formulated in the form a differential inclusion of the form (30) or (31). The function $f$ is smooth but the behavior of $\sigma_{J^\epsilon}^2$ and $\sigma_{K^\epsilon}^2$ is not of order $O(\epsilon^2)$, only of order $O(\epsilon)$. The second one (Figure 9) is the reflection of an integrated colored noise that can be formulated with a differential inclusion which is similar to (30) but the multivalued drift does not satisfy the condition (51). However, $\sigma_{J^\epsilon}^2$ and $\sigma_{K^\epsilon}^2$ behave as $O(\epsilon^2)$. To summarize, the numerical simulations indicate that the $O(\epsilon^2)$ behavior of the asymptotic variance is observed in the cases predicted by the theory and also slightly beyond. The smoothness of the function $f$ that appears in the quantity of interest is, however, an important condition to ensure the $O(\epsilon^2)$-behavior, otherwise one only observes a $O(\epsilon)$-behavior.
<figure><img src="image_7.png" /><figcaption>Figure 6: Example 7.1 (friction problem). The target is to estimate $I^\epsilon = \mathbb{E}[(X_T^\epsilon)^2]$ for $T = 1$ where $X^\epsilon$ satisfies (30) with $\varphi(x) \triangleq c_f|x|$ with $c_f > 0$. The expectation of the control variate $I^0 = \mathbb{E}[(X_T^0)^2]$, where $X^0$ satisfies (53), is obtained by solving the partial differential inclusion (68).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 7: Example 7.2 (elasto-plastic problem). In the top row the target is to estimate $I^\epsilon = \mathbb{E}[(X_T^\epsilon)^2 + (Z_T^\epsilon)^2]$ with $T = 1$ where $(X^\epsilon, Z^\epsilon)$ satisfies (31) with $\varphi(x) \triangleq 0$ and $\psi(x) \triangleq 0$ if $|x| \le c_{ep}$ and $\infty$ otherwise. Here $c_{ep} = 0.25$. The expectation of the control variate $I^0 = \mathbb{E}[(X_T^0)^2 + (Z_T^0)^2]$, where $(X^0, Z^0)$ satisfies (55), is obtained by using the PDE method of [28]. In the bottom row the target is to estimate $I^\epsilon = \mathbb{P}(|Z_T^\epsilon| = c_{ep})$.</figcaption></figure>

# 8 Concluding remark

When the expectation of the limit process $\mathbb{E}[f(X_T^0)]$ cannot be computed by a PDE method but is estimated by a massive Monte Carlo method, the control variate method with $\hat{J}_N^\epsilon$ (or $\hat{K}_N^\epsilon$) shares an important similarity with a two-level Monte Carlo method [10] in the sense that massive but cheap simulations are performed (samples of $X^0$ used to estimate $\mathbb{E}[f(X_T^0)]$) together with a few expensive simulations (samples of $(X^\epsilon, X^0)$ used to estimate $\mathbb{E}[f(X_T^\epsilon) - f(X_T^0)]$). More generally, Multilevel Monte Carlo (MLMC) methods rely on random samples taken on different levels of accuracy, when several approximations with different costs and accuracies are available. The overall idea of MLMC methods is to reduce the computational cost of standard Monte Carlo methods by taking most samples with a low accuracy and corresponding low cost, and by taking only few samples with a high accuracy and corresponding high cost [10].
<figure><img src="image_9.png" /><figcaption>Figure 8: Example 7.4 (impact problem). The target is to estimate $I^{\epsilon} = \mathbb{E}[(X_{2,T}^{\epsilon})^2]$ with $T = 1$ where $\mathbf{X}^{\epsilon}$ satisfies the impact problem with a colored noise forcing. Here $P_O = 0.25$. The expectation of the control variate $I^0 = \mathbb{E}[(X_{2,T}^0)^2]$, where $\mathbf{X}^0$ satisfies the impact problem with a colored noise forcing, is obtained by using the PDE method of [28].</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 9: Example 7.5 (reflection of an integrated colored noise). The target is to estimate $I^{\epsilon} = \mathbb{E}[X_T^{\epsilon}]$ for $T = 1$ where $X^{\epsilon}$ satisfies (72). The expectation of the control variate is $I^0 = \mathbb{E}[X_T^0] = \sqrt{2/\pi}$.</figcaption></figure>

In this two-level Monte Carlo framework, the total cost of computing the control variate estimator $\hat{J}_{N}^{\epsilon}$ (or $\hat{K}_{N}^{\epsilon}$) is $C_{J}^{\epsilon} = N_{0}C_{0} + N_{1}C_{1}$ where $C_{0}$ is the cost of computing one realization of $f(\mathbf{X}_{T}^{0})$, $C_{1}$ is the cost of computing one realization of $(f(\mathbf{X}_{T}^{\epsilon}) - f(\mathbf{X}_{T}^{0}))$, $N_{0}$ is the number of samples of $f(\mathbf{X}_{T}^{0})$, and $N_{1}$ is the number of samples of $(f(\mathbf{X}_{T}^{\epsilon}) - f(\mathbf{X}_{T}^{0}))$. The variance of the control variate estimator $\hat{J}_{N}^{\epsilon}$ is $V_{J}^{\epsilon} = N_{0}^{-1}V_{0} + N_{1}^{-1}V_{1}$, where $V_{0}$ is the variance of $f(\mathbf{X}_{T}^{0})$ and $V_{1}$ is the variance of $(f(\mathbf{X}_{T}^{\epsilon}) - f(\mathbf{X}_{T}^{0}))$. For a fixed total budget $C_{\text{tot}}$, the variance is minimized when $(N_{1}/N_{0})^{2} = (V_{1}/V_{0})(C_{0}/C_{1})$ and it is then equal to $V_{J}^{\epsilon} = (\sqrt{V_{0}C_{0}} + \sqrt{V_{1}C_{1}})^{2}/C_{\text{tot}}$. This can be compared to the brute force Monte Carlo method: the cost is $C_{I}^{\epsilon} = N_{I}C_{1}$ where $N_{I}$ is the number of samples of $f(\mathbf{X}_{T}^{\epsilon})$ (we neglect the difference of cost between $(f(\mathbf{X}_{T}^{\epsilon}) - f(\mathbf{X}_{T}^{0}))$ and $f(\mathbf{X}_{T}^{\epsilon})$, which is very small because $\mathbf{X}^{\epsilon}$ is more difficult to simulate than $\mathbf{X}^{0}$) and the variance is $V_{I}^{\epsilon} = V_{0}/N_{I}$ (by (11)), so that for the total budget $C_{\text{tot}}$, we have $V_{I}^{\epsilon} = V_{0}C_{1}/C_{\text{tot}}$. If $V_{1}/V_{0} = O(\epsilon^{2})$ (by Proposition 2.1) and $C_{0}/C_{1} = O(\epsilon^{2})$ (because the time step used to simulate $\mathbf{X}^{\epsilon}$ should be $O(\epsilon^{2})$ smaller than the one used to simulate $\mathbf{X}^{0}$), then we find that $N_{1}/N_{0} = O(\epsilon^{2})$ (hence the “massive Monte Carlo” strategy for the control variate) and the ratio of the variance of the control variate estimator over the one of the brute force Monte Carlo estimator is finally $V_{J}^{\epsilon}/V_{I}^{\epsilon} = O(\epsilon^{2})$. The control variate method is very advantageous in this context.
A Proof of Proposition 3.5

Let us first study the driving noise. The process $\eta^1$ is a Gaussian, Markov process. It has the form

$$ \eta_t^1 = e^{-\mathbf{A}t} \eta_0^1 + \int_0^t e^{-\mathbf{A}(t-s)} \mathbf{K} d\mathbf{W}_s. $$

Its infinitesimal generator is:

$$ Q = \frac{1}{2} \sum_{k,k'=1}^{d} \sum_{k''=1}^{d'} K_{kk''} K_{k'k''} \partial_{\eta_k \eta_{k'}}^2 - \sum_{k,k'=1}^{d} A_{kk'} \eta_{k'} \partial_{\eta_k}. \quad (78) $$

The properties of the matrix **A** show that the process $\eta^1$ is stationary and ergodic; its unique invariant probability measure is the Gaussian measure with mean zero and variance **C** given by (8). The process $(\eta_t^\varepsilon)_{t\ge0}$ has the same distribution as $(\eta_{t/\varepsilon^2}^1)_{t\ge0}$ because $(\varepsilon^{-1}W_t)_{t\ge0}$ has the same distribution as $(W_{t/\varepsilon^2})_{t\ge0}$. Therefore it is a Markov process with generator $\varepsilon^{-2}Q$. The process $(\mathbf{X}^\varepsilon, \mathbf{X}^0, \boldsymbol{\eta}^\varepsilon)$ is Markov with generator $\mathcal{L}^\varepsilon$ given by:

$$ \mathcal{L}^{\varepsilon} = \frac{1}{\varepsilon^2} Q + \frac{1}{\varepsilon} \left[ \sum_{j=1}^{n} \sum_{i=1}^{d} \sigma_{ji}(\mathbf{x}) \eta_i \partial_{x_j} + \sum_{j=1}^{n} \sum_{i=1}^{d} \sum_{k=1}^{d'} \Gamma_{jk}(\mathbf{x}^0) K_{ik} \partial_{\eta_i x_j^0}^2 \right. \\ \left. + \left[ \sum_{j=1}^{n} b_j(\mathbf{x}) \partial_{x_j} + \frac{1}{2} \sum_{i,j=1}^{n} \sum_{k=1}^{d'} \Gamma_{ik}(\mathbf{x}^0) \Gamma_{jk}(\mathbf{x}^0) \partial_{x_i^0 x_j^0}^2 + \sum_{j=1}^{n} \tilde{b}_j(\mathbf{x}^0) \partial_{x_j^0} \right], \right] \quad (79) $$

where $Q$ is the generator (78).

**Lemma A.1.** For any smooth and bounded test function $\phi: \mathbb{R}^n \times \mathbb{R}^n \to \mathbb{R}$, for any compact subset $K$ of $\mathbb{R}^{2n}$, there exists a test function $\phi^\varepsilon$ such that

$$ \sup_{(\mathbf{x}, \mathbf{x}^0) \in K} |\phi^\varepsilon(\mathbf{x}, \mathbf{x}^0, \boldsymbol{\eta}) - \phi(\mathbf{x}, \mathbf{x}^0)| \le C\varepsilon(1 + \|\boldsymbol{\eta}\|^2), \quad (80) $$

$$ \sup_{(\mathbf{x}, \mathbf{x}^0) \in K} |\mathcal{L}^{\varepsilon} \phi^{\varepsilon}(\mathbf{x}, \mathbf{x}^0, \boldsymbol{\eta}) - \mathcal{L} \phi(\mathbf{x}, \mathbf{x}^0)| \le C\varepsilon(1 + \|\boldsymbol{\eta}\|^3), \quad (81) $$

for any $\varepsilon \in (0, 1)$, where $\mathcal{L}$ is the generator defined by

$$ \mathcal{L} = \sum_{j=1}^{n} \tilde{b}_j(\mathbf{x}) \partial_{x_j} + \sum_{j=1}^{n} \tilde{b}_j(\mathbf{x}^0) \partial_{x_j^0} + \frac{1}{2} \sum_{i,j=1}^{n} [\Gamma(\mathbf{x}) \Gamma(\mathbf{x})^T]_{ij} \partial_{x_i x_j}^2 \\ + \frac{1}{2} \sum_{i,j=1}^{n} [\Gamma(\mathbf{x}^0) \Gamma(\mathbf{x}^0)^T]_{ij} \partial_{x_i^0 x_j^0}^2 + \sum_{i,j=1}^{n} [\Gamma(\mathbf{x}) \Gamma(\mathbf{x}^0)^T]_{ij} \partial_{x_i x_j^0}^2. \quad (82) $$

*Proof.* Let $\phi(\mathbf{x}, \mathbf{x}^0)$ be a smooth and bounded test function. We look for a perturbed test function $\phi^\varepsilon$ of the form

$$ \phi^{\varepsilon}(\mathbf{x}, \mathbf{x}^0, \boldsymbol{\eta}) = \phi(\mathbf{x}, \mathbf{x}^0) + \varepsilon\phi_1(\mathbf{x}, \mathbf{x}^0, \boldsymbol{\eta}) + \varepsilon^2\phi_2(\mathbf{x}, \mathbf{x}^0, \boldsymbol{\eta}). \quad (83) $$
Applying $\mathcal{L}^\epsilon$ (given by (79)) to this $\phi^\epsilon$ we get

$$ \mathcal{L}^{\epsilon} \phi^{\epsilon} = \frac{1}{\epsilon} \left[ Q \phi_1(\mathbf{x}, \mathbf{x}^0, \boldsymbol{\eta}) + \sum_{j=1}^{n} \sum_{i=1}^{d} \sigma_{ji}(\mathbf{x}) \eta_i \partial_{x_j} \phi(\mathbf{x}, \mathbf{x}^0) \right] \\ + \left[ Q \phi_2(\mathbf{x}, \mathbf{x}^0, \boldsymbol{\eta}) + \sum_{j=1}^{n} \sum_{i=1}^{d} \sigma_{ji}(\mathbf{x}) \eta_i \partial_{x_j} \phi_1(\mathbf{x}, \mathbf{x}^0, \boldsymbol{\eta}) + \sum_{j=1}^{n} \sum_{i=1}^{d} \sum_{k=1}^{d'} \Gamma_{jk}(\mathbf{x}^0) K_{ik} \partial_{\eta_i x_j^0}^2 \phi_1(\mathbf{x}, \mathbf{x}^0, \boldsymbol{\eta}) \right] \\ + \left[ \sum_{j=1}^{d} b_j(\mathbf{x}) \partial_{x_j} \phi(\mathbf{x}, \mathbf{x}^0) + \sum_{j=1}^{d} \tilde{b}_j(\mathbf{x}^0) \partial_{x_j^0} \phi(\mathbf{x}, \mathbf{x}^0) + \frac{1}{2} \sum_{i,j=1}^{n} \sum_{k=1}^{d'} \Gamma_{ik}(\mathbf{x}^0) \Gamma_{jk}(\mathbf{x}^0) \partial_{x_i^0 x_j^0}^2 \phi(\mathbf{x}, \mathbf{x}^0) \right] \\ + O(\epsilon). \qquad (84) $$

The term $O(\epsilon)$ depends on the first-order derivatives of $\phi_1$ and $\phi_2$ with respect to $\mathbf{x}$ and on the first- and second-order derivatives of $\phi_1$ and $\phi_2$ with respect to $\mathbf{x}^0$.

We define the first corrector $\phi_1$ to cancel the $\epsilon^{-1}$ term in (84). This gives a Poisson equation for $\phi_1$ as a function of $\boldsymbol{\eta}$ with $(\mathbf{x}, \mathbf{x}^0) \in \mathbb{R}^{2n}$ a frozen parameter. The Poisson equation

$$ Qg = -\boldsymbol{\eta} \qquad (85) $$

can be solved by Fredholm alternative because the process $\boldsymbol{\eta}^1$ has mean zero (with respect to its invariant probability measure) [8, Chapter 6]. We can write a solution in the form:

$$ g(\boldsymbol{\eta}) = \int_{0}^{\infty} \mathbb{E}[\boldsymbol{\eta}_{s}^{1} | \boldsymbol{\eta}_{0}^{1} = \boldsymbol{\eta}] ds, $$

which is here linear in $\boldsymbol{\eta}$:

$$ g(\boldsymbol{\eta}) = \mathbf{A}^{-1}\boldsymbol{\eta}. \qquad (86) $$

Therefore we set

$$ \phi_1(\mathbf{x}, \mathbf{x}^0, \boldsymbol{\eta}) = \sum_{j=1}^{n} \sum_{i=1}^{d} \sigma_{ji}(\mathbf{x}) g_i(\boldsymbol{\eta}) \partial_{x_j} \phi(\mathbf{x}, \mathbf{x}^0). \qquad (87) $$

We cannot define the second corrector $\phi_2$ so as to cancel the order-one terms in (84) because that would require solving a Poisson equation with a right-hand side that is not centered. To center this term we subtract its mean relative to the invariant distribution of $\boldsymbol{\eta}^1$. This gives the Poisson equation

$$ Q\phi_2(\mathbf{x}, \mathbf{x}^0, \boldsymbol{\eta}) + \sum_{j=1}^{n} \sum_{i=1}^{d} \sigma_{ji}(\mathbf{x}) \eta_i \partial_{x_j} \phi_1(\mathbf{x}, \mathbf{x}^0, \boldsymbol{\eta}) + \sum_{j=1}^{n} \sum_{i=1}^{d} \sum_{k=1}^{d'} \Gamma_{jk}(\mathbf{x}^0) K_{ik} \partial_{\eta_i x_j^0}^2 \phi_1(\mathbf{x}, \mathbf{x}^0, \boldsymbol{\eta}) \\ - \mathbb{E} \left[ \sum_{j=1}^{n} \sum_{i=1}^{d} \sigma_{ji}(\mathbf{x}) \eta_{i,0}^1 \partial_{x_j} \phi_1(\mathbf{x}, \mathbf{x}^0, \boldsymbol{\eta}_0^1) + \sum_{j=1}^{n} \sum_{i=1}^{d} \sum_{k=1}^{d'} \Gamma_{jk}(\mathbf{x}^0) K_{ik} \partial_{\eta_i x_j^0}^2 \phi_1(\mathbf{x}, \mathbf{x}^0, \boldsymbol{\eta}_0^1) \right] = 0, $$

where the expectation $\mathbb{E}$ is taken over $\boldsymbol{\eta}_0^1$ with respect to its invariant probability measure. Since $\phi_1$ is linear in $\boldsymbol{\eta}$, the third term of the left-hand side is independent of $\boldsymbol{\eta}$ and is equal to its expectation, so the Poisson equation can be reduced to

$$ Q\phi_2(\mathbf{x}, \mathbf{x}^0, \boldsymbol{\eta}) + \sum_{j=1}^{n} \sum_{i=1}^{d} \sigma_{ji}(\mathbf{x}) \eta_i \partial_{x_j} \phi_1(\mathbf{x}, \mathbf{x}^0, \boldsymbol{\eta}) - \mathbb{E} \left[ \sum_{j=1}^{n} \sum_{i=1}^{d} \sigma_{ji}(\mathbf{x}) \eta_{i,0}^1 \partial_{x_j} \phi_1(\mathbf{x}, \mathbf{x}^0, \boldsymbol{\eta}_0^1) \right] = 0. $$
This equation has a solution $\phi_2$ that is a smooth function in $(\boldsymbol{x}, \boldsymbol{x}^0)$ and that is a quadratic form in $\boldsymbol{\eta}$. Note that $\phi_1(\boldsymbol{x}, \boldsymbol{x}^0, \boldsymbol{\eta})$ and $\phi_2(\boldsymbol{x}, \boldsymbol{x}^0, \boldsymbol{\eta})$ depend only on $\boldsymbol{\sigma}(\boldsymbol{x})$ and its first-order derivatives, and not on $\boldsymbol{b}$. By assuming that $\boldsymbol{\sigma}$ belongs to $\mathcal{C}^2$ with bounded derivatives, we get the control of the $O(\varepsilon)$ term in (84). It follows that

$$ \mathcal{L}^{\varepsilon}\phi^{\varepsilon} = \mathbb{E}\left[\sum_{j=1}^{n} \sum_{i=1}^{d} \sigma_{ji}(\boldsymbol{x})\eta_{i,0}^{1}\partial_{x_j}\phi_{1}(\boldsymbol{x}, \boldsymbol{x}^{0}, \eta_{0}^{1}) + \sum_{j=1}^{n} \sum_{i=1}^{d} \sum_{k=1}^{d'} \Gamma_{jk}(\boldsymbol{x}^{0})K_{ik}\partial_{\eta_{i}x_{j}^{0}}^{2}\phi_{1}(\boldsymbol{x}, \boldsymbol{x}^{0}, \eta_{0}^{1})\right] \\ + \left[\sum_{j=1}^{n} b_{j}(\boldsymbol{x})\partial_{x_j}\phi(\boldsymbol{x}, \boldsymbol{x}^{0}) + \sum_{j=1}^{n} \tilde{b}_{j}(\boldsymbol{x}^{0})\partial_{x_{j}^{0}}\phi(\boldsymbol{x}, \boldsymbol{x}^{0}) + \frac{1}{2}\sum_{i,j=1}^{n} \sum_{k=1}^{d'} \Gamma_{ik}(\boldsymbol{x}^{0})\Gamma_{jk}(\boldsymbol{x}^{0})\partial_{\eta_{i}x_{j}^{0}}^{2}\phi(\boldsymbol{x}, \boldsymbol{x}^{0})\right] \\ + O(\varepsilon). $$

Using (87), the expectation takes the form

$$ \mathbb{E}\left[\sum_{j=1}^{n} \sum_{i=1}^{d} \sigma_{ji}(\boldsymbol{x}) \eta_{i,0}^{1} \partial_{x_j} \phi_1(\boldsymbol{x}, \boldsymbol{x}^0, \eta_0^1) + \sum_{j=1}^{n} \sum_{i=1}^{d} \sum_{k=1}^{d'} \Gamma_{jk}(\boldsymbol{x}^0) K_{ik} \partial_{\eta_i x_j^0}^2 \phi_1(\boldsymbol{x}, \boldsymbol{x}^0, \eta_0^1)\right] \\ = \sum_{j,j'=1}^{n} \sum_{i,i'=1}^{d} \mathbb{E}[\eta_{i,0}^1 g_{i'}(\eta_0^1)] \sigma_{ji}(\boldsymbol{x}) \partial_{x_j}(\sigma_{j'i'}(\boldsymbol{x}) \partial_{x_{j'}} \phi(\boldsymbol{x}, \boldsymbol{x}^0)) \\ + \sum_{j,j'=1}^{n} \sum_{i,i'=1}^{d} \sum_{k=1}^{d'} \mathbb{E}[\partial_{\eta_i} g_{i'}(\eta_0^1)] \Gamma_{jk}(\boldsymbol{x}^0) K_{ik} \partial_{x_j^0}^2 (\sigma_{j'i'}(\boldsymbol{x}) \partial_{x_{j'}} \phi(\boldsymbol{x}, \boldsymbol{x}^0)). $$

From the explicit form (86) of $\boldsymbol{g}$ we get

$$ \mathbb{E}[\eta_{i,0}^1 g_{i'}(\eta_0^1)] = \sum_{k=1}^{d} (\mathbf{A}^{-1})_{i'k} C_{ki} = (\mathbf{A}^{-1}\mathbf{C})_{i'i} \text{ and } \mathbb{E}[\partial_{\eta_i} g_{i'}(\eta_0^1)] = (\mathbf{A}^{-1})_{i'i}. $$
Therefore

$$
\begin{align*}
& \mathbb{E} \left[ \sum_{j=1}^{n} \sum_{i=1}^{d} \sigma_{ji}(\mathbf{x}) \eta_{i,0}^{1} \partial_{x_j} \phi_1(\mathbf{x}, \mathbf{x}^0, \eta_0^1) + \sum_{j=1}^{n} \sum_{i=1}^{d} \sum_{k=1}^{d'} \Gamma_{jk}(\mathbf{x}^0) K_{ik} \partial_{\eta_i x_j^0}^2 \phi_1(\mathbf{x}, \mathbf{x}^0, \eta_0^1) \right] \\
&= \sum_{j,j'=1}^{n} \sum_{i,i'=1}^{d} (\mathbf{A}^{-1}\mathbf{C})_{i'i} \sigma_{ji}(\mathbf{x}) \partial_{x_j}(\sigma_{j'i'}(\mathbf{x}) \partial_{x_{j'}} \phi(\mathbf{x}, \mathbf{x}^0)) \\
&\quad + \sum_{j,j'=1}^{n} \sum_{i,i'=1}^{d} \sum_{k=1}^{d'} \Gamma_{jk}(\mathbf{x}^0) K_{ik} \sigma_{j'i'}(\mathbf{x}) (\mathbf{A}^{-1})_{i'i} \partial_{x_j^0 x_{j'}} \phi(\mathbf{x}, \mathbf{x}^0) \\
&= \sum_{j,j'=1}^{n} (\boldsymbol{\sigma}(\mathbf{x}) \mathbf{A}^{-1} \mathbf{C} \boldsymbol{\sigma}(\mathbf{x})^T)_{j'j} \partial_{x_{j'},x_j} \phi(\mathbf{x}, \mathbf{x}^0) \\
&\quad + \sum_{j'=1}^{n} \left( \sum_{j=1}^{n} (\partial_{x_j} \boldsymbol{\sigma}(\mathbf{x})) \mathbf{A}^{-1} \mathbf{C} \boldsymbol{\sigma}(\mathbf{x})^T \right)_{j'j} \partial_{x_{j'}^0 x_{j'}} \phi(\mathbf{x}, \mathbf{x}^0) \\
&\quad + \sum_{j,j'=1}^{n} (\boldsymbol{\sigma}(\mathbf{x}) \mathbf{A}^{-1} \mathbf{K} \boldsymbol{\Gamma}(\mathbf{x}^0)^T)_{j'j} \partial_{x_j^0 x_{j'}} \phi(\mathbf{x}, \mathbf{x}^0) \\
&= \frac{1}{2} \sum_{j,j'=1}^{n} (\boldsymbol{\sigma}(\mathbf{x}) (\mathbf{A}^{-1}\mathbf{C} + \mathbf{C}\mathbf{A}^{T-1}) \boldsymbol{\sigma}(\mathbf{x})^T)_{jj'} \partial_{x_{j'},x_j} \phi(\mathbf{x}, \mathbf{x}^0) \\
&\quad + \sum_{j'=1}^{n} (\tilde{b}_{j'}(\mathbf{x}) - b_{j'}(\mathbf{x})) \partial_{x_{j'}} \phi(\mathbf{x}, \mathbf{x}^0) + \sum_{j,j'=1}^{n} (\boldsymbol{\Gamma}(\mathbf{x}) \boldsymbol{\Gamma}(\mathbf{x}^0)^T)_{j'j'} \partial_{x_j^0 x_{j'}} \phi(\mathbf{x}, \mathbf{x}^0).
\end{align*}
$$

We have from $d\eta^1 = \mathbf{K}d\mathbf{W}_t - \mathbf{A}\eta^1 dt$ and Itô's formula:

$$
d(\eta_i^1 \eta_j^1) = \sum_{k=1}^{d'} (K_{ik} \eta_j^1 + K_{jk} \eta_i^1) dW_k t + \sum_{k=1}^{d'} K_{ik} K_{jk} dt - \sum_{k=1}^{d} (A_{ik} \eta_k^1 \eta_j^1 + A_{jk} \eta_k^1 \eta_i^1) dt.
$$

Taking the expectation (under the invariant probability measure) gives the identity:

$$
\mathbf{K}\mathbf{K}^T - \mathbf{C}\mathbf{A}^T - \mathbf{A}\mathbf{C} = \mathbf{0}.
$$

By left-multiplying by $\mathbf{A}^{-1}$ and by right-multiplying by $\mathbf{A}^{T-1}$ we find

$$
\mathbf{A}^{-1}\mathbf{K}\mathbf{K}^T\mathbf{A}^{T-1} - \mathbf{A}^{-1}\mathbf{C} - \mathbf{C}\mathbf{A}^{T-1} = \mathbf{0},
$$

which gives

$$
\sigma(\mathbf{x})(\mathbf{A}^{-1}\mathbf{C} + \mathbf{C}\mathbf{A}^{T-1})\sigma(\mathbf{x})^T = \mathbf{\Gamma}(\mathbf{x})\mathbf{\Gamma}(\mathbf{x})^T,
$$

and we obtain the desired result:

$$
\mathcal{L}^{\epsilon}\phi^{\epsilon} = \mathcal{L}\phi + O(\epsilon).
$$

<figure><img src="image_11.png" /></figure>

We can then prove Proposition 3.5 as follows.

*Proof of Proposition 3.5.* By the perturbed test function method [8, Section 6.3], Lemma A.1 establishes that the continuous process $(\mathbf{X}^{\epsilon}, \mathbf{X}^0)$ converges in distribution to the Markov process
with infinitesimal generator $\mathcal{L}$ defined by (82). The infinitesimal generator $\mathcal{L}$ can be associated to a diffusion process $(\tilde{\mathbf{X}}, \tilde{\mathbf{X}}^0)$ that is solution of the coupled SDEs:

$$d\tilde{\mathbf{X}} = \tilde{b}(\tilde{\mathbf{X}})dt + \Gamma(\tilde{\mathbf{X}})d\tilde{W}_t,$$

$$d\tilde{\mathbf{X}}^0 = \tilde{b}(\tilde{\mathbf{X}}^0)dt + \Gamma(\tilde{\mathbf{X}}^0)d\tilde{W}_t,$$

where $\tilde{W}$ is a $d'$-dimensional Brownian motion. This shows that, if $\tilde{\mathbf{X}}_0 = \tilde{\mathbf{X}}_0^0$ almost surely, then $\tilde{\mathbf{X}}_t - \tilde{\mathbf{X}}_t^0 = 0$ for all $t$ almost surely. Therefore, if $\mathbf{X}_0^\epsilon = \mathbf{X}_0^0$, then the continuous process $(\mathbf{X}^\epsilon - \mathbf{X}^0)$ converges in distribution to 0, which implies convergence in probability. $\square$

## B Proof of Lemma 3.8

We define

$$\phi(\mathbf{x}, \mathbf{x}^0) = g(\mathbf{x}^0)(f(\mathbf{x}) - f(\mathbf{x}^0)) \quad \text{or} \quad \phi(\mathbf{x}, \mathbf{x}^0) = (f(\mathbf{x}) - f(\mathbf{x}^0))^2. \qquad (88)$$

Lemma A.1 applied to $\phi$ gives an estimate for (35) of order $\epsilon$, but the particular form of $\phi$ makes it possible to get $\epsilon^2$, as we show in the following. We prove Lemma 3.8 in four steps.

Step 1. There exist smooth functions $\phi_{1i}$, $\phi_{20}$, $\phi_{2ij}$, $\Lambda_{1i}$, $\Lambda_{1ijk}$ with bounded derivatives such that

$$\phi^{\epsilon}(\mathbf{x}, \mathbf{x}^0, \eta) = \phi(\mathbf{x}, \mathbf{x}^0) + \epsilon\phi_1(\mathbf{x}, \mathbf{x}^0, \eta) + \epsilon^2\phi_2(\mathbf{x}, \mathbf{x}^0, \eta), \qquad (89)$$

$$\phi_1(\mathbf{x}, \mathbf{x}^0, \eta) = \sum_{i=1}^{d} \phi_{1i}(\mathbf{x}, \mathbf{x}^0)\eta_i, \qquad (90)$$

$$\phi_2(\mathbf{x}, \mathbf{x}^0, \eta) = \phi_{20}(\mathbf{x}, \mathbf{x}^0) + \sum_{i,j=1}^{d} \phi_{2ij}(\mathbf{x}, \mathbf{x}^0)\eta_i\eta_j, \qquad (91)$$

$$\mathcal{L}^{\epsilon}\phi^{\epsilon}(\mathbf{x}, \mathbf{x}^0, \eta) = \epsilon\Lambda_1(\mathbf{x}, \mathbf{x}^0, \eta) + O(\epsilon^2), \qquad (92)$$

$$\Lambda_1(\mathbf{x}, \mathbf{x}^0, \eta) = \sum_{i=1}^{d} \Lambda_{1i}(\mathbf{x}, \mathbf{x}^0)\eta_i + \sum_{i,j,k=1}^{d} \Lambda_{1ijk}(\mathbf{x}, \mathbf{x}^0)\eta_i\eta_j\eta_k. \qquad (93)$$

Proof. We apply the perturbed test function method as described in the proof of Lemma A.1 and we get the result by keeping track of the $\eta$-dependence of the perturbed functions $\phi_1$ and $\phi_2$. $\square$

Step 2. For $s \le t$, the conditional distribution of $\eta_t^\epsilon$ given $\mathcal{F}_s = \sigma(\mathbf{W}_u, u \le s)$ is

$$\mathcal{N}\left(\exp\left(-\frac{\mathbf{A}(t-s)}{\epsilon^2}\right)\eta_s^\epsilon, \int_0^{(t-s)/\epsilon^2} e^{-\mathbf{A}u}\mathbf{K}\mathbf{K}^T e^{-\mathbf{A}^T u}du\right). \qquad (94)$$

There exists $\lambda, C > 0$ such that

$$|\mathbb{E}[\eta_{i,t}^{\epsilon}|\mathcal{F}_s]| \le e^{-\lambda(t-s)/\epsilon^2} \|\eta_s^{\epsilon}\|, \qquad (95)$$

$$|\mathbb{E}[\eta_{i,t}^{\epsilon}\eta_{j,t}^{\epsilon}\eta_{k,t}^{\epsilon}|\mathcal{F}_s]| \le C e^{-\lambda(t-s)/\epsilon^2} \|\eta_s^{\epsilon}\|(1 + \|\eta_s^{\epsilon}\|^2). \qquad (96)$$

Proof. We can integrate (2) from $s$ to $t$:

$$\eta_t^\epsilon = \exp\left(-\frac{\mathbf{A}(t-s)}{\epsilon^2}\right)\eta_s^\epsilon + \int_s^t \exp\left(-\frac{\mathbf{A}(t-u)}{\epsilon^2}\right)\frac{\mathbf{K}}{\epsilon}d\mathbf{W}_u,$$
which gives (94) and

$$ \mathbb{E}[\eta_t^\varepsilon | \mathcal{F}_s] = \exp\left(-\frac{\mathbf{A}(t-s)}{\varepsilon^2}\right) \eta_s^\varepsilon. $$

Eq. (95) is a straightforward consequence. Eq. (96) follows from (94) and Isserlis theorem for multivariate normal random vectors. □

Step 3. If $(\mathbf{x}, \mathbf{x}^0) \mapsto \psi(\mathbf{x}, \mathbf{x}^0)$ is a smooth function with bounded derivatives, then there exists $C > 0$ such that, for all $i, j, k = 1, \dots, d$ and $t \in [\varepsilon, T]$:

$$ |\mathbb{E}[\eta_{i,t}^{\varepsilon} \psi(\mathbf{X}_{t}^{\varepsilon}, \mathbf{X}_{t}^{0})]| \le C\varepsilon, \quad (97) $$

$$ |\mathbb{E}[\eta_{i,t}^{\varepsilon} \eta_{j,t}^{\varepsilon} \eta_{k,t}^{\varepsilon} \psi(\mathbf{X}_{t}^{\varepsilon}, \mathbf{X}_{t}^{0})]| \le C\varepsilon. \quad (98) $$

Proof. We have for any $t \in [\varepsilon, T]$ and $\delta < \varepsilon$:

$$ \begin{aligned} \psi(\mathbf{X}_t^\varepsilon, \mathbf{X}_t^0) ={}& \psi(\mathbf{X}_{t-\delta}^\varepsilon, \mathbf{X}_{t-\delta}^0) + \int_{t-\delta}^{\delta} \frac{1}{\varepsilon} \psi_1(\mathbf{X}_s^\varepsilon, \mathbf{X}_s^0, \eta_s^\varepsilon) ds \\ & + \int_{t-\delta}^{\delta} \psi_2(\mathbf{X}_s^\varepsilon, \mathbf{X}_s^0) ds + \sum_{j=1}^{d} \int_{t-\delta}^{\delta} \psi_{3j}(\mathbf{X}_s^\varepsilon, \mathbf{X}_s^0) dW_{js}, \end{aligned} $$

with

$$ \begin{aligned} \psi_1(\mathbf{x}, \mathbf{x}^0, \boldsymbol{\eta}) &= \boldsymbol{\sigma}(\mathbf{x}) \boldsymbol{\eta} \cdot \nabla_{\mathbf{x}} \psi(\mathbf{x}, \mathbf{x}^0), \\ \psi_2(\mathbf{x}, \mathbf{x}^0) &= [\mathbf{b}(\mathbf{x}) \cdot \nabla_{\mathbf{x}} \psi + \tilde{\mathbf{b}}(\mathbf{x}^0) \cdot \nabla_{\mathbf{x}^0} \psi + \nabla_{\mathbf{x}^0} \psi^T \mathbf{\Gamma}(\mathbf{x}^0)^T \mathbf{\Gamma}(\mathbf{x}^0) \nabla_{\mathbf{x}^0} \psi](\mathbf{x}, \mathbf{x}^0), \\ \psi_{3j}(\mathbf{x}, \mathbf{x}^0) &= \sum_{i=1}^{n} \partial_{\mathbf{x}_i^0} \psi(\mathbf{x}, \mathbf{x}^0) \Gamma_{ij}(\mathbf{x}^0). \end{aligned} $$

We have, by (95),

$$ |\mathbb{E}[\eta_{i,t}^{\varepsilon} \psi(\mathbf{X}_{t-\delta}^{\varepsilon}, \mathbf{X}_{t-\delta}^{0})]| = |\mathbb{E}[\mathbb{E}[\eta_{i,t}^{\varepsilon} | \mathcal{F}_{t-\delta}] \psi(\mathbf{X}_{t-\delta}^{\varepsilon}, \mathbf{X}_{t-\delta}^{0})]| \le C \exp(-\lambda \delta / \varepsilon^2). $$

Similarly, for any $s \in [t - \delta, t]$

$$ \begin{aligned} |\mathbb{E}[\eta_{i,t}^{\varepsilon} \psi_1(\mathbf{X}_s^{\varepsilon}, \mathbf{X}_s^0, \eta_s^{\varepsilon})]| &= |\mathbb{E}[\mathbb{E}[\eta_{i,t}^{\varepsilon} | \mathcal{F}_s] \psi_1(\mathbf{X}_s^{\varepsilon}, \mathbf{X}_s^0, \eta_s^{\varepsilon})]| \le C \exp(-\lambda(t-s)/\varepsilon^2), \\ |\mathbb{E}[\eta_{i,t}^{\varepsilon} \psi_2(\mathbf{X}_s^{\varepsilon}, \mathbf{X}_s^0)]| &\le C \exp(-\lambda(t-s)/\varepsilon^2), \end{aligned} $$

and for $j = 1, \dots, d'$ and for any positive $q$,

$$ \begin{aligned} & |\mathbb{E}[\eta_{i,t}^{\varepsilon} \int_{t-\delta}^{t} \psi_{3j}(\mathbf{X}_s^{\varepsilon}, \mathbf{X}_s^0) dW_{js}]| \\ & \le \sum_{k=0}^{q-1} |\mathbb{E}[\mathbb{E}[\eta_{i,t}^{\varepsilon} | \mathcal{F}_{t-k\delta/q}] \int_{t-(k+1)\delta/q}^{t-k\delta/q} \psi_{3j}(\mathbf{X}_s^{\varepsilon}, \mathbf{X}_s^0) dW_{js}]| \\ & \le \sum_{k=0}^{q-1} |\mathbb{E}[\mathbb{E}[\eta_{i,t}^{\varepsilon} | \mathcal{F}_{t-k\delta/q}]^2]^{1/2} \left[ \int_{t-(k+1)\delta/q}^{t-k\delta/q} \mathbb{E}[\psi_{3j}(\mathbf{X}_s^{\varepsilon}, \mathbf{X}_s^0)^2] ds \right]^{1/2} \\ & \le C \sum_{k=0}^{q-1} \exp(-\lambda k\delta/(q\varepsilon^2)) \sqrt{\delta/q}. \end{aligned} $$
Consequently

$$
\begin{align*}
& |\mathbb{E}[\eta_{i,t}^{\varepsilon} \psi(\mathbf{X}_t^{\varepsilon}, \mathbf{X}_t^0)]| \\
& \le C \exp\left(-\frac{\lambda \delta}{\varepsilon^2}\right) + \frac{C}{\varepsilon} \int_{t-\delta}^t \exp\left(-\frac{\lambda(t-s)}{\varepsilon^2}\right) ds + C \frac{\sqrt{\delta}}{\sqrt{q}} \sum_{k=0}^{q-1} \exp\left(-\frac{\lambda k \delta}{q \varepsilon^2}\right) \\
& \le C \exp\left(-\frac{\lambda \delta}{\varepsilon^2}\right) + \frac{C \varepsilon}{\lambda} + C \frac{\sqrt{\delta}}{\sqrt{q}} \sum_{k=0}^{q-1} \exp\left(-\frac{\lambda k \delta}{q \varepsilon^2}\right).
\end{align*}
$$

By taking $\delta = \varepsilon^2 |\ln \varepsilon|/\lambda$ and $q = [|\ln \varepsilon|]$ we finally get

$$
|\mathbb{E}[\eta_{i,t}^{\epsilon} \psi(\mathbf{X}_t^{\epsilon}, \mathbf{X}_t^0)]| \le C' \epsilon + C' \epsilon \sum_{k=0}^{\infty} e^k \le C'' \epsilon,
$$

which gives the first desired result. The calculations with the third-order product of coefficients $\eta$
are similar and use (96). $\square$

Step 4. *Proof of Lemma* 3.8.

For any $t \in [\varepsilon, T]$, we have

$$
\begin{align*}
\mathbb{E}[\phi(\mathbf{X}_t^\varepsilon, \mathbf{X}_t^0)] &= \mathbb{E}[\phi^\varepsilon(\mathbf{X}_t^\varepsilon, \mathbf{X}_t^0, \eta_t^\varepsilon)] - \varepsilon\mathbb{E}[\phi_1(\mathbf{X}_t^\varepsilon, \mathbf{X}_t^0, \eta_t^\varepsilon)] + O(\varepsilon^2) \\
&= \mathbb{E}[\phi^\varepsilon(\mathbf{X}_t^\varepsilon, \mathbf{X}_t^0, \eta_t^\varepsilon)] + O(\varepsilon^2),
\end{align*}
$$

because (90) and (97) give $\mathbb{E}[\phi_1(\mathbf{X}_t^\varepsilon, \mathbf{X}_t^0, \eta_t^\varepsilon)] = O(\varepsilon)$.

We have

$$
\mathbb{E}[\phi^{\epsilon}(\mathbf{x}_0, \mathbf{x}_0, \eta_0^{\epsilon})] = \epsilon \mathbb{E}[\phi_1(\mathbf{x}_0, \mathbf{x}_0, \eta_0^{\epsilon})] + O(\epsilon^2) = O(\epsilon^2),
$$

because $\mathbb{E}[\phi_{1j}(\mathbf{x}_0, \mathbf{x}_0)\eta_{j,0}^{\epsilon}] = \phi_{1j}(\mathbf{x}_0, \mathbf{x}_0)\mathbb{E}[\eta_{j,0}^{\epsilon}] = 0$ for all $j = 1, \dots, d$.

Therefore

$$
\begin{align*}
\mathbb{E}[\phi(\mathbf{X}_t^\varepsilon, \mathbf{X}_t^0)] &= \int_0^t \mathbb{E}[\mathcal{L}^\varepsilon \phi^\varepsilon(\mathbf{X}_s^\varepsilon, \mathbf{X}_s^0, \boldsymbol{\eta}_s^\varepsilon)] ds + O(\varepsilon^2) \\
&= \varepsilon \int_0^t \mathbb{E}[\Lambda_1(\mathbf{X}_s^\varepsilon, \mathbf{X}_s^0, \boldsymbol{\eta}_s^\varepsilon)] ds + O(\varepsilon^2) \\
&= O(\varepsilon^2),
\end{align*}
$$

because (93), (97) and (98) give $\mathbb{E}[\Lambda_1(\mathbf{X}_s^\varepsilon, \mathbf{X}_s^0, \boldsymbol{\eta}_s^\varepsilon)] = O(\varepsilon)$ for any $s \in [\varepsilon, t]$ and $\mathbb{E}[\Lambda_1(\mathbf{X}_s^\varepsilon, \mathbf{X}_s^0, \boldsymbol{\eta}_s^\varepsilon)] = O(1)$ for any $s \in [0, \varepsilon]$. This completes the proof of Lemma 3.8. $\square$

C Proofs of existence and uniqueness of (30) and (31)

**Proposition C.1.** Fix $T > 0$, $n \in \mathbb{N}^*$. Suppose that $\mathbf{f} \in C([0, T]; \mathbb{R}^n)$, $\mathbf{b}$ is Lipschitz, and $\varphi$ is a l.s.c. convex function satisfying (51). Then there exists a unique solution $\mathbf{x} \in C([0, T]; \mathbb{R}^n)$ to the following differential inclusion

$$
\mathbf{x}(0) = \mathbf{x}_0 \in \mathbb{R}^n, \quad \dot{\mathbf{x}}(t) + \partial\varphi(\mathbf{x}(t)) \ni \mathbf{b}(\mathbf{x}(t)) + \mathbf{f}(t), \quad t > 0. \quad (99)
$$
Proof. Let $\varphi_p$ be the Moreau-Yosida regularization of $\varphi$. For each $p \ge 1$, we consider the penalized problem

$$ \boldsymbol{x}^p(0) = \boldsymbol{x}_0 \in \mathbb{R}^n, \quad \dot{\boldsymbol{x}}^p(t) + \nabla \varphi_p(\boldsymbol{x}^p(t)) = \boldsymbol{b}(\boldsymbol{x}^p(t)) + \boldsymbol{f}(t), \ t > 0. $$

This is a standard ODE with Lipschitz coefficients, so $\boldsymbol{x}^p \in C([0, T]; \mathbb{R}^n)$ is well-defined. Now, we show that $\boldsymbol{x}^p$ is a Cauchy sequence in $C([0, T]; \mathbb{R}^n)$. Fix $p, q \in \mathbb{N}^*$ and $t \in [0, T]$. We have the following expansion

$$ \frac{1}{2} \|\boldsymbol{x}^p(t) - \boldsymbol{x}^q(t)\|^2 = \int_0^t (\boldsymbol{x}^p(s) - \boldsymbol{x}^q(s)) \cdot (\boldsymbol{b}(\boldsymbol{x}^p(s)) - \boldsymbol{b}(\boldsymbol{x}^q(s))) \mathrm{d}s \\ - \int_0^t (\boldsymbol{x}^p(s) - \boldsymbol{x}^q(s)) \cdot (\nabla \varphi_p(\boldsymbol{x}^p(s)) - \nabla \varphi_q(\boldsymbol{x}^q(s))) \mathrm{d}s $$

which, using the properties of $\boldsymbol{b}$, $\varphi_p$, $\varphi_q$, leads to the following inequality

$$ \frac{1}{2} \|\boldsymbol{x}^p(t) - \boldsymbol{x}^q(t)\|^2 \le C \int_0^t \|\boldsymbol{x}^p(s) - \boldsymbol{x}^q(s)\|^2 \mathrm{d}s + \left(\frac{1}{p} + \frac{1}{q}\right) \int_0^t \nabla \varphi_p(\boldsymbol{x}^p(s)) \cdot \nabla \varphi_q(\boldsymbol{x}^q(s)) \mathrm{d}s. \quad (100) $$

Under the assumption $\sup_{p \ge 1} \|\nabla \varphi_p(\boldsymbol{x})\| < \infty$, we deduce from the inequality above that

$$ \frac{1}{2} \|\boldsymbol{x}^p(t) - \boldsymbol{x}^q(t)\|^2 \le C \int_0^t \|\boldsymbol{x}^p(s) - \boldsymbol{x}^q(s)\|^2 \mathrm{d}s + \left(\frac{1}{p} + \frac{1}{q}\right) Ct. $$

Thus, we can apply Gronwall's inequality to obtain

$$ \sup_{0 \le t \le T} \|\boldsymbol{x}^p(t) - \boldsymbol{x}^q(t)\|^2 \le \left(\frac{1}{p} + \frac{1}{q}\right) C_T. \quad (101) $$

Therefore $\boldsymbol{x}^p$ is a Cauchy sequence and there exists a function $\boldsymbol{x} \in C([0, T]; \mathbb{R}^n)$ such that $\boldsymbol{x}^p \to \boldsymbol{x}$, as $p \to \infty$ in $C([0, T]; \mathbb{R}^n)$. Next we verify that $\boldsymbol{x}$ satisfies the differential inclusion. Define $\forall t \in [0, T]$, $\boldsymbol{\Delta}^p(t) \triangleq \int_0^t \nabla \varphi_p(\boldsymbol{x}^p(s)) \mathrm{d}s$ and denote $\boldsymbol{\Delta}(t) \triangleq \lim_{p \to \infty} \boldsymbol{\Delta}^p(t)$. We then have

$$ \forall t \in [0, T], \quad \boldsymbol{x}(t) + \boldsymbol{\Delta}(t) = \boldsymbol{x}_0 + \int_0^t \boldsymbol{b}(\boldsymbol{x}(s)) \mathrm{d}s + \int_0^t \boldsymbol{f}(s) \mathrm{d}s. $$

Moreover, since $\sup_p \int_0^T \|\nabla \varphi_p(\boldsymbol{x}^p(s))\|^2 \mathrm{d}s < \infty$, there exists a function $\boldsymbol{\delta} \in L^2(0, T)$ such that

$$ \forall \boldsymbol{h} \in L^2(0, T), \quad \lim_{p \to \infty} \int_0^T \nabla \varphi_p(\boldsymbol{x}^p(s)) \cdot \boldsymbol{h}(s) \mathrm{d}s = \int_0^T \boldsymbol{\delta}(s) \cdot \boldsymbol{h}(s) \mathrm{d}s. $$

As a consequence, we must have $\boldsymbol{\Delta}(t) = \int_0^t \boldsymbol{\delta}(s) \mathrm{d}s$. Now, to finally check the differential inclusion, we want to show that $\forall \boldsymbol{v} \in C([0, T]; \mathbb{R}^n)$,

$$ \forall 0 \le t < t + h \le T, \quad \int_t^{t+h} \boldsymbol{\delta}(s) \cdot (\boldsymbol{v}(s) - \boldsymbol{x}(s)) + \varphi(\boldsymbol{x}(s)) \mathrm{d}s \le \int_t^{t+h} \varphi(\boldsymbol{v}(s)) \mathrm{d}s. $$

We exploit property 2 as listed in section 5.1

$$ \forall t \in [0, T], \quad \nabla \varphi_p(\boldsymbol{x}^p(t)) \in \partial \varphi(J_p \boldsymbol{x}^p(t)) $$
which implies $\forall \mathbf{v} \in C([0, T]; \mathbb{R}^n)$,

$$ \forall t \in [0, T], \quad \nabla \varphi_p(\mathbf{x}^p(t)) \cdot (\mathbf{v}(t) - J_p \mathbf{x}^p(t)) + \varphi(J_p \mathbf{x}^p(t)) \le \varphi(\mathbf{v}(t)). $$

We integrate the inequality above on $[t, t + h]$ and obtain $\forall \mathbf{v} \in C([0, T]; \mathbb{R}^n)$,

$$ \forall 0 \le t < t + h \le T, \quad \int_{t}^{t+h} \nabla \varphi_p(\mathbf{x}^p(s)) \cdot (\mathbf{v}(s) - J_p \mathbf{x}^p(s)) + \varphi(J_p \mathbf{x}^p(s)) \mathrm{d}s \le \int_{t}^{t+h} \varphi(\mathbf{v}(s)) \mathrm{d}s. $$

On the one hand, by using Fatou's lemma ($\varphi \ge 0$), the convergence $J_p \mathbf{x}^p \to \mathbf{x}$ in $C([0, T]; \mathbb{R}^n)$ as $p \to +\infty$, and the fact that $\varphi$ is l.s.c, we get

$$ \forall 0 \le t < t + h \le T, \quad \liminf_{p \to \infty} \int_{t}^{t+h} \varphi(J_p \mathbf{x}^p(s)) \mathrm{d}s \ge \int_{t}^{t+h} \liminf_{p \to \infty} \varphi(J_p \mathbf{x}^p(s)) \mathrm{d}s \ge \int_{t}^{t+h} \varphi(\mathbf{x}(s)) \mathrm{d}s. $$

On the other hand, we have

$$ \int_{t}^{t+h} \nabla \varphi_p(\mathbf{x}^p(s)) \cdot (\mathbf{v}(s) - J_p \mathbf{x}^p(s)) \mathrm{d}s = \int_{t}^{t+h} \nabla \varphi_p(\mathbf{x}^p(s)) \cdot (\mathbf{v}(s) - \mathbf{x}(s)) \mathrm{d}s + \int_{t}^{t+h} \nabla \varphi_p(\mathbf{x}^p(s)) \cdot (\mathbf{x}(s) - J_p \mathbf{x}^p(s)) \mathrm{d}s. $$

By using $\nabla \varphi_p(\mathbf{x}^p(s)) \to \delta$ in $L^2(0, T)$ weak, the first term of the right-hand side goes to zero. By using $\|\nabla \varphi_p(\mathbf{x}^p(s))\|$ is uniformly bounded in $L^2(0, T)$ and $\|J_p \mathbf{x}^p - \mathbf{x}\| \to 0$ in $L^2(0, T)$, the second term of the right-hand side goes to zero. Therefore, we have

$$ \liminf_{p \to \infty} \int_{t}^{t+h} \nabla \varphi_p(\mathbf{x}^p(s)) \cdot (\mathbf{v}(s) - J_p \mathbf{x}^p(s)) + \varphi(J_p \mathbf{x}^p(s)) \mathrm{d}s \ge \int_{t}^{t+h} \delta(s) \cdot (\mathbf{v}(s) - \mathbf{x}(s)) + \varphi(\mathbf{x}(s)) \mathrm{d}s, $$

which proves that $\mathbf{x}$ is solution of (99).

We now show uniqueness. Assume that $\mathbf{x}$ and $\mathbf{y}$ satisfy (99) with $\mathbf{x}(0) = \mathbf{x}_0$ and $\mathbf{y}(0) = \mathbf{y}_0$. Then

$$ \mathbf{x}(t) + \Delta_x(t) = \mathbf{x}_0 + \int_0^t \mathbf{b}(\mathbf{x}(s)) \mathrm{d}s + \int_0^t \mathbf{f}(s) \mathrm{d}s $$

and

$$ \mathbf{y}(t) + \Delta_y(t) = \mathbf{y}_0 + \int_0^t \mathbf{b}(\mathbf{y}(s)) \mathrm{d}s + \int_0^t \mathbf{f}(s) \mathrm{d}s. $$

With $\int_0^t \delta_x(s) \mathrm{d}s = \Delta_x(t)$ and $\int_0^t \delta_y(s) \mathrm{d}s = \Delta_y(t)$, we have the following inequalities

$$ \forall 0 \le t < t + h \le T, \quad \int_{t}^{t+h} \delta_x(s) \cdot (\mathbf{y}(s) - \mathbf{x}(s)) + \varphi(\mathbf{x}(s)) \mathrm{d}s \le \int_{t}^{t+h} \varphi(\mathbf{y}(s)) \mathrm{d}s $$

and

$$ \forall 0 \le t < t + h \le T, \quad \int_{t}^{t+h} \delta_y(s) \cdot (\mathbf{x}(s) - \mathbf{y}(s)) + \varphi(\mathbf{y}(s)) \mathrm{d}s \le \int_{t}^{t+h} \varphi(\mathbf{x}(s)) \mathrm{d}s, $$

which give

$$ \forall 0 \le t < t + h \le T, \quad \int_{t}^{t+h} (\delta_x(s) - \delta_y(s)) \cdot (\mathbf{y}(s) - \mathbf{x}(s)) \mathrm{d}s \le 0. $$
That means

$$ \forall 0 \le t < t + h \le T, \quad \int_t^{t+h} (\mathbf{b}(\mathbf{x}(s)) - \mathbf{b}(\mathbf{y}(s)) - (\dot{\mathbf{x}}(s) - \dot{\mathbf{y}}(s))) \cdot (\mathbf{y}(s) - \mathbf{x}(s)) \mathrm{d}s \le 0, $$

which turns into

$$ \forall 0 \le t < t + h \le T, \quad \|\mathbf{x}(t+h) - \mathbf{y}(t+h)\|^2 \le \|\mathbf{x}(t) - \mathbf{y}(t)\|^2 + C \int_t^{t+h} \|\mathbf{x}(s) - \mathbf{y}(s)\|^2 \mathrm{d}s. $$

Gronwall inequality gives the desired result and it is clear that if $\mathbf{x}_0 = \mathbf{y}_0$ then $\forall t \in [0, T], \mathbf{x}(t) = \mathbf{y}(t)$. □

**Proposition C.2.** Fix $T > 0, n, m \in \mathbb{N}^*$. Suppose that $\mathbf{f} \in C([0, T]; \mathbb{R}^n)$, $\mathbf{b}^x : \mathbb{R}^n \to \mathbb{R}$ and $\mathbf{b}^z : \mathbb{R}^m \to \mathbb{R}$ are Lipschitz, and $\varphi : \mathbb{R}^n \to \mathbb{R}$, $\psi : \mathbb{R}^m \to \mathbb{R}$ are l.s.c. convex functions, with $\varphi$ satisfying (51) and $\psi$ satisfying (52). Then there exists a unique solution $(\mathbf{x}, \mathbf{z}) \in C([0, T]; \mathbb{R}^n \times \mathbb{R}^m)$ to the following differential inclusion $(\mathbf{x}(0), \mathbf{z}(0)) = (\mathbf{x}_0, \mathbf{z}_0) \in \mathbb{R}^n \times \mathbb{R}^m$,

$$ \dot{\mathbf{x}}(t) + \partial\varphi(\mathbf{x}(t)) \ni \mathbf{b}^x(\mathbf{x}(t), \mathbf{z}(t)) + \mathbf{f}(t), \quad \dot{\mathbf{z}}(t) + \partial\psi(\mathbf{z}(t)) \ni \mathbf{b}^z(\mathbf{x}(t), \mathbf{z}(t)), \quad t > 0. \quad (102) $$

*Proof.* Let $\varphi_p$ and $\psi_p$ be the Moreau-Yosida regularization of $\varphi$ and $\psi$. We consider the penalized problems $(\mathbf{x}^p(0), \mathbf{z}^p(0)) = (\mathbf{x}_0, \mathbf{z}_0) \in \mathbb{R}^n \times \mathbb{R}^m$,

$$ \left. \begin{array}{l} \dot{\mathbf{x}}^p(t) + \nabla\varphi_p(\mathbf{x}^p(t)) = \mathbf{b}^x(\mathbf{x}^p(t), \mathbf{z}^p(t)) + \mathbf{f}(t), \\ \dot{\mathbf{z}}^p(t) + \nabla\psi_p(\mathbf{z}^p(t)) = \mathbf{b}^z(\mathbf{x}^p(t), \mathbf{z}^p(t)), \end{array} \right\} t > 0. \quad (103) $$

It can be shown than $(\mathbf{x}^p, \mathbf{z}^p)$ is a Cauchy sequence in $C([0, T]; \mathbb{R}^n \times \mathbb{R}^m)$. The proof follows similar steps as in the proof of Theorem C.1, except that we need a bound of the form

$$ \int_0^t \nabla\psi_p(\mathbf{z}^p(s)) \cdot \nabla\psi_q(\mathbf{z}^q(s)) \le Ct, \quad (104) $$

where $C$ does not depend on $p, q$, without using $\sup_{p \ge 1} \sup_{\mathbf{z} \in \mathbb{R}^m} \|\nabla\psi_p(\mathbf{z})\| < \infty$, which we do not assume. We proceed first with the following expansion

$$ \psi_p(\mathbf{z}^p(t)) = \psi_p(\mathbf{z}_0) + \int_0^t \nabla\psi_p(\mathbf{z}^p(s)) \cdot (\mathbf{b}^z(\mathbf{x}^p(s), \mathbf{z}^p(s)) - \nabla\psi_p(\mathbf{z}^p(s))) \mathrm{d}s, $$

which implies

$$ \psi_p(\mathbf{z}^p(t)) + \int_0^t \|\nabla\psi_p(\mathbf{z}^p(s))\|^2 \mathrm{d}s = \psi_p(\mathbf{z}_0) + \int_0^t \nabla\psi_p(\mathbf{z}^p(s)) \cdot \mathbf{b}^z(\mathbf{x}^p(s), \mathbf{z}^p(s)) \mathrm{d}s. $$

We get

$$ \psi_p(\mathbf{z}^p(t)) + \frac{1}{2} \int_0^t \|\nabla\psi_p(\mathbf{z}^p(s))\|^2 \mathrm{d}s \le \psi_p(\mathbf{z}_0) + \frac{1}{2} \int_0^t \|\mathbf{b}^z(\mathbf{x}^p(s), \mathbf{z}^p(s))\|^2 \mathrm{d}s $$

which implies

$$ \psi_p(\mathbf{z}^p(t)) + \frac{1}{2} \int_0^t \|\nabla\psi_p(\mathbf{z}^p(s))\|^2 \mathrm{d}s \le C \left( 1 + \int_0^t \{\|\mathbf{x}^p(s)\|^2 + \|\mathbf{z}^p(s)\|^2\} \mathrm{d}s \right). \quad (105) $$
Besides,

$$ \begin{aligned} \|\boldsymbol{x}^p(t)\|^2 + \|\boldsymbol{z}^p(t)\|^2 = \|\boldsymbol{x}^p(0)\|^2 + \|\boldsymbol{z}^p(0)\|^2 &+ 2 \int_0^t \boldsymbol{x}^p(s) \cdot (\boldsymbol{b}^x(\boldsymbol{x}^p(s), \boldsymbol{z}^p(s)) - \nabla \varphi_p(\boldsymbol{x}^p(s)) + \boldsymbol{f}(s)) \mathrm{d}s \\ &+ 2 \int_0^t \boldsymbol{z}^p(s) \cdot (\boldsymbol{b}^z(\boldsymbol{x}^p(s), \boldsymbol{z}^p(s)) - \nabla \psi_p(\boldsymbol{z}^p(s))) \mathrm{d}s. \end{aligned} $$

Using $\int_0^T \|\boldsymbol{f}(s)\|^2 \mathrm{d}s < \infty$, $\forall \boldsymbol{x} \in \mathbb{R}^n$, $\boldsymbol{x} \cdot \nabla \varphi_p(\boldsymbol{x}) \ge 0$, and $\forall \boldsymbol{z} \in \mathbb{R}^m$, $\boldsymbol{z} \cdot \nabla \psi_p(\boldsymbol{z}) \ge 0$, we get

$$ \|\boldsymbol{x}^p(t)\|^2 + \|\boldsymbol{z}^p(t)\|^2 \le C \left( 1 + \int_0^t \|\boldsymbol{x}^p(s)\|^2 + \|\boldsymbol{z}^p(s)\|^2 \mathrm{d}s \right). $$

Thus, Gronwall inequality yields

$$ \sup_{s \le t} [\|\boldsymbol{x}^p(s)\|^2 + \|\boldsymbol{z}^p(s)\|^2] \le C \exp(Ct). $$

Substituting into (105) and using the fact that $\psi_p \ge 0$ gives (104). Next we can use the same arguments as in the proof of Theorem C.1 to show that $\boldsymbol{x}$ satisfies the differential inclusion involving $\varphi$. We finally discuss the case of $\boldsymbol{z}$. Proceeding as in the proof of Theorem C.1, we have $\forall \boldsymbol{v} \in C([0, T]; \mathbb{R}^n)$,

$$ \forall 0 \le t < t + h \le T, \quad \int_t^{t+h} \nabla \psi_p(\boldsymbol{z}^p(s)) \cdot (\boldsymbol{v}(s) - J_p \boldsymbol{z}^p(s)) + \psi(J_p \boldsymbol{z}^p(s)) \mathrm{d}s \le \int_t^{t+h} \psi(\boldsymbol{v}(s)) \mathrm{d}s. $$

We do not assume $\sup_{p \ge 1} \|\nabla \psi_p(\boldsymbol{z})\| < \infty$, so we cannot claim that $J_p \boldsymbol{z}^p \to \boldsymbol{z}$ in $C([0, T]; \mathbb{R}^m)$, in contrast with the proof of Theorem C.1. However, since $\nabla \psi_p(\boldsymbol{z}^p(s)) = p(\boldsymbol{z}^p(s) - J_p^z \boldsymbol{z}^p(s))$ and $\sup_p \int_0^T \|\nabla \psi_p(\boldsymbol{z}^p(s))\|^2 \mathrm{d}s < \infty$, we get

$$ \lim_{p \to \infty} \int_0^T \|\boldsymbol{z}^p(s) - J_p^z \boldsymbol{z}^p(s)\|^2 \mathrm{d}s = 0. $$

We can extract a subsequence $p_k$ to get

$$ \boldsymbol{z}^{p_k} - J_{p_k}^z \boldsymbol{z}^{p_k} \to 0, \text{ a.e in } (0, T). $$

Therefore

$$ J_{p_k}^z \boldsymbol{z}^{p_k} \to \boldsymbol{z}, \text{ a.e in } (0, T). $$

Observe that by the l.s.c. property of $\psi$, we have

$$ \psi(\boldsymbol{z}(t)) \le \liminf_{k \to \infty} \psi(J_{p_k}^z \boldsymbol{z}^{p_k}(t)), \text{ a.e in } (0, T). $$

Thus using Fatou's lemma, we obtain

$$ \int_t^{t+h} \psi(\boldsymbol{z}(s)) \mathrm{d}s \le \int_t^{t+h} \liminf_{k \to \infty} \psi(J_{p_k}^z \boldsymbol{z}^{p_k}(s)) \mathrm{d}s \le \liminf_{k \to \infty} \int_t^{t+h} \psi(J_{p_k}^z \boldsymbol{z}^{p_k}(s)) \mathrm{d}s. $$
Also, we have

$$ \int_t^{t+h} \nabla \psi_{p_k}(\boldsymbol{z}^{p_k}(s)) \cdot (\boldsymbol{v}(s) - J_{p_k}^z \boldsymbol{z}^{p_k}(s)) \mathrm{d}s = \int_t^{t+h} \nabla \psi_{p_k}(\boldsymbol{z}^{p_k}(s)) \cdot (\boldsymbol{v}(s) - \boldsymbol{z}(s)) \mathrm{d}s \\ + \int_t^{t+h} \nabla \psi_{p_k}(\boldsymbol{z}^{p_k}(s)) \cdot (\boldsymbol{z}(s) - J_{p_k}^z \boldsymbol{z}^{p_k}(s)) \mathrm{d}s. $$

In the right-hand side, as $k \uparrow \infty$, the first term goes to $\int_t^{t+h} \delta_z(s) \cdot (\boldsymbol{v}(s) - \boldsymbol{z}(s)) \mathrm{d}s$ because $\nabla \psi_{p_k}(\boldsymbol{z}^{p_k})$ weakly converges to $\delta_z$ and the second term goes to 0 because $\|\nabla \psi_{p_k}(\boldsymbol{z}^{p_k}(s))\|$ is uniformly bounded in $L^2(0, T)$ with respect to $k$ and $\boldsymbol{z} - J_{p_k}^z(\boldsymbol{z}^{p_k})$ goes to 0 in $L^2(0, T; \mathbb{R}^m)$. Therefore

$$ \forall 0 \le t < t + h \le T, \quad \int_t^{t+h} \delta_z(s) \cdot (\boldsymbol{v}(s) - \boldsymbol{z}(s)) + \psi(\boldsymbol{z}(s)) \mathrm{d}s \le \int_t^{t+h} \psi(\boldsymbol{v}(s)) \mathrm{d}s. $$

The proof is complete.

<figure><img src="image_12.png" /></figure>

# D Existence and uniqueness for (53)

Let us first take a look at the case where we remove the multivalued operator $\partial\varphi$ from the drift in (53). The problem becomes the same as (5) where $\sigma$ is constant and in particular it does not involve a stochastic integral. Thus, as pointed out in page 294 of [18], the proof of existence and uniqueness of a solution (still in [18], Theorem 2.9 page 289) can be simplified in a way that makes no use of probabilistic tools. We consider the Wiener space $(\Omega \triangleq C([0, T]; \mathbb{R}^d), \mathcal{F} \triangleq \mathcal{B}(C([0, T]; \mathbb{R}^d)), \mathbb{P})$ here $\Omega$ is the space of $\mathbb{R}^d$-valued continuous functions on $[0, T]$ endowed with the norm $\forall \boldsymbol{\omega} \in \Omega, \|\boldsymbol{\omega}\| \triangleq \sup_{0 \le t \le T} \|\boldsymbol{\omega}(t)\|, \mathcal{B}(C([0, T]; \mathbb{R}^d))$ is the Borel $\sigma$-algebra on $\Omega, \mathbb{P}$ is the Wiener measure;

the mappings indexed by $t \in [0, T]$, $W_t(\cdot): \Omega \to \mathbb{R}^d, \boldsymbol{\omega} \mapsto W_t(\boldsymbol{\omega}) \triangleq \boldsymbol{\omega}(t)$, the sequence of $\sigma$-algebras $\mathcal{F}^t \triangleq \sigma\{W_s, 0 \le s \le t\}$ and the map $\mathbf{X}: \Omega \to C([0, T]; \mathbb{R}^n), \boldsymbol{\omega} \mapsto \mathbf{X}(\boldsymbol{\omega}) \triangleq \mathbf{x}$ where $\forall 0 \le t \le T, \mathbf{x}(t) = \mathbf{x}(0) + \int_0^t \mathbf{b}(\mathbf{x}(s))\mathrm{d}s + \Gamma\boldsymbol{\omega}(t)$. Under $\mathbb{P}$, $\mathbf{W}$ is a Wiener process and $\mathbf{X}(\mathbf{W})$ solves (5) where $\sigma$ is constant. In this approach, the key ingredient is the mapping $\mathbf{X}$. For obtaining the existence and uniqueness of the solution to (53) with the multivalued operator $\partial\varphi$, we discuss below the properties of a similar mapping to $\mathbf{X}$ which involves the multivalued operator. This is done via the so-called “Generalized Skorokhod Problem”. The discussion follows [30] from page 245 to page 252. We use the notation $BV[0, T]$ for the space of functions with bounded variation on $[0, T]$.

**Definition D.1 (Generalized Convex Skorokhod Problem).** If a pair of functions $(\mathbf{x}, \mathbf{\Delta})$ satisfies the following conditions

1. $\mathbf{x}, \mathbf{\Delta}: [0, T] \to \mathbb{R}^n$ are continuous, $\mathbf{x}(0) = \mathbf{x}_0$ and $\mathbf{\Delta}(0) = 0$,

2. $\forall 0 \le t \le T, \mathbf{x}(t) \in \overline{\mathrm{Dom}(\partial\varphi)}, \mathbf{\Delta} \in BV([0, T]; \mathbb{R}^n)$,

3. $\forall 0 \le t \le T, \mathbf{x}(t) + \mathbf{\Delta}(t) = \mathbf{x}_0 + \int_0^t \mathbf{b}(\mathbf{x}(s))\mathrm{d}s + \Gamma\boldsymbol{\omega}(t)$,

4. $\forall 0 \le s \le t \le T, \forall \mathfrak{z} \in \mathbb{R}^n, \int_s^t (\mathbf{x}(r) - \mathfrak{z}) \cdot d\mathbf{\Delta}(r) + \int_s^t \varphi(\mathbf{x}(r))\mathrm{d}r \le (t - s)\varphi(\mathfrak{z})$,

then we say that $\mathbf{x}$ solves the generalized Skorokhod problem with parameters $\partial\varphi, \mathbf{x}_0, \mathbf{b}$ and $\boldsymbol{\omega}$ and we use the notation $\mathbf{x} = \mathcal{G}\mathcal{S}\mathcal{P}(\partial\varphi, \mathbf{x}_0, \mathbf{b}, \boldsymbol{\omega})$.
Existence and uniqueness of a solution for the Generalized Skorokhod Problem can be found in Theorem 4.17 page 252. This is obtained under the following conditions : $\varphi$ is a l.s.c convex function and int(Dom($\varphi$)) $\neq \emptyset$; $\mathbf{b}$ is Lipschitz, $x_0 \in \overline{\mathrm{Dom}}(\partial\varphi)$ and $\omega: [0, T] \mapsto \mathbb{R}^n$ is continuous with $\omega(0) = 0$. The continuity of the mapping $\mathbf{X}: \Omega \to C([0, T]; \mathbb{R}^n)$, $\omega \mapsto \mathbf{X}(\omega) \triangleq \mathbf{x}$ where $\mathbf{x} = \mathcal{G}\mathcal{S}\mathcal{P}(\partial\varphi, \mathbf{x}_0, \mathbf{b}, \omega)$ is shown in proposition 4.16 page 247. We use the notation $S_n^0[0, T]$ for the space of progressively measurable continuous stochastic processes (p.m.c.s.p.) from $\Omega \times [0, T]$ to $\mathbb{R}^n$,

$$S_n^2[0, T] \triangleq \left\{ \mathbf{Z} \in S_n^0[0, T], \mathbb{E}\left[ \sup_{0 \le t \le T} \|\mathbf{Z}(t)\|^2 \right] < \infty \right\}.$$

Within the framework of the aforementioned Wiener space, $\mathbf{X}(\mathbf{W}) \in S_n^0[0, T]$ solves (53). Furthermore, it can be shown that $\mathbf{X}(\mathbf{W}) \in S_n^2[0, T]$.

## E Proof of Proposition 5.1

*Proof of the first item of Proposition 5.1.* Using (47), we get

$$\begin{aligned} \frac{1}{2} \|\mathbf{X}_t^{\varepsilon,p} - \mathbf{X}_t^{\varepsilon,q}\|^2 &= \int_0^t (\mathbf{X}_s^{\varepsilon,p} - \mathbf{X}_s^{\varepsilon,q}) \cdot (\mathbf{b}(\mathbf{X}_s^{\varepsilon,p}) - \mathbf{b}(\mathbf{X}_s^{\varepsilon,q}))ds \\ &\quad - \int_0^t (\mathbf{X}_s^{\varepsilon,p} - \mathbf{X}_s^{\varepsilon,q}) \cdot (\nabla\varphi_p(\mathbf{X}_s^{\varepsilon,p}) - \nabla\varphi_q(\mathbf{X}_s^{\varepsilon,q}))ds \\ &\le C \int_0^t \|\mathbf{X}_s^{\varepsilon,p} - \mathbf{X}_s^{\varepsilon,q}\|^2 ds \\ &\quad + \left(\frac{1}{p} + \frac{1}{q}\right) \int_0^t \nabla\varphi_p(\mathbf{X}_s^{\varepsilon,p}) \cdot \nabla\varphi_q(\mathbf{X}_s^{\varepsilon,q})ds. \end{aligned}$$

Then under the condition (51) and from an application of Gronwall inequality, we obtain

$$\|\mathbf{X}_t^{\varepsilon,p} - \mathbf{X}_t^{\varepsilon,q}\|^2 \le C_{b,\varphi,t} \left(\frac{1}{p} + \frac{1}{q}\right).$$

Here the constant $C_{b,\varphi,t}$ depends on $b, \varphi$ and $t$. This implies the result. $\square$

<figure><img src="image_13.png" /></figure>

*Proof of the second item of Proposition 5.1.* Using the same arguments as in the proof of the first item, we get that

$$\forall p \ge 1, \mathbb{E} \left[ \sup_{t \le T} \| \mathbf{X}_t^{0,p} - \mathbf{X}_t^0 \|^2 \right] \le \frac{C_T}{p}, \quad (106)$$

where $\mathbf{X}^{0,p}$ is an approximation of $\mathbf{X}^0$ in the following sense:

$$d\mathbf{X}^{0,p} + \nabla\varphi_p(\mathbf{X}^{0,p})dt = \mathbf{b}(\mathbf{X}^{0,p})dt + \Gamma d\mathbf{W}_t.$$

Therefore, for any $\delta > 0$ and $p$,

$$\begin{aligned} \mathbb{P} \left( \sup_{t \le T} \| \mathbf{X}_t^{\varepsilon} - \mathbf{X}_t^0 \| \ge \delta \right) &\le \mathbb{P} \left( \sup_{t \le T} \| \mathbf{X}_t^{\varepsilon} - \mathbf{X}_t^{\varepsilon,p} \| \ge \frac{\delta}{3} \right) + \mathbb{P} \left( \sup_{t \le T} \| \mathbf{X}_t^{\varepsilon,p} - \mathbf{X}_t^{0,p} \| \ge \frac{\delta}{3} \right) \\ &\quad + \mathbb{P} \left( \sup_{t \le T} \| \mathbf{X}_t^{0,p} - \mathbf{X}_t^0 \| \ge \frac{\delta}{3} \right) \\ &\le \frac{18C_T}{p\delta^2} + \mathbb{P} \left( \sup_{t \le T} \| \mathbf{X}_t^{\varepsilon,p} - \mathbf{X}_t^{0,p} \| \ge \frac{\delta}{3} \right), \end{aligned}$$
by Markov inequality. From Proposition 3.5, we have

$$ \limsup_{\epsilon \to 0} \mathbb{P} \left( \sup_{t \in [0, T]} \| \mathbf{X}_t^{\epsilon} - \mathbf{X}_t^0 \| \ge \delta \right) \le \frac{18C_T}{p\delta^2}, $$

which holds for any $p$, hence the desired result. $\square$

## F Proof of Proposition 5.2

*Proof of the first item of Proposition 5.2.* We first note that $\psi$ satisfies (48). Using (47) for $\varphi$ and $\psi$, we get

$$ \begin{aligned} & \frac{1}{2} \|\mathbf{X}_t^{\epsilon,p} - \mathbf{X}_t^{\epsilon,q}\|^2 + \frac{1}{2} \|\mathbf{Z}_t^{\epsilon,p} - \mathbf{Z}_t^{\epsilon,q}\|^2 \\ &= \int_0^t (\mathbf{X}_s^{\epsilon,p} - \mathbf{X}_s^{\epsilon,q}) \cdot (\mathbf{b}^X(\mathbf{X}_s^{\epsilon,p}, \mathbf{Z}_s^{\epsilon,q}) - \mathbf{b}^X(\mathbf{X}_s^{\epsilon,q}, \mathbf{Z}_s^{\epsilon,q})) \mathrm{d}s \\ & \quad + \int_0^t (\mathbf{Z}_s^{\epsilon,p} - \mathbf{Z}_s^{\epsilon,q}) \cdot (\mathbf{b}^Z(\mathbf{X}_s^{\epsilon,p}, \mathbf{Z}_s^{\epsilon,q}) - \mathbf{b}^Z(\mathbf{X}_s^{\epsilon,q}, \mathbf{Z}_s^{\epsilon,q})) \mathrm{d}s \\ & \quad - \int_0^t (\mathbf{X}_s^{\epsilon,p} - \mathbf{X}_s^{\epsilon,q}) \cdot (\nabla\varphi_p(\mathbf{X}_s^{\epsilon,p}) - \nabla\varphi_q(\mathbf{X}_s^{\epsilon,q})) \mathrm{d}s \\ & \quad - \int_0^t (\mathbf{Z}_s^{\epsilon,p} - \mathbf{Z}_s^{\epsilon,q}) \cdot (\nabla\psi_p(\mathbf{Z}_s^{\epsilon,p}) - \nabla\psi_q(\mathbf{Z}_s^{\epsilon,q})) \mathrm{d}s \\ & \le C \int_0^t \{\|\mathbf{X}_s^{\epsilon,p} - \mathbf{X}_s^{\epsilon,q}\|^2 + \|\mathbf{Z}_s^{\epsilon,p} - \mathbf{Z}_s^{\epsilon,q}\|^2\} \mathrm{d}s \\ & \quad + \left(\frac{1}{p} + \frac{1}{q}\right) \left\{ \int_0^t \|\nabla\varphi_p(\mathbf{X}_s^{\epsilon,p})\| \|\nabla\varphi_q(\mathbf{X}_s^{\epsilon,q})\| \mathrm{d}s + \int_0^t \|\nabla\psi_p(\mathbf{Z}_s^{\epsilon,p})\| \|\nabla\psi_q(\mathbf{Z}_s^{\epsilon,q})\| \mathrm{d}s \right\}. \end{aligned} $$

From (51),

$$ \int_0^t \|\nabla\varphi_p(\mathbf{X}_s^{\epsilon,p})\| \|\nabla\varphi_q(\mathbf{X}_s^{\epsilon,q})\| \mathrm{d}s \le Ct. $$

From now on we focus on $\int_0^t \|\nabla\psi_p(\mathbf{Z}_s^{\epsilon,p})\| \|\nabla\psi_q(\mathbf{Z}_s^{\epsilon,q})\| \mathrm{d}s$. Let us proceed with the following expansion

$$ \psi_p(\mathbf{Z}_T^{\epsilon,p}) = \psi_p(\mathbf{z}_0) + \int_0^T \nabla\psi_p(\mathbf{Z}_s^{\epsilon,p}) \cdot (\mathbf{b}^Z(\mathbf{X}_s^{\epsilon,p}, \mathbf{Z}_s^{\epsilon,p}) - \nabla\psi_p(\mathbf{Z}_s^{\epsilon,p})) \mathrm{d}s, $$

which implies

$$ \psi_p(\mathbf{Z}_T^{\epsilon,p}) + \frac{1}{2} \int_0^T \|\nabla\psi_p(\mathbf{Z}_s^{\epsilon,p})\|^2 \mathrm{d}s \le \psi_p(\mathbf{z}_0) + \frac{1}{2} \int_0^T \|\mathbf{b}^Z(\mathbf{X}_s^{\epsilon,p}, \mathbf{Z}_s^{\epsilon,p})\|^2 \mathrm{d}s. $$

Therefore

$$ \psi_p(\mathbf{Z}_T^{\epsilon,p}) + \frac{1}{2} \int_0^T \|\nabla\psi_p(\mathbf{Z}_s^{\epsilon,p})\|^2 \mathrm{d}s \le \psi_p(\mathbf{z}_0) + \frac{C}{2} \int_0^T \{1 + \|\mathbf{X}_s^{\epsilon,p}\|^2 + \|\mathbf{Z}_s^{\epsilon,p}\|^2\} \mathrm{d}s. \quad (107) $$
We study the term in the integral of the right-hand side of the inequality (107):

$$
\begin{align}
\|\mathbf{X}_s^{\varepsilon,p}\|^2 + \|\mathbf{Z}_s^{\varepsilon,p}\|^2 &= \|\mathbf{X}_0^{\varepsilon,p}\|^2 + \|\mathbf{Z}_0^{\varepsilon,p}\|^2 + 2 \int_0^s \mathbf{X}_r^{\varepsilon,p} \cdot (\mathbf{b}^X(\mathbf{X}_r^{\varepsilon,p}, \mathbf{Z}_r^{\varepsilon,p}) - \nabla \phi_p(\mathbf{X}_r^{\varepsilon,p}) + \sigma \frac{\eta_r^\varepsilon}{\varepsilon}) \mathrm{d}r \nonumber \\
&\quad + 2 \int_0^s \mathbf{Z}_r^{\varepsilon,p} \cdot (\mathbf{b}^Z(\mathbf{X}_r^{\varepsilon,p}, \mathbf{Z}_r^{\varepsilon,p}) - \nabla \psi_p(\mathbf{Z}_r^{\varepsilon,p})) \mathrm{d}r \nonumber \\
&\le \|\mathbf{X}_0^{\varepsilon,p}\|^2 + \|\mathbf{Z}_0^{\varepsilon,p}\|^2 + 2C \int_0^s \{1 + \|\mathbf{X}_r^{\varepsilon,p}\|^2 + \|\mathbf{Z}_r^{\varepsilon,p}\|^2\} \mathrm{d}r \nonumber \\
&\quad + \int_0^s \mathbf{X}_r^{\varepsilon,p} \cdot \left(\sigma \frac{\eta_r^\varepsilon}{\varepsilon}\right) \mathrm{d}r, \tag{108}
\end{align}
$$

where we have used (46) to get the last inequality. We want to estimate the last term of the
right-hand side of (108). If we introduce the function $\phi(\mathbf{x}, \eta) = \mathbf{x} \cdot (\boldsymbol{\sigma}\mathbf{A}^{-1}\eta)$, then we get by
(85-86)

$$
\mathbb{E}[\varepsilon\phi(\mathbf{X}_s^{\varepsilon,p}, \eta_s^\varepsilon) - \varepsilon\phi(\mathbf{X}_0^{\varepsilon,p}, \eta_0^\varepsilon)] = -\mathbb{E}\left[\int_0^s \mathbf{X}_r^{\varepsilon,p} \cdot \left(\boldsymbol{\sigma} \frac{\eta_r^\varepsilon}{\varepsilon}\right) \mathrm{d}r\right] \\
+ \mathbb{E}\left[\int_0^s \left(\varepsilon \mathbf{b}^X(\mathbf{X}_r^{\varepsilon,p}, \mathbf{Z}_r^{\varepsilon,p}) - \varepsilon \nabla \varphi_p(\mathbf{X}_r^{\varepsilon,p}) + \boldsymbol{\sigma} \eta_r^\varepsilon\right) \cdot \left(\boldsymbol{\sigma} \mathbf{A}^{-1} \eta_r^\varepsilon\right) \mathrm{d}r\right].
$$

As $\eta^\epsilon$ is stationary:

$$
\mathbb{E} \left[ \int_{0}^{s} \left( \boldsymbol{\sigma} \eta_{r}^{\epsilon} \right) \cdot \left( \boldsymbol{\sigma} \mathbf{A}^{-1} \eta_{r}^{\epsilon} \right) \mathrm{d}r \right] = C_0 s, \quad C_0 = \mathbb{E} \left[ \left( \boldsymbol{\sigma} \eta_{0}^{1} \right) \cdot \left( \boldsymbol{\sigma} \mathbf{A}^{-1} \eta_{0}^{1} \right) \right].
$$

As $\nabla\varphi_p$ is bounded and $b^X$ is Lipschitz, we get

$$
\left| \mathbb{E} \left[ \int_0^s \mathbf{X}_r^{\varepsilon,p} \cdot \left( \boldsymbol{\sigma} \frac{\eta_r^\varepsilon}{\varepsilon} \right) \mathrm{d}r \right] \right| \le C_s + C_\varepsilon \left( \|\mathbf{X}_0^{\varepsilon,p}\| + \mathbb{E}[\|\mathbf{X}_s^{\varepsilon,p}\|^2]^{1/2} + \int_0^s \{\mathbb{E}[\|\mathbf{X}_r^{\varepsilon,p}\|^2] + \mathbb{E}[\|\mathbf{Z}_r^{\varepsilon,p}\|^2]\}^{1/2} \mathrm{d}r \right).
$$

Therefore, by substituting into (108), we can deduce that

$$
\sup_{s \le t} \mathbb{E} [\|\mathbf{X}_s^{\varepsilon,p}\|^2 + \|\mathbf{Z}_s^{\varepsilon,p}\|^2] \le C \left( 1 + \int_0^t \sup_{s \le r} \mathbb{E} [\|\mathbf{X}_s^{\varepsilon,p}\|^2 + \|\mathbf{Z}_s^{\varepsilon,p}\|^2] \, dr \right),
$$

which yields by Gronwall’s inequality

$$
\sup_{s \le t} \mathbb{E} \left[ \| \mathbf{X}_s^{\varepsilon,p} \|^2 + \| \mathbf{Z}_s^{\varepsilon,p} \|^2 \right] \le C \exp(Ct).
$$

The constant $C$ does not depend on $\varepsilon, p$. Substituting into (107) and using (52) and $\psi_p \ge 0$ (by (45)) gives

$$
\sup_{p \le \varepsilon} \sup \mathbb{E} \left[ \int_0^T \| \nabla \psi_p(\mathbf{Z}_s^{\varepsilon,p}) \|^2 \mathrm{d}s \right] < \infty. \quad (109)
$$

Finally, combining inequalities above, we obtain

$$
\mathbb{E} \left[ \sup_{s \le t} \left\{ \|\mathbf{X}_s^{\varepsilon,p} - \mathbf{X}_s^{\varepsilon,q}\|^2 + \|\mathbf{Z}_s^{\varepsilon,p} - \mathbf{Z}_s^{\varepsilon,q}\|^2 \right\} \right] \le C \int_0^t \mathbb{E} \left[ \sup_{s \le r} \left\{ \|\mathbf{X}_s^{\varepsilon,p} - \mathbf{X}_s^{\varepsilon,q}\|^2 + \|\mathbf{Z}_s^{\varepsilon,p} - \mathbf{Z}_s^{\varepsilon,q}\|^2 \right\} \right] \mathrm{d}r \\
+ C \left( \frac{1}{p} + \frac{1}{q} \right),
$$
which in turn provides

$$ \mathbb{E} \left[ \sup_{s \le t} \left\{ \| \mathbf{X}_s^{\varepsilon,p} - \mathbf{X}_s^{\varepsilon,q} \|^2 + \| \mathbf{Z}_s^{\varepsilon,p} - \mathbf{Z}_s^{\varepsilon,q} \|^2 \right\} \right] \le C \left( \frac{1}{p} + \frac{1}{q} \right). $$

The constant C does not depend on $\varepsilon, p, q$. The proof is complete.

*Proof of the second item of Proposition 5.2.*

Using the same arguments as in the proof of Proposition 5.1 (note that $\mathbb{E}[\int_0^s \mathbf{X}_r^{0,p} \cdot (\Gamma d\mathbf{W}_r)] = 0$), we get that

$$ \forall p \ge 1, \quad \mathbb{E} \left[ \sup_{t \le T} \left\{ \| \mathbf{Z}_t^{0,p} - \mathbf{X}_t^0 \|^2 + \| \mathbf{Z}_t^{0,p} - \mathbf{Z}_t^0 \|^2 \right\} \right] \le \frac{C_T}{p}, \quad (110) $$

where $(\mathbf{X}^{0,p}, \mathbf{Z}^{0,p})$ is an approximation of $(\mathbf{X}^0, \mathbf{Z}^0)$ in the following sense

$$ \begin{cases} \mathrm{d}\mathbf{X}^{0,p} + \nabla\varphi_p(\mathbf{X}^{0,p})\mathrm{d}t = \mathbf{b}^X(\mathbf{X}^{0,p}, \mathbf{Z}^{0,p})\mathrm{d}t + \Gamma\mathrm{d}\mathbf{W}_t, \\ \mathrm{d}\mathbf{Z}^{0,p} + \nabla\psi_p(\mathbf{Z}^{0,p})\mathrm{d}t = \mathbf{b}^Z(\mathbf{X}^{0,p}, \mathbf{Z}^{0,p})\mathrm{d}t. \end{cases} \quad (111) $$

The proof is then similar as the one of Proposition 5.1.

# G Proof of Lemma 6.1

Let $\mathbf{X}^\varepsilon$ satisfy (30). We address the evaluation of $\mathbb{E}[(f(\mathbf{X}_t^\varepsilon) - f(\mathbf{X}_t^0))^2]$. For any $p > 0$, by using Eqs. (54) and (106), we have

$$ \begin{aligned} \mathbb{E}[(f(\mathbf{X}_t^\varepsilon) - f(\mathbf{X}_t^0))^2] & \le 4\mathbb{E}[(f(\mathbf{X}_t^\varepsilon) - f(\mathbf{X}_t^{0,p}))^2] + 4\mathbb{E}[(f(\mathbf{X}_t^{0,p}) - f(\mathbf{X}_t^0))^2] \\ & \le \frac{8C_T \|\nabla f\|_\infty^2}{p} + 4\mathbb{E}[(f(\mathbf{X}_t^{0,p}) - f(\mathbf{X}_t^0))^2]. \end{aligned} $$

In order to get an estimate of the last term, we can follow the steps of the proof of Lemma 3.8 in the same way, because $\nabla\phi_p$, that appears only in $\Lambda_1$, is bounded uniformly in $p$. We get

$$ \begin{aligned} \mathbb{E}[(f(\mathbf{X}_t^\varepsilon) - f(\mathbf{X}_t^0))^2] & \le 4\mathbb{E}[(f(\mathbf{X}_t^\varepsilon) - f(\mathbf{X}_t^{0,p}))^2] + 4\mathbb{E}[(f(\mathbf{X}_t^{0,p}) - f(\mathbf{X}_t^0))^2] \\ & \le \frac{8C_T \|\nabla f\|_\infty^2}{p} + 4C\varepsilon^2. \end{aligned} $$

As this holds true for any $p$, this gives the desired result.

Let $(\mathbf{X}^\varepsilon, \mathbf{Z}^\varepsilon)$ satisfy (31). We address the evaluation of $\mathbb{E}[(f(\mathbf{X}_t^\varepsilon, \mathbf{Z}_t^\varepsilon) - f(\mathbf{X}_t^0, \mathbf{Z}_t^0))^2]$. For any $p$, by using (56) and (110), we have

$$ \begin{aligned} \mathbb{E}[(f(\mathbf{X}_t^\varepsilon, \mathbf{Z}_t^\varepsilon) - f(\mathbf{X}_t^0, \mathbf{Z}_t^0))^2] & \le 4\mathbb{E}[(f(\mathbf{X}_t^\varepsilon, \mathbf{Z}_t^\varepsilon) - f(\mathbf{X}_t^{0,p}, \mathbf{Z}_t^{0,p}))^2] \\ & + 4\mathbb{E}[(f(\mathbf{X}_t^{0,p}, \mathbf{Z}_t^{0,p}) - f(\mathbf{X}_t^0, \mathbf{Z}_t^0))^2] \\ & + 4\mathbb{E}[(f(\mathbf{X}_t^{0,p}, \mathbf{Z}_t^{0,p}) - f(\mathbf{X}_t^0, \mathbf{Z}_t^0))^2] \\ & \le \frac{8C_T \|\nabla f\|_\infty^2}{p} + 4\mathbb{E}[(f(\mathbf{X}_t^{0,p}, \mathbf{Z}_t^{0,p}) - f(\mathbf{X}_t^0, \mathbf{Z}_t^0))^2]. \end{aligned} $$
In order to get an estimate of the last term, we can follow the steps of the proof of Lemma 3.8 by keeping track of the bound (48) on $\nabla \psi_p$ (that appears only in $\Lambda_1$), and we get

$$ \mathbb{E}[(f(\mathbf{X}_t^\varepsilon, \mathbf{Z}_t^\varepsilon) - f(\mathbf{X}_t^0, \mathbf{Z}_t^0))^2] \le \frac{8C_T \|\nabla f\|_{\infty}^2}{p} + 4C(1+p)\varepsilon^2. $$

By optimizing in $p$ we get the desired result.

## H Proof of (74)

We consider the case of an OU noise that satisfies the equation $d\eta^\varepsilon = -\varepsilon^{-2}\eta^\varepsilon dt + \varepsilon^{-1}dW$ and $\eta_0 \sim \mathcal{N}(0, 1/2)$. The idea remains the same for the case of a Langevin noise. Define $X^\varepsilon_\star$ and $X^0_\star$ as follows:

$$ X_{\star,t}^{\varepsilon} \triangleq x_0 + \frac{1}{\varepsilon} \int_{0}^{t} \eta_{s}^{\varepsilon} ds \quad \text{and} \quad X_{\star,t}^{0} \triangleq x_0 + W_t. \qquad (112) $$

We first show that

$$ \mathbb{E}\left[\sup_{0 \le t \le T} |X_{\star,t}^{\varepsilon} - X_{\star,t}^{0}|^2\right] = O(\varepsilon^2 |\log \varepsilon|). \qquad (113) $$

From the equation for $\eta^\varepsilon$, it can be seen that $X^\varepsilon_{\star,t} = X^0_{\star,t} + \varepsilon(\eta_0 - \eta^\varepsilon_t)$. Thus, we have

$$ \mathbb{E}\left[\sup_{0 \le t \le T} |X_{\star,t}^{\varepsilon} - X_{\star,t}^{0}|^2\right] = \varepsilon^2 \mathbb{E}\left[\sup_{0 \le t \le T} |\eta_0 - \eta^\varepsilon_t|^2\right] \le \varepsilon^2 + 2\varepsilon^2 \mathbb{E}\left[\sup_{0 \le t \le T} |\eta^\varepsilon_t|^2\right]. $$

The process $\eta^\varepsilon_t$ is a stationary centered Gaussian process with covariance function $\mathbb{E}[\eta^\varepsilon_t \eta^\varepsilon_{t'}] = (1/2) \exp(-|t - t'|/\varepsilon^2)$. By the maximal inequality for the OU process [14] we get $\mathbb{E}[\sup_{0 \le t \le T} |\eta^\varepsilon_t|] \le C\sqrt{\log(1 + T/\varepsilon^2)}$ and by [26, Proposition 3.19] we obtain

$$ \mathbb{E}\left[\sup_{0 \le t \le T} |X_{\star,t}^{\varepsilon} - X_{\star,t}^{0}|^2\right] \le \varepsilon^2 + \varepsilon^2 + C^2 \varepsilon^2 \log(1 + T/\varepsilon^2) \le C' \varepsilon^2 (1 + |\log \varepsilon|). $$

Next, we show that

$$ \mathbb{E}\left[\sup_{0 \le t \le T} |X_t^\varepsilon - X_t^0|^2\right] = O(\varepsilon^2 |\log \varepsilon|). \qquad (114) $$

We can use an explicit formula for $X^\varepsilon$ (resp. $X^0$) that involves $X^\varepsilon_\star$ (resp. $X^0_\star$). Indeed, $X^\varepsilon_t = \mathcal{M}_t(X^\varepsilon_\star)$ and $X^0_t = \mathcal{M}_t(X^0_\star)$ where $\mathcal{M}$ is the self map on the set of continuous functions defined by $\mathcal{M}_t(f) \triangleq f(t) - \min_{0 \le s \le t} \min(0, f(s))$. This leads to

$$ |X_t^\varepsilon - X_t^0| \le |X_{\star,t}^\varepsilon - X_{\star,t}^0| + \left| \min_{0 \le s \le t} \min(0, X_{\star,s}^\varepsilon) - \min_{0 \le s \le t} \min(0, X_{\star,s}^0) \right|. $$

The second term in the right-hand side can be bounded by using the following inequalities i) $|\min(0, a) - \min(0, b)| \le |a - b|$ for all $a, b$ and ii) $|\min_{0 \le s \le t} f(s) - \min_{0 \le s \le t} g(s)| \le \max_{0 \le s \le t} |f(s) - g(s)|$, for all continuous functions $f, g$. Therefore,

$$ \mathbb{E}\left[\sup_{0 \le t \le T} |X_t^\varepsilon - X_t^0|^2\right] \le 4\mathbb{E}\left[\sup_{0 \le t \le T} |X_{\star,t}^\varepsilon - X_{\star,t}^0|^2\right] = O(\varepsilon^2 |\log \varepsilon|), $$

which gives (74).
## Acknowledgments.

LM expresses his sincere gratitude to Prof. Jean Michel Coron and the Sino-French International Associated Laboratory for Applied Mathematics for being supported for travels and housing at Ecole Polytechnique.

## References

[1] A. Bensoussan, C. Féau, L. Mertz, and S. C. P. Yam, An analytical approach for the growth rate of the variance of the deformation related to an elasto-plastic oscillator excited by a white noise, Appl. Math. Res. Express 2015 (2015), 99-128. 2

[2] A. Bensoussan, L. Mertz, and S. C. P. Yam, Long cycle behavior of the plastic deformation of an elasto-perfectly-plastic oscillator with noise, C. R. Math. Acad. Sci. Paris 350 (2012), 853-859. 2

[3] F. Bernardin, Equations différentielles multivoques : aspects théoriques et numériques - Applications, PhD thesis, Université Claude Bernard - Lyon I, 2004. 19

[4] H. Brézis, Opérateurs maximaux monotones et semi-groupes de contractions dans les espaces de Hilbert. (French) North-Holland Mathematics Studies, No. 5. Notas de Matemática (50), North-Holland, Amsterdam-London; Elsevier, New York, 1973. 7, 16, 17

[5] A. Cordle and J. Jonkman, State of the art design tools for floating offshore wind turbines, 21st International Offshore and Polar Engineering Conference, Maui, Hawaii, June 19-24, 2011. available at https://www.nrel.gov/docs/fy12osti/50543.pdf 1

[6] W. E. D. Liu, and E. Vanden-Eijnden, Analysis of multiscale methods for stochastic differential equations, Commun. Pure Appl. Math. 58 (2005), 1544-1585. 1, 2

[7] C. Féau, M. Laurière, and L. Mertz, Asymptotic formulae for the risk of failure related to an elasto-plastic problem with noise, Asymptot. Anal. 106 (2018), 47-60. 2

[8] J.-P. Fouque, J. Garnier, G. Papanicolaou, and K. Sølna, Wave Propagation and Time Reversal in Randomly Layered Media, Springer, New York, 2007. 1, 2, 4, 11, 27, 29

[9] J. Garnier and L. Mertz. Online extension of the present paper. arXiv:1906.01225

[10] M. B. Giles, Multilevel Monte Carlo methods, Acta Numer. 24 (2015), 259-328. 24

[11] D. Givon, I. G. Kevrekidis, and R. Kupferman, Strong convergence of projective integration schemes for singularly perturbed stochastic differential systems, Commun. Math. Sci. 4 (2006), 707-729. 2

[12] P. Glasserman, Monte Carlo Methods in Financial Engineering, Springer, New York, 2004. 1, 5, 6

[13] J. B. Goodman and K. K. Lin, Coupling control variates for Markov chain Monte Carlo, J. Comput. Phys. 228 (2009), 7127-7136. 2
[14] S. E. Graversen and G. Peskir, Maximal inequalities for the Ornstein-Uhlenbeck process, Proc. Amer. Math. Soc. 128 (2000), 3035-3041. 22, 42

[15] K. Hasselmann, T. Barnett, E. Bouws, H. Carlson, D. Cartwright, K. Enke, J. Ewing, H. Gien-app, D. Hasselmann, P. Kruseman, et al. Measurements of wind-wave growth and swell decay during the Joint North Sea Wave Project (JONSWAP), Tech. rep., Deutsches Hydrographisches Institut, 1973. 1

[16] J. M. Huang, J.-Q. Zhong, J. Zhang, and L. Mertz, Stochastic dynamics of fluid-structure interaction in turbulent thermal convection, J. Fluid Mech. 854 (2018), R5. 2

[17] S. Kameshwar and J. E. Padgett, Storm surge fragility assessment of above ground storage tanks, Structural Safety 70 (2018), 48-58. 1

[18] I. Karatzas and S. E. Shreve, Brownian motion and stochastic calculus, Springer-Verlag, New York, 1991. 37

[19] Y. Kifer, Stochastic versions of Anosov's and Neistadt's theorems on averaging, Stoch. Dyn. 1 (2001), 1-21. 2

[20] P. E. Kloeden and E. Platen, Numerical Solution of Stochastic Differential Equations, Springer, Berlin, 1992. 4, 5, 11

[21] W. Kohler and G. Papanicolaou, Asymptotic analysis of mixing stochastic ordinary differential equations, Comm. Pure and Appl. Math. 27 (1974), 641-668. 1

[22] H. J. Kushner, Weak Convergence Methods and Singularly Perturbed Stochastic Control and Filtering Problems, Birkhäuser, Boston, 1990. 11

[23] M. Laurière and L. Mertz, Penalization of non-smooth dynamical systems with noise: ergodicity and asymptotic formulae for threshold crossings probabilities, SIAM J. Appl. Dyn. Syst. 18 (2019), 853-880. 2

[24] Y. K. Lin and Y. Yong, Evolutionary Kanai-Tajimi earthquake models, J. Eng. Mech. 113 (1987), 1119-1137. 3, 10

[25] D. Liu, Strong convergence of principle of averaging for multiscale stochastic dynamical systems, Commun. Math. Sci. 8 (2010), 999-1020. 2

[26] P. Massart, Concentration inequalities and model selection, Ecole d'Eté de Probabilités de Saint-Flour XXXIII, Springer, Berlin, 2007. 22, 42

[27] L. Mertz, Stochastic variational inequalities for random mechanics, CIRM, Audiovisual resource, doi:10.24350/CIRM.V.19217703, available at http://dx.doi.org/10.24350/CIRM.V.19217703 20

[28] L. Mertz, G. Stadler, and J. Wylie, A backward Kolmogorov equation approach to compute means, moments and correlations of non-smooth stochastic dynamical systems, Physica D 397 (2019), 25-38. 20, 22, 24, 25
[29] MIT open courseware, https://ocw.mit.edu/courses/physics/8-01sc-classical-mechanics-fall-2016/week-2-newtons-laws/6.1-contact-forces/
https://ocw.mit.edu/courses/physics/8-01sc-classical-mechanics-fall-2016/week-2-newtons-laws/6.2-static-friction-lesson/ 19

[30] E. Pardoux and A. Răşcanu, Stochastic Differential Equations, Backward SDEs, Partial Differential Equations, Springer, Cham, 2014. 16, 37

[31] A. Patil, S. Jung, and O.-S. Kown, Structural performance of a parked wind turbine tower subjected to strong ground motions, Eng. Struct. 120 (2016), 92–102. 1

[32] G. A. Pavliotis and A. M. Stuart, Multiscale Methods Averaging and Homogenization, Springer, New York, 2008. 4

[33] A. Quilligan, A. O’Connor, and V. Pakrashi, Fragility analysis of steel and concrete wind turbine towers, Eng. Struct. 36 (2012), 270–282. 1

[34] J. B. Roberts and P. D. Spanos, Random Vibration and Statistical Linearization, Dover Publications, New York, 2003. 8, 14

[35] O. Vasicek, An equilibrium characterisation of the term structure, J. Financ. Econ. 5 (1977), 177–188. 3

[36] E. Wong and M. Zakai, On the convergence of ordinary integrals to stochastic integrals, Ann. Math. Statist. 36 (1965), 1560–1564. 4
This figure "TDOF.png" is available in "png" format from:

http://arxiv.org/ps/1906.01225v2