Jimmy Lominé
email: jlomine@rockwellcollins.com

Christophe Morlaas
email: morlaas@recherche.enac.fr

Hervé Aubert
email: haubert@laas.fr

Method for vector sensor design based on a spherical mode approach for 3D DoA estimation

Keywords: vector sensor, 3D direction of arrival finding, spherical mode, efficiency criterion

published or not. The documents may come

Method for vector sensor design based on a spherical mode approach for 3D DoA estimation

Jimmy Lominé, Christophe Morlaas, Hervé Aubert

I. INTRODUCTION

To determine the direction of arrival (DoA) of incoming electromagnetic signals, two parameters can be used: the spatial diversity and/or the polarization diversity. The major part of direction finding algorithms uses the spatial diversity such as interferometry which extracts the phase differential between each radiating element of the direction finding antenna in order to compute the DoA. It can be applied for 3D direction finding but the antenna design becomes too complex and cumbersome to compute the DoA in all space. The use of the polarization diversity has the advantage to require smaller sensor size. The DoA is estimated through the components measurement of the incoming electromagnetic field. It can be performed with 6 co-located antennas [START_REF] Hatke | Conditions for unambiguous source location using polarization diverse arrays[END_REF]. For instance, a combination of 3 elementary orthogonal electric dipoles and 3 elementary orthogonal magnetic dipoles can measure the six components of the electromagnetic field [START_REF] Nehorai | Vector-sensor array processing for electromagnetic source localization[END_REF]. Generally, to reduce the coupling effects, a spatial separation of antenna elements can be applied. A method has been proposed to estimate the DoA with a vector sensor with distributed elements using a spherical mode approach [START_REF] Chabory | Efficiency characterisation of vector-sensor antennas with distributed elements for 3d direction finding[END_REF]. For co-located antenna elements, only the 6 first spherical harmonics are necessary. But when a spatial separation is applied between antenna elements, higher order spherical modes radiate and thus the DoA estimation with only 6 antenna elements becomes inaccurate.

In this article, a criterion is firstly introduced to evaluate the DoA estimation efficiency. Secondly, with a vector sensor of 6 theoretical elements, the efficiency criterion is related to the accuracy of the DoA according to the spatial distribution of the antenna elements. Thus it allows a rapid determination of estimation performances trend. Finally, the previous results lead to study a vector sensor composed of 3 co-located elements by considering the radiation perturbations by means of a full wave electromagnetic simulation (FEKO). The aim of these numerical experimentations is to check that the efficiency criterion considers the effects of the radiation pertubations on the DoA accuracy. To achieve it, different configurations of the realistic antenna are simulated: one with an infinite metallic plate and others with a finite plate of different sizes.

II. PRINCIPLE

The number of significant spherical harmonics for accurate estimation can be related to the number of antenna elements [START_REF] Chabory | Reduction of the sensor number in distributed vector-antennas for 3d direction finding[END_REF]. But when the number of relevant spherical modes becomes higher than the 𝑁 𝑎 antenna elements due to a spatial separation applied to antenna elements or radiation perturbations, the DoA estimation becomes inaccurate. Thus it is important to evaluate the trend of the DoA estimation performances according to different vector sensor configurations and to the surrounding perturbations.

A. DoA estimation

According to the method presented by Chabory et al. [START_REF] Chabory | Efficiency characterisation of vector-sensor antennas with distributed elements for 3d direction finding[END_REF], the Poynting vector ⃗ 𝑃 (1) can be estimated from an inversion process (2) between the spherical harmonic decomposition of the outgoing electromagnetic fields and the incoming electromagnetic field computed at the center of sensor from the spherical harmonic recomposition.

⃗ 𝑃 = 1 2 𝑅𝑒{ ⃗ 𝐸 0 ∧ ⃗ 𝐻 0 * } (1) [ ⃗ 𝐸 0 ⃗ 𝐻 0 ] = [ 𝑒 𝑠- 0 , ℎ 𝑠- 0 ] (𝑀 + 𝑇 ) -1 𝐴 - (2)
with ( ⃗ 𝐸 0 , ⃗ 𝐻 0 ) the recomposed electromagnetic field at the center of the sensor, (𝑒 𝑠- 0 , ℎ 𝑠- 0 ) the electromagnetic fields at the center of the sensor of incoming spherical harmonics, 𝐴 - the vector of the wave amplitudes received at the antenna ports and 𝑀 + the matrix of the coefficients associated with the outgoing spherical harmonics 𝑠 + 𝑚,𝑛 [START_REF] Harrington | Time-Harmonic Electromagnetic Fields[END_REF], with 𝑚 the azimuth index and 𝑛 the elevation index.

𝑛 ∈ [1, ..., 𝑁 𝑠 ] 𝑚 ∈ [-𝑁 𝑠 , ..., 𝑁 𝑠 ] (3)
where 𝑁 𝑠 is the spherical mode order.

B. Criterion of the estimation efficiency

Thanks to the reaction theorem, the knowledge of 𝑀 + allows to predict the trend of the DoA estimation performances of a vector sensor. An efficiency criterion 𝜂𝑒 (5) is introduced from the singular values of 𝑀 + (4). Thus coupling effects and the surrounding perturbations for a realistic antenna can be considered.

𝑀 + = 𝑈 Σ𝑉 = [ 𝑢 1 . . . 𝑢 𝑁𝑠 ] ⎡ ⎢ ⎣ 𝜎 1 0 . . . 0 𝜎 𝑁𝑎 ⎤ ⎥ ⎦ ⎡ ⎢ ⎣ 𝑣 1 . . . 𝑣 𝑁𝑠 ⎤ ⎥ ⎦ (
𝜂𝑒 = Σ 𝑐𝑜𝑛𝑑(Σ) ( 5
)
where Σ and 𝑐𝑜𝑛𝑑(Σ) are the average of the singular values of 𝑀 + and its condition number respectively. To normalize the criterion 𝜂𝑒 according to another antenna configuration, a reference criterion of the estimation efficiency 𝜂𝑒 𝑅𝑒𝑓 can be used (6).

𝜂𝑒 𝑁 = 𝜂𝑒 𝜂𝑒 𝑅𝑒𝑓 (6) III. SIX DISTRIBUTED ANTENNA ELEMENTS (a) (b) (c) (d)
Fig. 1: Simulated antenna configurations with 3 elementary orthogonal electric dipoles and 3 elementary orthogonal magnetic dipoles.

Four geometrical configurations which are represented in Fig. 1 have been simulated. Each one is a combination of 3 elementary orthogonal electric dipoles and 3 elementary orthogonal magnetic dipoles spatially distributed. The elementary magnetic dipoles are considered from elementary electric loops. No coupling effects between antenna elements is considered. The reference configuration is the six co-located antenna elements. It is the ideal case with the spherical mode approach. The vector sensor performances are analyzed from the angular distance of the estimation error Δ𝑎 between estimated and real DoA as defined in (7). Δ𝑎 = cos -1 (cos(𝜃 𝑟 ) cos(𝜃 𝑒 ) + sin(𝜃 𝑟 ) sin(𝜃 𝑒 ) cos(𝜙 𝑟 -𝜙 𝑒 ))

(7) where (𝜙 𝑟 , 𝜃 𝑟 ) are the coordinates of the real DoA and (𝜙 𝑒 , 𝜃 𝑒 ) are the coordinates of estimated DoA. In Fig. 2 and in Fig. 3 the estimation efficiency (𝜂𝑒) and the estimation error (Δ 𝑎 ) are displayed according to the normalized separation distance 𝑑/𝜆 for the 3 antenna configurations represented in Fig. 1, with 𝜆 the wavelength. The incident plane waves have a circular polarization. The estimation error corresponds to the maximum estimation error in full 3D space. The criterion 𝜂𝑒 decreases according to the distance 𝑑/𝜆 due to the radiation of the higher order spherical modes and follows the behavior of the error estimation. The less the efficiency, the more the maximum DoA error. The most efficient sensor is the antenna configuration (𝑎) in Fig. 1. The ratio of the computation time between the simulation of the maximum DoA error and the computation of the estimation efficiency is 17 for an angular resolution of 1 ∘ . Thus it allows a rapid determination of the best antenna configuration.

IV. THREE CO-LOCATED ANTENNA ELEMENTS -{EZ, HX, HY}

The previous study of the efficiency criterion according to only the separation distance between antenna elements has shown that a lower distance of separation yields better performances. But the spatial distribution is not the only cause of DoA estimation errors. Indeed the coupling effects and the surrounding perturbations can disturb the radiation pattern of the antenna and so the direction finding accuracy. To check that these perturbations are taken into account through the efficiency criterion, a vector sensor of 3 co-located antenna elements is modeled. Only the half-space domain is considered by using a metallic plate (Fig. 4) in the aim of an easier manufacturing and to ensure a good reproducibility of DoA estimation whatever the antenna support. Full wave simulations are performed to take into account realistic antennas with the coupling effects and the metallic plate. Several antenna configurations are simulated in order to have a comparative study: one with an infinite plate and others with a finite plate with different sizes. The components of the electromagnetic field {Ex, Ey, Hz} are shorted due to the metallic plate at the center of the sensor. That's why only incident plane waves in TM mode have been simulated.

A. Adaptation of the DoA algorithm

The DoA algorithm must be adapted to a sensor of 3 antenna elements by judiciously selecting the relevant spherical harmonics in 𝑀 + . Indeed, with this spherical mode approach, only 𝑁 𝑎 spherical harmonics can be used to recompose the field at the center of the sensor.

Moreover the infinite metallic plate modifies the amplitudes distribution of the outgoing spherical harmonics coefficients as shown in Fig. 5 and in Fig. 6. Without metallic plate, Fig. 5: Theoretical amplitudes distribution of the outgoing spherical harmonics coefficients associated with the electric field without (blue) and with (green) an infinite plate.

only the first spherical mode radiates, 𝑁 𝑠 = 1. But with a plate, some higher spherical modes radiate. These higher Fig. 6: Theoretical amplitudes distribution of the outgoing spherical harmonics coefficients associated with the magnetic field without (blue) and with (green) an infinite plate. spherical harmonics allow only to cancel the electromagnetic field in the directions below the plate by combination with the fundamental mode. These modes are not necessary to recompose the electric and magnetic field [ ⃗ 𝐸 0 , ⃗ 𝐻 0 ] from the incoming waves received in the upper half-space domain. Thus 𝑀 + is adjusted by conserving the coefficients associated with the outgoing spherical harmonics 𝑠 + -1,1 and 𝑠 + 1,1 of the magnetic field and 𝑠 + 0,1 of the electric field.

Moreover the Poynting vector computation can not be performed without all field components. The DoA computation must be also adapted as defined in (8).

𝜃 = sin -1 ⎛ ⎝ |𝐸 𝑧 | √ (𝐻 2 𝑥 + 𝐻 2 𝑦 )𝜁 ⎞ ⎠ 𝜙 = tan -1 ( 𝐻 𝑦 𝐻 𝑥 ) (8)
with 𝜁 the wave impedance.

This adjustment of 𝑀 + has been validated by numerical simulations for elementary dipoles with an infinite metallic plate. No DoA error is generated by removing these higher spherical modes because all necessary information for an accurate DoA estimation are carried by the first spherical mode. This antenne will be the reference for the following study.

B. Realistic antenna design

The realistic vector sensor, represented in Fig. 7, is composed of two orthogonal and horizontal magnetic dipoles and a virtual vertical electric dipoles resulting from the combination of the two loops response by setting [1, 1, 1, 1] to the four ports. The operating frequency is set to 1GHz.

C. Full wave simulations

The purpose of these simulations is to check that the computation of the efficiency criterion takes into account the radiation perturbations and the metallic plate effect on the Fig. 7: Simulated antenna with a full wave electromagnetic simulation software. DoA accuracy. The fields radiated by the realistic antenna elements and simulated from full wave electromagnetic simulations are thus used to compute 𝑀 + . Only the half of the upper half-space domain is simulated because there is a symmetry in the results between [0 -180 ∘ ] and [180 -360 ∘ ] azimuth angles. The reference configuration used to normalize the criterion of the estimation efficiency is the theoretical sensor {Ez, Hx, Hy} with an infinite metallic plate previously studied in part IV-A. For each simulation, the DoA accuracy is computed separately for the azimuth and the elevation angles as shown for two configurations in Fig. 8 to 11. Then for each one, the efficiency criterion is computed and compared to the DoA estimation accuracy in Table I. The trend of DoA accuracy is well evaluated with the criterion of the estimation efficiency. It shows that the direction finding performances are sensitive to the size of the metallic plate. Indeed, the larger the plate, the better the estimation accuracy.

However with the infinite metallic plate the DoA error is degraded in comparison to the theoretical case with a maximum error of 12 ∘ in elevation. This DoA error is caused by the realistic antennas which generate radiation patterns different to elementary dipoles. This difference can be observed in the amplitudes distribution of the outgoing spherical harmonics coefficients for both configurations as shown in Fig. 12 and in Fig. 13. The coefficients amplitude associated with the higher spherical harmonics are modified. Thus with a realistic vector sensor, the removal of higher spherical modes for the adjustment of 𝑀 + prevents to get back all useful information for an accurate DoA estimation. But regarding the azimuth angles, the estimation is accurate because the ratio between the coefficients of the outgoing spherical harmonics associated with 𝐻 𝑥 and 𝐻 𝑦 and their orthogonality are respected. For the configurations with a finite metallic plate, this phenomenon is accentuated as shown in Fig. 14 and Fig. 15 for a plate of 1.5𝜆. This is the reason of the accuracy degradation increase.

V. CONCLUSION

An efficiency criterion of 3D DoA estimation from a singular values decomposition has been introduced to characterize the trend of the estimation accuracy. This criterion has been firstly evaluated according to different spatial distributions of antenna elements with a theoretical vector sensor of 6 colocated elementary dipoles. The results have confirmed that the efficiency parameter is related to the DoA estimation accuracy. For an angular resolution of 1 ∘ , the ratio of the computation time is 17 between the DoA algorithm with the spherical mode approach and the computation of the efficiency criterion. Thus this method of efficiency evaluation can be used as a design tool to allow a rapid determination of the most efficient antenna configuration among several sensors. Secondly, a vector sensor of 3 co-located realistic dipoles with different sizes of metallic plate has been simulated. These full wave simulations have allowed to validate that the computation of the estimation efficiency takes into account the radiation perturbations. However the 3D DoA estimation accuracy of a realistic vector sensor with this spherical mode approach is sensitive to the antenna elements radiation and to the surrounding perturbations. But for 2D direction finding, the sensor retains good performances. In future work, the prototype of the sensor {Ez, Hx, Hy} will be developed to validate the computation of the criterion efficiency in a real environment. According to the measured DoA accuracy, an optimization of the vector sensor design could be performed to improve its performances.

  [START_REF] Chabory | Reduction of the sensor number in distributed vector-antennas for 3d direction finding[END_REF] where 𝑢 𝑖 and 𝑣 𝑖 correspond to the 𝑖-th left and right singular vectors respectively, and 𝜎 1 ⩾ 𝜎 2 ⩾ . . . ⩾ 𝜎 𝑁𝑎 are the singular values of 𝑀 + .

Fig. 2 :

 2 Fig. 2: Efficiency of DoA estimation normalized with respect to the co-located configuration.

Fig. 3 :

 3 Fig. 3: Maximum DoA estimation error -(Δ𝑎) 𝑚𝑎𝑥 .

Fig. 4 :

 4 Fig. 4: Antenna configuration -{Ez, Hx, Hy}.

Fig. 8 :

 8 Fig. 8: DoA estimation error in elevation with an infinite metallic plate.

Fig. 9 :

 9 Fig. 9: DoA estimation error in azimuth with an infinite metallic plate.

Fig. 10 :

 10 Fig. 10: DoA estimation error in elevation with a finite metallic plate of 1.5𝜆.

Fig. 11 :

 11 Fig. 11: DoA estimation error in azimuth with a finite metallic plate of 1.5𝜆.

Fig. 12 :

 12 Fig. 12: Amplitudes distribution of the outgoing spherical harmonics coefficients associated with the electric field -Infinite plate: Theoretical (blue) and FEKO (green).

Fig. 13 :

 13 Fig. 13: Amplitudes distribution of the outgoing spherical harmonics coefficients associated with the magnetic field -Infinite plate: Theoretical (blue) and FEKO (green).

Fig. 14 :

 14 Fig. 14: Amplitudes distribution of the outgoing spherical harmonics coefficients associated with the electric field: Infinite plate with theoretical simulation (blue) and Finite plate of 1.5𝜆 with FEKO (green).

Fig. 15 :

 15 Fig. 15: Amplitudes distribution of the outgoing spherical harmonics coefficients associated with the magnetic field: Infinite plate with theoretical simulation (blue) and Finite plate of 1.5𝜆 with FEKO (green).

TABLE I :

 I Criterion of the estimation efficiency computed from full wave simulations according to maximal DoA error.

	Metallic plate size	𝜂𝑒 𝑁	Elevation error	Azimuth error
	Infinite	-6.05dB	12 ∘	0.3 ∘
	8𝜆	-6.323dB	27 ∘	1.7 ∘
	1.5𝜆	-6.340dB	30.6 ∘	2.8 ∘


Assyr Abdulle

Gilles Vilmart

Konstantinos C Zygalakis

Long time accuracy of Lie-Trotter splitting methods for Langevin dynamics

Keywords: stochastic differential equations, splitting method, Langevin dynamics, weak convergence, modified differential equations, backward error analysis, invariant measure, ergodicity AMS subject classifications. 65C30, 60H35, 37M25

A new characterization of sufficient conditions for the Lie-Trotter splitting to capture the numerical invariant measure of nonlinear ergodic Langevin dynamics up to an arbitrary order is discussed. Our characterization relies on backward error analysis and needs weaker assumptions than assumed so far in the literature. In particular, neither high weak order of the splitting scheme nor symplecticity are necessary to achieve high order approximation of the invariant measure of the Langevin dynamics. Numerical experiments confirm our theoretical findings.

Introduction

In many applications (see [START_REF] Milstein | Computing ergodic limits for Langevin equations[END_REF][START_REF] Leimkuhler | Rational construction of stochastic numerical methods for molecular sampling[END_REF][START_REF] Bou-Rabee | Pathwise accuracy and ergodicity of metropolized integrators for SDEs[END_REF][START_REF] Li | Effectiveness of implicit methods for stiff stochastic differential equations[END_REF]] and references therein), one is interested in simulating the invariant measure of a stochastic differential equation (SDE) by running a numerical scheme that approximates its time dynamics. In particular, one uses the numerical trajectories to construct an empirical measure either by averaging one single long trajectory or by averaging over many realisations to obtain a finite ensemble average (see for example [START_REF] Milstein | Computing ergodic limits for Langevin equations[END_REF]). One immediately is faced with the question regarding the quality of such an approximation. Another issue is the ergodicity of the numerical approximation itself. This second question is still an area of active research [START_REF] Mattingly | Ergodicity for SDEs and approximations: locally Lipschitz vector fields and degenerate noise[END_REF][START_REF] Talay | Stochastic Hamiltonian systems: exponential convergence to the invariant measure, and discretization by the implicit Euler scheme[END_REF][START_REF] Roberts | Exponential convergence of Langevin distributions and their discrete approximations[END_REF][START_REF] Shardlow | A perturbation theory for ergodic markov chains and application to numerical approximations[END_REF], but is not addressed here.

In this paper, we focus on the case of second order stochastic dynamics generated by the Langevin equation, of the form

dq(t) = M -1 p(t)dt, (1a)
dp(t) = (-∇V (q(t))γp(t))dt + 2β -1 γM 1/2 dW (t),

where the scalar function V : R d → R denotes the potential, the mass matrix M is symmetric positive definite, γ > 0 (the friction), β > 0 (related to the temperature) are fixed constants, p(t), q(t) ∈ R d , and W (t) denotes a standard d-dimensional Wiener process. For simplicity, we assume that the mass matrix is the identity matrix M = I, but we emphasize that our analysis would apply straightforwardly to general positive definite matrices M . Under appropriate smoothness and growth assumptions on the Hamiltonian energy (see e.g. [START_REF] Mattingly | Convergence of numerical time-averaging and stationary measures via Poisson equations[END_REF][START_REF] Kopec | Weak backward error analysis for Langevin process[END_REF]),

H(p, q) = 1 2 p T p + V (q), (2)
one can show that the dynamics generated by [START_REF] Abdulle | Improved stabilized Multilevel Monte Carlo method for stiff stochastic differential equations[END_REF] are ergodic, i.e. for arbitrary initial conditions q(0), p(0) (assumed deterministic for simplicity), and for all smooth test function φ ∈ C ∞ P (R 2d , R) with derivatives of all orders with polynomial growth, the time average of the trajectories of (1) satisfy with probability 1, lim

T →∞ 1 T T 0 φ(p(t), q(t))dt = R d ×R d φ(p, q)dµ(p, q).
The invariant measure µ is characterized by the Gibbs density function

ρ ∞ (p, q) = Ze -βH(p,q) , (3)
where Z is the normalization constant, which ensures that R d ×R d ρ ∞ (p, q)dpdq = 1. Similarly to [START_REF] Bou-Rabee | Long-run accuracy of variational integrators in the stochastic context[END_REF], we consider the following class of schemes based on the Lie-Trotter splitting,

X n+1 = Φ h • Θ h,n (X n ), (4)
where X n = (p n , q n ) T . The integrator Φ h approximates the exact flow with time h of the deterministic Hamiltonian part dq(t) = p(t)dt, dp(t) = -∇V (q(t))dt, [START_REF] Bou-Rabee | Pathwise accuracy and ergodicity of metropolized integrators for SDEs[END_REF] while Θ h,n is a numerical integrator for the stochastic part dq(t) = 0, (6a)

dp(t) = -γp(t)dt + 2β -1 γdW (t). (6b)
The exact solution of ( 6) is given by the variation of constants formula, q(t n+1 ) = q(t n ), p(t n+1 ) = e -γh p(t n ) + 2β -1 γ t n+1 tn e -γ(h-s) dW (s).

The trajectories of the above stochastic integral itself need not be simulated. Indeed, the properties of the stochastic integral [START_REF]Stochastic differential equations[END_REF] imply that the solution p(t) of (6) given p(0) is Gaussian and thus characterized by its mean and covariance. Using this, we obtain the following flow (q n+1 , p n+1 ) = Θ h,n (q n , p n ) for [START_REF] Burrage | Numerical methods for second-order stochastic differential equations[END_REF] which is exact in law, 1

q n+1 = q n , p n+1 = e -γh p n + β -1 (1 -e -2γh )ξ n ,
where ξ n ∼ N (0, I) are independent d-dimensional Gaussian random variables. The aim of this paper is to analyse splitting schemes such as (4) for a variety of different deterministic integrators Φ h used for solving [START_REF] Bou-Rabee | Pathwise accuracy and ergodicity of metropolized integrators for SDEs[END_REF], for the numerical approximation of the invariant measure, i.e. we quantify the error e(φ, h) := lim

N →∞ 1 N N n=1 φ(p n , q n ) - R d ×R d φ(p, q)ρ ∞ (p, q)dpdq, (7)
for all smooth test function φ ∈ C ∞ P (R 2d , R). We say that the method has order r for sampling the invariant measure if |e(φ, h)| ≤ Ch r for all h small enough, where C is independent of h, but depends on the smooth test function φ. In particular, using recent advances in the theory of modified integrators for SDEs and ODEs [START_REF] Abdulle | High order weak methods for stochastic differential equations based on modified equations[END_REF][START_REF] Chartier | Numerical integrators based on modified differential equations[END_REF] and backward error analysis for SDEs [START_REF] Zygalakis | On the existence and the applications of modified equations for stochastic differential equations[END_REF][START_REF] Debussche | Weak backward error analysis for SDEs[END_REF][START_REF] Kopec | Weak backward error analysis for Langevin process[END_REF][START_REF] Kopec | Weak backward error analysis for overdamped Langevin processes[END_REF][START_REF] Abdulle | High order numerical approximation of the invariant measure of ergodic SDEs[END_REF], we provide a new characterization of sufficient conditions for the Lie-Trotter splitting to capture the numerical invariant measure of nonlinear ergodic Langevin dynamics up to an arbitrary order.

It is well known that controlling the weak error |E(φ(X n )) -E(φ(X(t n ))| up to a given order allows to control the accuracy error |e(φ, h)| for the invariant measure up to the same order [START_REF] Talay | Expansion of the global error for numerical schemes solving stochastic differential equations[END_REF], [START_REF] Milstein | Weak approximation of solutions of systems of stochastic differential equations[END_REF]. However, the converse is not true, and new sufficient conditions independent of the weak order of accuracy for a numerical method to approximate with high order of accuracy the invariant measure of ergodic SDEs have recently been derived in [START_REF] Abdulle | High order numerical approximation of the invariant measure of ergodic SDEs[END_REF]. The accuracy of the numerical computation of the invariant measure for Langevin dynamics has recently been studied in a number of papers [START_REF] Bou-Rabee | Long-run accuracy of variational integrators in the stochastic context[END_REF][START_REF] Leimkuhler | Rational construction of stochastic numerical methods for molecular sampling[END_REF][START_REF] Leimkuhler | The computation of averages from equilibrium and nonequilibrium Langevin molecular dynamics[END_REF]. In [START_REF] Bou-Rabee | Long-run accuracy of variational integrators in the stochastic context[END_REF], the splitting method ( 4) is proposed and analyzed, while several other splitting methods with high order for the invariant measure are considered and studied in [START_REF] Leimkuhler | Rational construction of stochastic numerical methods for molecular sampling[END_REF][START_REF] Leimkuhler | The computation of averages from equilibrium and nonequilibrium Langevin molecular dynamics[END_REF]. It is shown in [START_REF] Bou-Rabee | Long-run accuracy of variational integrators in the stochastic context[END_REF] that if Φ h is a variational integrator (equivalently a symplectic scheme) which preserves the energy up to order r, then the Lie-Trotter splitting (4) has order r of accuracy for the invariant measure, although the standard weak error is only of order one. This is another example illustrating the fact that a high weak order of accuracy is not necessary to achieve high order for the invariant measure.

In this paper, we provide weaker assumptions to achieve order r in the approximation of the numerical invariant measure using the splitting (4). Our results show that to achieve order r for the invariant measure by the splitting (4), neither an order r of accuracy in the weak sense nor an order r of accuracy in the deterministic sense for Φ h in (4) are needed. This is also related to the results in [START_REF] Burrage | Numerical methods for second-order stochastic differential equations[END_REF][START_REF] Burrage | Accurate stationary densities with partitioned numerical methods for stochastic differential equations[END_REF] which provide a detailed analysis of the accuracy of the numerical invariant measure in the case of a linear drift, i.e. a quadratic potential V (q). It is shown there that among the class of stochastic Runge-Kutta methods applied to Langevin dynamics with a linear drift (which does not include the Lie-Trotter splitting (4)), the (stochastic) implicit midpoint rule is the only (stochastic) Runge-Kutta method with a nonsingular tableau matrix that reproduces the exact stationary distribution for all values of damping. See also related results in [START_REF] Li | Effectiveness of implicit methods for stiff stochastic differential equations[END_REF]. Our conclusion for Lie-Trotter splitting methods (4) is different: we prove in particular that it samples exactly the invariant measure for a linear drift if the deterministic integrator Φ h is a consistent symmetric Runge-Kutta method.

Our new characterization is expressed in terms of functions arising in the modified differential equation for the Hamiltonian system. We thus heavily rely on backward error analysis for our finding. The first error term in the error estimate for the numerical invariant measure can also be expressed in term of functions arising in the modified differential equation. This is also a new expression for error estimates à la Talay-Tubaro [START_REF] Talay | Expansion of the global error for numerical schemes solving stochastic differential equations[END_REF].

To simplify the presentation of our results, we focus on the globally Lipschitz case and we assume that the potential V is a C ∞ function where ∇V has bounded derivatives of any order (which implies the global Lipschitzness of ∇V ) and satisfies the following standard growth condition, which guaranties that (1) is ergodic,

q T ∇V (q) ≥ C 1 q T q -C 2 , for all q ∈ R d , (8)
for some constants C 1 , C 2 > 0. We further assume that the numerical flow Φ h is globally

Lipschitz in R 2d , Φ h (x 1 ) -Φ h (x 2 ) ≤ L x 1 -x 2 , (9)
for all x 1 , x 2 ∈ R 2d and all stepsize h small enough. In a more general non-globally Lipschitz setting, important in most applications, two implicit schemes are studied in [START_REF] Kopec | Weak backward error analysis for Langevin process[END_REF], and the implicitness is used to guarantee the boundedness of the numerical moments along time. Using implicit schemes for the Hamiltonian part is one way to extend our results in the non-globally Lipschitz setting. Alternatively, notice that the class of integrators considered in this paper (which are not assumed implicit) could still be applied rigorously in the non-globally Lipschitz setting by following the methodology of rejecting exploding trajectories proposed in [START_REF] Milstein | Numerical integration of stochastic differential equations with nonglobally Lipschitz coefficients[END_REF] and applied to Langevin-type ergodic systems in [START_REF] Milstein | Computing ergodic limits for Langevin equations[END_REF].

We close this introduction by mentioning the we focus here only on the discrepancy error [START_REF] Burrage | Accurate stationary densities with partitioned numerical methods for stochastic differential equations[END_REF], and do not discuss the Monte-Carlo error that can't be avoided in practice. One could use the so-called multi-level Monte-Carlo method (MLMC) [START_REF] Giles | Multilevel Monte Carlo path simulation[END_REF] which is a popular technique for reducing the variance. Indeed, it applies not only to first weak order integrators for SDEs but also to higher orders weak schemes, as shown in [START_REF] Giles | Antithetic Multilevel Monte Carlo estimation for multidimensional SDEs without Lévy area simulation[END_REF], where the antithetic MLMC has been introduced. Also in the context of stiff SDEs, in was shown in [START_REF] Abdulle | Improved stabilized Multilevel Monte Carlo method for stiff stochastic differential equations[END_REF] that applying a weak second order method at the finer level of the multilevel Monte-Carlo method permits to significantly improve the error constant. Finally we mention the schemes based on Markov Chain Monte-Carlo methods, see e.g. the survey [START_REF] Sanz-Serna | Markov Chain Monte Carlo and Numerical Differential Equation[END_REF], for which the bias error is reduced to zero, thanks to an appropriate acceptance/rejection criteria at each step.

The paper is organized as follows. In section 2, we present the basic ingredients of backward error analysis and the framework in [START_REF] Abdulle | High order numerical approximation of the invariant measure of ergodic SDEs[END_REF] that allows us to characterise the long time behaviour of numerical integrators of ergodic SDEs. In Section 3, we prove the main result of this paper, namely sufficient conditions for a Lie-Trotter splitting scheme to capture the numerical invariant measure of nonlinear ergodic Langevin dynamics up to an arbitrary order. We also relate this characterization to geometric properties of the deterministic integrator and present in Section 4 various numerical investigations, both for linear and nonlinear problems that corroborate such theoretical findings.

Preliminaries

In Section 2.1, we discuss some important concepts related to backward error analysis for ODEs, while in Section 2.2 we discuss the standard framework of backward Kolmogorov equation and weak Taylor expansions. Finally, in Section 2.3, we recall the order conditions framework from [START_REF] Abdulle | High order numerical approximation of the invariant measure of ergodic SDEs[END_REF] for the approximation of the invariant measure.

Deterministic backward error analysis

Backward error analysis is a powerful tool for the analysis of numerical integrators for differential equations [START_REF] Sanz-Serna | Numerical Hamiltonian problems[END_REF][START_REF] Leimkuhler | Simulating Hamiltonian Dynamics[END_REF][START_REF] Hairer | Geometric Numerical Integration. Structure-Preserving Algorithms for Ordinary Differential Equations[END_REF]. In particular, it is the main ingredient for the proof of the good energy conservation (without drift) of symplectic Runge-Kutta methods when applied to deterministic Hamiltonian systems over exponentially long time intervals. In our context it is useful to characterize the "generator" of the deterministic method Φ h of (5). Given a consistent integrator y n+1 = Φ h (y n ) for a system of ODEs

dy(t) dt = f (y(t)), (10)
the idea of backward error analysis is to search for a modified differential equation written as a formal series in powers of the stepsize h,

d y dt = f ( y) + hf 1 ( y) + h 2 f 2 ( y) + . . . , y(0) = y 0 (11)
such that (formally) y n = y(t n ), where t n = nh (in the above differential equation, we omit the time variable for brevity). This means that the numerical solution can be interpreted as the exact solution of a modified ODE. The vector fields f j in the series can be constructed for all reasonable integrators by induction on j [START_REF] Leimkuhler | Simulating Hamiltonian Dynamics[END_REF][START_REF] Hairer | Geometric Numerical Integration. Structure-Preserving Algorithms for Ordinary Differential Equations[END_REF], setting f 0 = f . Notice however that this series diverges in general for nonlinear systems and needs to be truncated. Considering the truncated modified ODE at order s

d y dt = f ( y) + hf 1 ( y) + h 2 f 2 ( y) + . . . + h s f s ( y), y(0) = y 0 (12)
we have y n = y(t n ) + O(h s+1 ) for h → 0 for bounded times t n = nh ≤ T . We note that the flow Φ h (y) of the modified differential equation ( 12) can be expressed as

Φ h = M k=0 h k L k D k! I+O(h M +1 ), L D = F 0 + hF 1 + h 2 F 2 + . . . + h s F s , (13)
for all M ≥ 0, where I is the identity map, F j φ = f j • ∇φ, j = 0, . . . , s and f 0 = f . This is a particular case (φ(y) = Iy) of the general result [13, Sect. III.5.1]

φ • Φ h = M k=0 h k L k D k! φ+O(h M +1 ), (14)
for all smooth test function φ. The constants symbolized by O in ( 13), ( 14), are independent of h → 0 but depend on M, s and φ.2

Order conditions for the numerical invariant measure

A numerical method applied to (1) is called ergodic if it has a unique invariant law µ h with finite moments of any order and lim

N →∞ 1 N N n=1 φ(p n , q n ) = R d ×R d φ(p, q)dµ h (p, q).
The aim of this subsection is to briefly describe the conditions that the numerical method applied to (1) should satisfy in order to approximate the invariant measure µ with an error (7) of size O(h r ). These conditions relate directly to the expansion of one-step numerical expectations in powers of h. In particular, let X 1 = (p 1 , q 1 ), x = (p, q) and

U (x, h) = E(φ(X 1 )|X 0 = x)
We assume that it is possible to expand U in powers of h to obtain the following expansion

U (x, h) = φ(x) + hA 0 φ(x) + h 2 A 1 φ(x) + h 3 A 2 φ(x) + • • • (15)
where A i , i = 0, 1, • • • are linear differential operators with coefficients depending on the choice of the numerical integrator and depending smoothly on the Hamiltonian function H(p, q) and its derivatives. For a consistent method (i.e. of weak order at least one), we have A 0 = L where L denotes the generator of (1) given by

L := p • ∇ q -(∇V (q) + γp) • ∇ p + β -1 γ∆ p . ( 16
)
where ∇ p , ∇ q denote the gradient differential operators with respect to p and q, and ∆ p is the Laplacian with respect to the p variable. Truncating the expansion (15) yields rigorous estimates thanks to the following result (see for instance [START_REF] Talay | Discrétisation d'une équation différentielle stochastique et calcul approché d'espérances de fonctionnelles de la solution[END_REF]).

Example 2.1. Consider the Euler-Maruyama method for (1) defined as

q n+1 = q n + hp n , p n+1 = p n -h(∇V (q n ) + γp n ) + 2β -1 γhξ n ,
where ξ n ∼ N (0, I) are independent d-dimensional Gaussian random variables. A straightforward calculation yields that A 1 is a differential operator of order 4,

A 1 φ = 1 2 φ ′′ (f, f ) + β -1 γ(∆ p φ) ′ f + β -2 γ 2 2 ∆ 2 p φ with f = -∇V (q) -γp p ,
where φ ′′ (•, •) denotes the second differential of φ (bilinear form), and φ ′ is the first derivative (linear form). We refer to [START_REF] Abdulle | High order weak methods for stochastic differential equations based on modified equations[END_REF][START_REF] Zygalakis | On the existence and the applications of modified equations for stochastic differential equations[END_REF] for other examples of such weak Taylor series calculations.

We can now state the conditions which ensure that µ h approximates µ to order r. We have the following theorem which involves the adjoint operators A * j of the differential operators A j involved in [START_REF] Kopec | Weak backward error analysis for overdamped Langevin processes[END_REF], defined as

R 2d A * j φ 1 (x)φ 2 (x)dx = R 2d φ 1 (x)A j φ 2 (x)dx
for all smooth test functions φ 1 , φ 2 with compact support on R 2d . The following result is proved in [START_REF] Abdulle | High order numerical approximation of the invariant measure of ergodic SDEs[END_REF].

Theorem 2.2. For the ergodic system (1), consider a consistent ergodic numerical method satisfying [START_REF] Kopec | Weak backward error analysis for overdamped Langevin processes[END_REF] and

A * j ρ ∞ = 0, for j = 1, • • • , r -1. ( 17
)
Then one obtains that the error in [START_REF] Burrage | Accurate stationary densities with partitioned numerical methods for stochastic differential equations[END_REF] satisfies for all test functions φ ∈ C

∞ P (R 2d , R), e(φ, h) = h r ∞ 0 R d ×R d A r u(p, q, t)ρ ∞ (p, q)dpdqdt + O(h r+1 ) (18)
where ρ ∞ is the invariant measure density given in (3), h is assumed small enough, and u(x, t) solves the backward Kolmogorov equation

∂u ∂t (p, q, t) = Lu(p, q, t), t > 0, p, q ∈ R d , (19a)
u(p, q, 0) = φ(p, q), ( 19b
)
and L is the generator of (1) defined in [START_REF] Leimkuhler | Rational construction of stochastic numerical methods for molecular sampling[END_REF].

Theorem 2.2 generalizes a standard result by Talay [START_REF] Talay | Expansion of the global error for numerical schemes solving stochastic differential equations[END_REF] and Milstein [START_REF] Milstein | Weak approximation of solutions of systems of stochastic differential equations[END_REF] (see e.g. [23, Chap. 2.2, 2.3] for a proof) for methods which have weak order r. Here, the order r assumption is replaced by the weaker condition [START_REF] Leimkuhler | The computation of averages from equilibrium and nonequilibrium Langevin molecular dynamics[END_REF]). An interpretation of ( 17) is the following: the invariant measure µ with density ρ ∞ of ( 1) is preserved up to a O(h r ) error by a single time step X 0 → X 1 of the integrator. Precisely, if X 0 is a random variable with law µ, then for all test functions φ ∈ C ∞ P (R 2d , R) and h → 0,

E(φ(X 1 )) = R 2d φ(x)ρ ∞ (x)dx + O(h r ),
where the constant symbolized by O is independent of h assumed small enough.

Remark 2.3. Assuming that ∇V is smooth with bounded derivatives at all orders, (19) has a unique solution for all test function φ ∈ C ∞ P (R 2d , R). It admits a rigorous Taylor expansion up to arbitrary order p for all h small enough

u(x, h) = φ(x) + hLφ(x) + h 2 2 L 2 φ(x) + . . . + h p p! L p φ(x) + O(h p+1 ),
where the function symbolized in O is smooth with respect to x with its all derivatives having a polynomial growth, independently of the smallness of the stepsize h. Notice that applying this result to the deterministic problem [START_REF] Hairer | Energy-preserving variant of collocation methods[END_REF], we recover the statement [START_REF] Kopec | Weak backward error analysis for Langevin process[END_REF].

3 Accuracy of Lie-Trotter splitting methods in sampling the invariant measure

In this section, we derive our main results. We first present in Section 3.1 sufficient order conditions for the Lie-Trotter splitting to achieve order r for the invariant measure. We discuss the relation of our new characterization with geometric properties such as energy and volume conservation and show in Section 3.2 that high order for the invariant measure can be achieved without high order conservation of energy and volume.

3.1 Order conditions to capture the invariant measure to a given order

Using arguments from backward error analysis, we prove the following result for the accuracy of the sampling the invariant measure (3) of Langevin dynamics of the form (1).

Theorem 3.1. Consider the ergodic system (1). Assume that Φ h is a consistent method for (5) satisfying (9) with modified equation from backward error analysis [START_REF] Giles | Antithetic Multilevel Monte Carlo estimation for multidimensional SDEs without Lévy area simulation[END_REF]. If the vector fields in [START_REF] Giles | Antithetic Multilevel Monte Carlo estimation for multidimensional SDEs without Lévy area simulation[END_REF] satisfy

div(f j ρ ∞ ) = 0, j = 1, . . . , r -1,
then, assuming ergodicity, the Lie-Trotter splitting (4) has order r of accuracy for the invariant measure. Precisely, we have for all φ ∈ C ∞ P (R d , R) and h → 0,

e(φ, h) = -h r ∞ 0 R d ×R d u(p, q, t) div(f r (p, q)ρ ∞ (p, q))dpdqdt + O(h r+1 ), ( 20
)
where e(φ, h) is defined in [START_REF] Burrage | Accurate stationary densities with partitioned numerical methods for stochastic differential equations[END_REF] and u(x, t) is defined in [START_REF] Li | Effectiveness of implicit methods for stiff stochastic differential equations[END_REF]. The constant symbolized in O is independent of the stepsize h assumed small enough.

Proof. Consider the method ( 4) and its corresponding one step expansion in powers of h

E(φ(X 1 )|X 0 = x) = φ(x) + hLφ(x) + • • • + h r A r-1 φ(x) + h r+1 A r φ(x) + O(h r+1 ).
In order to prove our theorem, it is enough to show that

A * j ρ ∞ = 0 for j = 1, • • • r -1, A * r ρ ∞ = div(f r ρ ∞ ). (21)
The result then follows immediately from Theorem 2.2 using the identity

R 2d A r u(x, t)ρ ∞ (x)dx = - R 2d u(x, t) div(f r (x)ρ ∞ (x))dx. (22)
We now start with the calculation of A j . In particular, given a smooth test function φ and x ∈ R 2d , using the semi-group property of the Markov process we have

E(φ(X 1 )|X 0 = x) = E(φ(Φ h • Θ h,n )(X 0 )|X 0 = x) = e hL S (φ • Φ h )(x), (23)
where e hL S φ denotes the exact flow of the Kolmogorov backward equation corresponding to [START_REF] Burrage | Numerical methods for second-order stochastic differential equations[END_REF], with generator L S given by

L S := -γp • ∇ p + β -1 γ∆ p .
We next recall the generator (13) of the truncated modified equation ( 12) of the integrator

Φ h L D = F 0 + hF 1 + . . . + h r F r ,
where we define the differential operators F j φ = f j • ∇φ (with f 0 = f ). We have, using [START_REF] Milstein | Stochastic numerics for mathematical physics[END_REF], applying Remark 2.3 for L S , and using ( 14) with M = s = r,

E(φ(X 1 )|X 0 = x) = r k=0 h k L k S k! r k=0 h k L k D k! φ(x) + O(h r+1 ) = φ(x) + hLφ(x) + r k=1 h k+1 A k φ(x) + O(h r+1 )
where we denote

A k = k+1 j=0 L k+1-j S 1 ≤ i ≤ j n 1 + n 2 + • • • + n i = j -i 1 i!(k + 1 -j)! F n 1 • • • F n i ,
where the second sum above is over integers n 1 , . . . , n i ≥ 0 and is equal to the identity I when j = 0. We obtain for all k ≥ 1,

A * k ρ ∞ = k+1 j=0 1 ≤ i ≤ j n 1 + n 2 + • • • + n i = j -i 1 i!(k + 1 -j)! F * n i • • • F * n 1 (L * S ) k+1-j ρ ∞ . Using L * S ρ ∞ = 0 and F * i ρ ∞ = 0, i = 1, • • • r-1
we see that for k ≤ r, the only possibly nonzero term in the above sum is obtained for j = r + 1, k = r, i = 1, i.e., F * r ρ ∞ = div(f r ρ ∞ ). Hence, we deduce (21) which permits to conclude the proof. Remark 3.2. Theorem 3.1 remains valid for the analogous scheme

X n+1 = Θ h,n •Φ h (X n )
a similar analysis.

An immediate consequence of Theorem 3.1 is the following corollary using the identity

div(f j ρ ∞ ) = div(f j ) -βf j • ∇H ρ ∞ . Corollary 3.3. Under the assumptions of Theorem 3.1, if div(f j ) -βf j • ∇H = 0, j = 1, . . . , r -1, (24)
then, assuming ergodicity, the Lie-Trotter splitting (4) has order r of accuracy for the invariant measure.

In addition to the estimates on the bias [START_REF] Burrage | Accurate stationary densities with partitioned numerical methods for stochastic differential equations[END_REF], as discussed in [START_REF] Abdulle | High order numerical approximation of the invariant measure of ergodic SDEs[END_REF] for Theorem 2.2, we also obtain the following exponential convergence estimate [START_REF] Milstein | Computing ergodic limits for Langevin equations[END_REF] which is classical for weak methods [23, Chap. 2] (see also [START_REF] Talay | Stochastic Hamiltonian systems: exponential convergence to the invariant measure, and discretization by the implicit Euler scheme[END_REF][START_REF] Mattingly | Ergodicity for SDEs and approximations: locally Lipschitz vector fields and degenerate noise[END_REF][START_REF] Milstein | Computing ergodic limits for Langevin equations[END_REF] in the context of non-globally Lipchitz vector fields). Remark 3.4. Under the assumptions of Theorem 2.2 or Theorem 3.1 or Corollary 3.3, exists constants λ, C, K(x) > 0 such that for all n ≥ 0,

E(φ(X n )) - R d ×R d φ(p, q)ρ(p, q)dpdq ≤ K(x)e -λtn + Ch p ( 25
)
where t n = nh, the constants C, K(x), λ are independent of n, h, and h is assumed small enough (C, K(x) depend on φ and K(x) depends on X 0 = x).

The simplest way for condition [START_REF] Milstein | Numerical integration of stochastic differential equations with nonglobally Lipschitz coefficients[END_REF] to be satisfied is to consider a method Φ h of order r, since in this case f j = 0, j = 1, • • • r -1. Indeed, in this case, standard results from backward error analysis show that the right-hand side of the modified equation ( 12) reads f ( y) + h r f r (y) + . . . and the f j = 0, j = 1 . . . , r -1.

Relation with results of [START_REF] Bou-Rabee | Long-run accuracy of variational integrators in the stochastic context[END_REF] In the paper [START_REF] Bou-Rabee | Long-run accuracy of variational integrators in the stochastic context[END_REF] the authors show that if a variational (equivalently symplectic) integrator Φ h preserving the energy with order r is used for solving the deterministic dynamics (5), then the Lie-Trotter splitting (4) has order r of accuracy for the invariant measure. Now a symplectic method applied to a Hamiltonian system is volume preserving hence div(f j ) = 0 for all j ≥ 0. Furthermore, if the method Φ h has order r, then it also preserves energy up to order r, i.e. for all vectors y,

H(Φ h (y)) = H(y) + O(h r+1 ) ( 26
)
this condition is equivalent to f j • ∇H = 0 for all 1 ≤ j r -1. If both terms in [START_REF] Milstein | Numerical integration of stochastic differential equations with nonglobally Lipschitz coefficients[END_REF] are required to vanish independently, we see that we obtain order r of accuracy for the invariant measure under weaker assumptions than in [START_REF] Bou-Rabee | Long-run accuracy of variational integrators in the stochastic context[END_REF]. It is not necessary for the method to be symplectic, only volume preservation up to to order r is required i.e., det ∂Φ h (y)

∂y = 1 + O(h r+1 ) ( 27
)
which is equivalent to div(f j ) = 0 for all 1 ≤ j ≤ r-1. Note that any deterministic method Φ h of order r automatically satisfies ( 26) and ( 27), and thus has the desired properties. We will further see in Section 3.2 that energy and volume conservation up to order r of the deterministic integrator Φ h in ( 4) are not necessary to get order r of accuracy for the invariant measure, as one can construct deterministic methods that satisfy [START_REF] Milstein | Numerical integration of stochastic differential equations with nonglobally Lipschitz coefficients[END_REF] with j = r + 1, while conserving energy and volume only up to order r.

Remark 3.5. In general, it is difficult to construct a numerical scheme which is simultaneously energy and volume preserving, so a choice has to be made. Indeed, already for problems with only one degree of freedom (d = 1), energy preserving means that the numerical method has an exact trajectory {y ∈ R 2d ; H(y) = H(y 0 )}, and thus is a time transformation of the exact solution. Such scheme can be also volume preserving if and only if this time transformation has divergence zero, which is not satisfied by known energy preserving methods such as the AVF method [START_REF] Quispel | A new class of energy-preserving numerical integration methods[END_REF] or energy preserving variants of collocation methods [START_REF] Hairer | Energy-preserving variant of collocation methods[END_REF].

Order of accuracy in the linear case The following corollary of Theorem 3.1 states that for (non-partitioned) Runge-Kutta methods and quadratic potentials, the order of accuracy for the invariant measure of the Lie-Trotter splitting (4) is always odd, i.e. if the Runge-Kutta method has order r, the order of accuracy for the invariant measure of the Lie-Trotter splitting is r + 1 if r is even, and r if r is odd. Moreover, if the Runge-Kutta method is time-symmetric, i.e. Φ -1 h = Φ -h then the Lie-Trotter splitting (4) samples exactly the invariant measure. We recall [START_REF] Hairer | Geometric Numerical Integration. Structure-Preserving Algorithms for Ordinary Differential Equations[END_REF] that the method Φ h is symmetric if and only if [START_REF] Giles | Antithetic Multilevel Monte Carlo estimation for multidimensional SDEs without Lévy area simulation[END_REF] has an expansion in even powers of h, i.e. f 2j+1 = 0 for all j.

(28)

Corollary 3.6. In the linear case (i.e. V (q) is quadratic), assume that Φ h is a consistent (non-partitioned) Runge-Kutta method for [START_REF] Bou-Rabee | Pathwise accuracy and ergodicity of metropolized integrators for SDEs[END_REF]. Then the order of accuracy r for sampling the invariant measure of the Lie-Trotter splitting (4) is odd. Moreover, if Φ h is symmetric, then the error (7) of ( 4) is zero for all h small enough.

Proof. Since the problem is linear, the Hamiltonian vector field in (5) has the form f (y) = Ay = J -1 Sy where S is a constant symmetric matrix and

J = 0 I -I 0 ( 29
)
is the 2d × 2d symplectic matrix. For a Runge-Kutta method applied to a linear problem, then we have the identity Φ(x) = R(hA)x where R(z) is a rational function. The modified differential equation for backward error analysis of Φ h then has the form (see [13, Sect.

d y dt = h -1 log(R(hA)) y = (A + α 1 hA 2 + α 2 h 2 A 3 + α 3 h 3 A 4 + . . .) y, (30)
with constants α 1 , α 2 , α 3 , . . .. Using (28), we have α j = 0 for all odd j in the case of a symmetric method. Since f j (y) = α j A j+1 , to conclude the proof using Corollary 3.3, it is sufficient to show that div(f j ) = f j • ∇H = 0 for all even integer j. Indeed, for Φ h of order r, we already have f j = 0 for all 1 ≤ j < r, while for symmetric Φ h , we already have f j = 0 for all odd j. Since j is even, we have that JA j+1 is a symmetric matrix and thus f j is a Hamiltonian vector field (f j = J -1 ∇H j with H j (y) = 1 2 y T JA j+1 y), which yields div(f j ) = 0. In addition, using the skew-symmetry of SA j+1 , we obtain f j (y) • ∇H(y) = α j y T SA j+1 y = 0. This concludes the proof.

Energy and volume conservation are sufficient, but not necessary conditions

Corollary 3.3 shows that energy and volume conservation up to order r of the deterministic integrator Φ h in (4) are sufficient conditions for the splitting method (4) to have order r for the invariant measure. These are however not necessary, as shown in the following theorem which proves that applying an appropriate time transformation to a symplectic integrator of order r permits to increase to order r + 1 the accuracy for the invariant measure, while the order of accuracy for the conservation of both the energy (26) and volume (27) remains r. We emphasize that this time transformed numerical method is only used to illustrate a theoretical finding and we do not advocate its use in practice.

Here, we use the fact that for a symplectic (possibly partitioned) Runge-Kutta method Φ h applied to a (deterministic) Hamiltonian system ẏ = f (y) = J -1 ∇H(y), the modified equation for backward error analysis (11) also possesses an Hamiltonian structure of the form [13, Sect. IX.3]

d y dt = J -1 ∇H h ( y), H h = H + h r H r + h r+1 H r+1 + . . . , (31)
where J in defined in [START_REF] Sanz-Serna | Markov Chain Monte Carlo and Numerical Differential Equation[END_REF], and the formal series H h is the modified Hamiltonian.

Theorem 3.7. Assume that Φ h is a consistent symplectic integrator of order r for (5) with modified Hamiltonian [START_REF] Shardlow | A perturbation theory for ergodic markov chains and application to numerical approximations[END_REF]. Consider the modified integrator Φ h defined via a time transformation by

Φ h (y) := Φ α h (y)h (y), α h (y) := 1 -h r βH r (y).
Then, assuming ergodicity, the splitting method

X n+1 = Φ h • Θ h,n (X n ) (32)
has order r + 1 for the invariant measure of (1).

Proof. We denote f h = f + h r f r (f r = J -1 ∇H r ) and f h = f + h r f r the vector fields [START_REF] Hairer | Energy-preserving variant of collocation methods[END_REF] of backward error analysis of Φ h and Φ h , repectively, truncated at order r. Using

α h = 1 + O(h r ), we have f h = α h f h + O(h r+1 ). We deduce div(f h ) -βf h • = α h div(f h ) + ∇α h • f h -βα h f h • ∇H + O(h r+1 ) = α h div(f h ) -βα h (∇H + h r ∇H r ) • f h + O(h r+1 ) = O(h r+1 )
where we use div(f h ) = 0 and an easy computation shows that (∇H + h r ∇H r ) • f h = 0. This yields div(f r )-βf r •∇H = 0. This shows that Φ h fulfils the assumptions of Corollary 3.3 with r replaced by r + 1 in [START_REF] Milstein | Numerical integration of stochastic differential equations with nonglobally Lipschitz coefficients[END_REF], which concludes the proof.

It is known for deterministic geometric integrators [13, Sect. VIII.2] that a non constant time transformation destroys in general the geometric properties of symplecticity and volume conservation. Thus, it is not surprising that the scheme ( 32) is energy and volume preserving up to order r but not up to order r + 1. Indeed, we have that div

(f h ) = div(α h f h ) + O(h p+1 ) = ∇α h • f h + O(h r+1 ) = -h r βα h ∇H r • f h + O(h r+1 ) is O(h r ) in general because ∇H r • f h = 0 (
otherwise, Φ h would be energy preserving up to order r + 1). This implies that the idea of Theorem 3.7 cannot be applied repeatedly to achieve arbitrarily high order.

Remark 3.8. The time transformation α h (y) in Theorem 3.7 is non-globally Lipschitz in general (see e.g. the formula (36) for the symplectic Euler method). This may lead to a non-globally Lipschitz numerical flow. One possibility to make α h (y) Lipschitz and bounded is to consider a modification, e.g., α h = 1 -(1 + h|H r | 2 ) -1 βh r H r ), which does not affect f, fr in the modified equation of Φ h .

Numerical experiments

In this section we illustrate the various theoretical results of this paper. We start in Section 4.1 by an example with a quadratic potential, that illustrates the results of Corollary 3.6. We then illustrate in Section 4.2 the finding of Theorem 3.1 by an example with quartic potential with various deterministic integrators Φ h . This section also corroborates that symplecticity is not a necessary geometric property to capture the invariant measure of Langevin dynamics with high order.

Illustration for quadratic potential

In Figure 1, we consider the Langevin dynamics (1) in R 2 with one degree of freedom (d = 1) with quadratic potential V (q) = 1 2 q 2 , which reduces (5) to the harmonic oscillator. We observe that (4) has the form X n+1 = U h X n +V h ξ n where U h , V h are constant matrices. This yields that X n is a Gaussian variable for all n, and a calculation shows that the covariance matrix Σ h of the corresponding numerical invariant measure is given by the Lyapunov equation

U h Σ h U T h -Σ h = V h V T h .
We compare in Figure 1 the covariance matrix Σ h with the covariance matrix Σ = β -1 I of the exact Gaussian invariant measure and plot the 2-norm of the difference. As predicted 1) in the linear case (V (q) = 1 2 q 2 ). The lines corresponds to explicit deterministic integrators of orders p = 1, 2, 3, . . . , 9 (from top to bottom), respectively. The orders of accuracy for the invariant measure are always odd. by Corollary 3.6, we observe the expected lines of with odd slopes 1, 3, 3, 5, 5, 7, 7, 9, 9 for the explicit deterministic integrators Φ h (y) = (I + A + A 2 2 + . . . + A p p! )y (with A = J -1 ) of order p, where p = 1, 2, 3, . . . , 9, respectively.

Illustration for quartic potential

Although our analysis applies only to globally Lipschitz vector fields, in Figure 2, we consider the Langevin dynamics (1) with one degree of freedom with quartic potential V (q) = (1q 2 ) 2 -1 2 q. We compare the invariant measure error (7) for φ(p, q) = p 2 + q 2 of the Lie-Trotter splitting method (4) versus the stepsize h for various choices of the deterministic integrator Φ h :

• the explicit Euler method,

p n+1 = p n -h∇V (q n ), q n+1 = q n + hp n , (33)
• the Heun method,

p n+1 = p n -h∇V q n + h 2 p n , q n+1 = q n + h p n - h 2 ∇V (q n ) , (34)
• the standard symplectic Euler method,

p n+1 = p n -h∇V (q n ), q n+1 = q n + hp n+1 , (35)
• A time transformed version of the symplectic Euler method, 7) for φ(p, q) = p 2 + q 2 of the Lie-Trotter splitting method (4) versus the stepsize h for Langevin dynamics (1) with quartic potential V (q) = (1q 2 ) 2 -1 2 q (d = 1), with γ = 4, β = 2 and initial condition p(0) = q(0) = -1.5. We compare various choices of the deterministic integrator Φ h : explicit Euler (33), symplectic Euler (36), Heun explicit RK2 [START_REF] Talay | Expansion of the global error for numerical schemes solving stochastic differential equations[END_REF], and the modified symplectic Euler method (36). The curves are obtained as the averages over 10 independent long trajectories of time length T = 10 8 . The vertical bars measure the standard deviation intervals due to the Monte-Carlo error among these 10 trajectories.

p n+1 = p n -α n h∇V (q n ), q n+1 = q n + α n hp n+1 , α n = 1 + h 2 βp n • ∇V (q n ). ( 36
The explicit Euler method, the standard symplectic Euler method and its time transformed version have deterministic order 1 for (5), while the Heun method has deterministic order 2. In Figure 2, we plot the error (7) for φ(p, q) = p 2 + q 2 for the methods considered above. We consider for each method the average over 10 independent long time trajectories with length T = 10 8 . The vertical bars represent the standard deviation intervals due to the Monte-Carlo error among these 10 trajectories. As predicted by Corollary 3.3, the expected convergence slopes 1 and 2 can be observed in Figure 2. Precisely, the method (36) has deterministic order 1 and preserves the energy and the volume only up to order r = 1 in ( 26), [START_REF] Quispel | A new class of energy-preserving numerical integration methods[END_REF], but it has order 2 of accuracy for the invariant measure, as predicted by Theorem 3.7. We observe in Figure 2 that the order 2 schemes have an accuracy improved by a factor of about 100 compared to the order 1 methods. Again, as mentioned in Section 3.2, the time transformed symplectic Euler method (1) is only presented to illustrate Theorem 3.7 and we do not claim that it is the best method at hand to solve problem [START_REF] Abdulle | Improved stabilized Multilevel Monte Carlo method for stiff stochastic differential equations[END_REF].

Geometric properties of symplectic and time transformed symplectic methods.

In this experiment we take again the quartic potential V (q) = (1q 2 ) 2 -1 2 q, temperature β = 2 and initial condition p(0) = q(0) = -1.5. In Figure 3, we plot the evolution of the Hamiltonian energy (2) along time for the standard symplectic Euler method [START_REF] Zygalakis | On the existence and the applications of modified equations for stochastic differential equations[END_REF] and the time transformation version (36) as described in Theorem 3.7 for two values of γ.

In the top picture of Figure 3, we consider the deterministic case (γ = 0). In contrast to the standard symplectic Euler method for which no drift in the energy can be observed, we obtain a drift in the energy for the time transformed version (36) which is non symplectic. In the bottom picture, we consider the stochastic case (γ = 4) and plot sample trajectories for both methods using the same sets of generated random variables. We observe no linear [START_REF] Zygalakis | On the existence and the applications of modified equations for stochastic differential equations[END_REF] and the time transformed version (36) (non symplectic) applied to (1) with V (q) = (1q 2 ) 2 -1 2 q (quartic potential). Top picture: deterministic case (γ = 0). Bottom picture: stochastic case (γ = 4). Stepsize is h = 0.02 and β = 2. drift in the energy for both methods in this stochastic case. It can be observed in Figure 2 that the modified symplectic Euler method yields order 2 of accuracy compared to the standard version of order 1, as predicted by Theorem 3.7. This last numerical example shows that among the standard symplectic Euler method [START_REF] Zygalakis | On the existence and the applications of modified equations for stochastic differential equations[END_REF] and the time transformed version (36) (applied in the Lie-Trotter splitting (4)), the most accurate scheme is not the same in the deterministic and stochastic contexts, as predicted by Theorem 3.7.

Conclusion

We have presented an analysis supported by numerical experiments of Lie-Trotter splitting methods applied to nonlinear Langevin dynamics. In particular, we obtained a new characterization of sufficient conditions for the Lie-Trotter splitting to capture the numerical invariant measure of nonlinear ergodic Langevin dynamics up to arbitrary order. While in previous results in the literature symplecticity and energy conservation (up to a given order) were required to prove high order accuracy for the numerical invariant measure, we showed that although sufficient, these conditions are not necessary. More precisely, we showed that the order of convergence of the deterministic integrator provides a lower bound for the order of convergence for the numerical invariant measure of the Lie-trotter splitting method. We illustrated this finding by constructing a first order non symplectic deterministic integrator that, despite its poor geometric properties, captures the invariant measure of Langevin dynamics to second order.

9 Figure 1 :

 91 Figure 1: Accuracy of the numerical invariant measure (covariance matrix error) of the Lie-Trotter splitting for Langevin dynamics (1) in the linear case (V (q) = 1 2 q 2 ). The lines corresponds to explicit deterministic integrators of orders p = 1, 2, 3, . . . , 9 (from top to bottom), respectively. The orders of accuracy for the invariant measure are always odd.

2 ref ere nce slo pe 1 Figure 2 :

 212 Figure2: Invariant measure error[START_REF] Burrage | Accurate stationary densities with partitioned numerical methods for stochastic differential equations[END_REF] for φ(p, q) = p 2 + q 2 of the Lie-Trotter splitting method (4) versus the stepsize h for Langevin dynamics (1) with quartic potential V (q) = (1q 2 ) 2 -1 2 q (d = 1), with γ = 4, β = 2 and initial condition p(0) = q(0) = -1.5. We compare various choices of the deterministic integrator Φ h : explicit Euler (33), symplectic Euler (36), Heun explicit RK2[START_REF] Talay | Expansion of the global error for numerical schemes solving stochastic differential equations[END_REF], and the modified symplectic Euler method (36). The curves are obtained as the averages over 10 independent long trajectories of time length T = 10 8 . The vertical bars measure the standard deviation intervals due to the Monte-Carlo error among these 10 trajectories.

Figure 3 :

 3 Figure 3: Evolution of the Hamiltonian energy (2) along time for the standard symplectic Euler method[START_REF] Zygalakis | On the existence and the applications of modified equations for stochastic differential equations[END_REF] and the time transformed version (36) (non symplectic) applied to (1) with V (q) = (1q 2 ) 2 -1 2 q (quartic potential). Top picture: deterministic case (γ = 0). Bottom picture: stochastic case (γ = 4). Stepsize is h = 0.02 and β = 2.

Mathematics Section, École Polytechnique Fédérale de Lausanne, Station 8, 1015 Lausanne, Switzerland, Assyr.Abdulle@epfl.ch

2 Université de Genève, Section de mathématiques, 2-4 rue du Lièvre, CP 64, 1211 Genève 4, Switzerland, Gilles.Vilmart@unige.

ch[START_REF] Abdulle | High order numerical approximation of the invariant measure of ergodic SDEs[END_REF] Mathematical Sciences, University of Southampton, Southampton SO17 1BJ, UK, k.zygalakis@soton.ac.uk

This means that the exact solution (p(h), q(h)) and the numerical solution Θ h,n (p(0), q(0)) have the same law of probability for a given state (p(0), q(0)).

For all h small enough, the sum in (13) (and (14)) can be shown to converge for M → ∞ in the case of analytic vector fields fj (and analytic test functions φ), which permits to remove the O remainder.

Acknowledgements. The work of AA and GV was partially supported by the Fonds National Suisse, project No. 200021 140692. and No. 200020 144313/1, respectively.


Mauro Napoletano
email: mauro.napoletano@sciencespo.fr

Jean-Luc Gaffard
email: jeanluc.gaffard@sciencespo.fr

Andrea Roventini
email: andrea.roventini@univr.it

Outside the corridor: fiscal multipliers and business cycles into an agent-based model with liquidity constraints

Keywords: Keynesian economics, Fiscal Multipliers, Corridor Effects, Agent-based models, Liquidity constraints JEL classification: E63, E21, C63

We build an agent-based model to study how fiscal multipliers can change over the business cycle. Our approach considers the economy as a complex evolving system. In that, fiscal state-dependent multipliers are emergent disequilibrium phenomenon stemming from the interaction among an ecology of heterogeneous agents. We study fiscal multipliers in response to different microeconomic shocks hitting the economy. We show that deficit-spending fiscal policy dampens the effect of a shock and lowers its persistence. Moreover, we show that the size and dynamics of the fiscal multiplier is inversely related to the evolution of credit rationing in the aftermath of the shock. We also investigate the effects of two different balanced budget rules. In the first type of such experiments, government expenditure is constrained to be equal to tax revenues of each period. In the second one the tax rate is eventually raised to balance a given level of government expenditure. We show that fiscal multipliers are very low with both balanced-budget rules. Finally, we show that fiscal multipliers are higher into more leveraged economies.

Introduction

In the economics of Keynes, agents who are unable to sell as much as they would like at prevailing prices restrict demands in other markets. Unemployed workers cut their consumption demand. Demand-constrained firms restrict their demand for labor. Then supplies and demands depend on the current income and and multiplier effects matter.

The recent crisis has revealed that multiplier values can be significantly higher in economic downturns than in normal times. For instance, in 2009 the IMF had estimated that multipliers in developed countries averaged about 0.5, while it now calculates that they range from 0.9 to 1.7 during the Great Recession (see [START_REF] Blanchard | Growth forecast errors and fiscal multipliers[END_REF]. The notion of state-dependent multipliers, corroborated by a growing body of recent empirical research (see e.g. [START_REF] Auerbach | Measuring the Output Responses to Fiscal Policy[END_REF][START_REF] Bachmann | Confidence and the transmission of government spending shocks[END_REF], is reminiscent of [START_REF] Leijonhufvud | Effective demand failures[END_REF] corridor theory. According to Leijonhufvud, "inside the corridor, multiplier-repercussions are weak and dominated by neoclassical market adjustments, while outside the corridor, they should be strong enough for effects of shocks to the prevailing state to be endogenously amplified" [START_REF] Leijonhufvud | Effective demand failures[END_REF]. One explanation for the corridor effect is the tightness of the liquidity constraints faced by households. During booms liquidity constraints are non-binding, and households engage into standard permanent-income consumption smoothing. Consequently, aggregate consumption depend less on current income and the size of the fiscal multiplier is low. In contrast, during economic downturns liquidity constraints are binding, and consumption is more sensitive to current income variations. Recent DSGE models (e.g. [START_REF] Eggertsson | Debt, Deleveraging, and the Liquidity Trap: A Fisher-Minsky-Koo Approach*[END_REF] have modeled consumers debt overhangs as exogenous, or have investigated corridor effects in frameworks with limited agent heterogeneity. In our model agents have heterogeneous financial conditions and debt overhangs may emerge endogenously. Moreover, the size of fiscal multipliers affected by monetary policy as well as by structural conditions (e.g. distribution of household leverage). Our approach considers the economy as a complex evolving system. In that, fiscal state-dependent multipliers are emergent disequilibrium phenomenon stemming from the interaction among an ecology of heterogeneous agents. We study fiscal multipliers in response to different microeconomic shocks hitting the economy. We show that deficit-spending fiscal policy dampens the effect of a shock and lowers its persistence. Moreover, we show that the size and dynamics of the fiscal multiplier is inversely related to the evolution of credit rationing in the aftermath of the shock. We also investigate the effects of two different balanced budget rules. In the first type of such experiments, government expenditure is constrained to be equal to tax revenues of each period. We find, that fiscal multipliers are still time varying, but they are lower and converge towards a zero value. Similar results are observed under the second type of balanced budget rule, according to which the tax rate is eventually raised to balance a given level of government expenditure. The paper is organized as follows. Section 2 describes the model. Section 3 discusses the steady state conditions of the model. Section 4 presents some preliminary simulation results. Finally, Section 5 concludes.

A simple endowment economy

There are N heterogeneous households. Each household i has access to a given amount of a homogeneous input good ('wheat') that she can gather at no cost, but that she cannot use for consumption . In order to consume, each household has to sell her endowment good to homogeneous firms ('the mills') that use it to produce an output good ('flour'). The output good is purchased by households and used for their consumption Each mill produces the consumption good on demand by using a constant returns to scale technology. The production function of a mill j is Y jt = L jt . This implies that total output is simply

Y t = L t
where L t is the amount of wheat purchased. Mills always buys the wheat from each households in proportion to the household share of total wheat available. In addition, each firm can buy all the input good she needs at the price P l up to to the total amount of input available (L max ). Furthermore, the mill produces at zero profit so that the final price of the output is simply P o = P l . Finally, notice that the above assumptions imply that the level of output of the mill will (and thus household incomes) will be determined by the level of overall consumption demand (up to the level where the input constraint L max becomes binding). Each household has a desired level of consumption Z i , that we assume to be constant over time Once desired consumption is determined it is possible to identify two classes of agents: savers and borrowers. Savers are households whose current liquid wealth W it is larger than desired consumption (i.e. Z i ≤ W it ). At the other extreme, borrowers have liquid wealth W it that is lower than desired consumption Z it (i.e. Z i > W it ). Savers can always finance their consumption with their own wealth. Accordingly, consumption of this class of agents is always equal to their desired level. In contrast, borrowers need financing from the credit sector to satisfy their consumption plans.

In the credit sector there is a representative bank that stocks the wealth of all agents and grants credit to borrowers. Total credit supply is set as a multiple of the net worth of the bank E B t (e.g. Delli Gatti et al., 2005)

T S t = kE B t , k > 0 (1)
where k is the credit multiplier. We assume an endogenous money framework, so that → k > 1. Furthermore, credit is allocated to agents using a pecking order (see [START_REF] Dosi | Income distribution, credit and fiscal policies in an agent-based Keynesian model[END_REF] that depends on the on the ratio between household's wealth W it and her credit demand

CD it . W it CD it
where CD it is equal to the difference between the borrower's desired consumption and her wealth

CD it = Z it -W it
The above assumptions imply that total credit demand can be lower or higher than total credit supply. In the former case all borrowers are able to get credit from the bank. In the latter case some borrowers are partially or totally rationed in the credit market. The lower is the position of the household in the pecking order, the higher is her probability of being credit-rationed. Borrowers who are denied credit are not able to satisfy their consumption plan. Their actual consumption is equal to their wealth W it .

The bank sets the interest rate on loans r b by applying a mark-up µ b on the baseline interest rate r set by the central bank. Likewise, the interest rate paid on deposits, r s is determined by applying a mark-down µ s on the central bank interest rate. We get

r b = r(1 + µ b ) r s = r(1 + µ s ) (2)
Bank liabilities, L B t are determined as the difference between assets (equal to total credit supply) and net-worth of the bank

L B t = kE B t -E B t = (k -1)E B t
Bank profits, π B t are simply equal to

π B t = r b t (kE B t ) -r s (k -1)E B t = [r s + k(r b -r s )]E B t (3)
The law of motion of bank's net worth is

E B t = E B t-1 + π B t - N i=1 BD it
Bank profits are distributed to a homogeneous class of agents ("the bankers") if bad debt is zero and are not distributed otherwise. The bankers fully consume their income.

Let us now turn to describing fiscal policy. We assume that taxation is proportional to income at the tax rate t so that disposable income of household i is simply y D it = (1-τ )y it . The government sets consumption level and the tax rate according to different fiscal rules. Government debt (if any) is purchased by the central bank. Aggregate demand AD t determined as the sum of households and government consumption, respectively C t and G t , plus the consumption of bankers if any ( and equal to bank profits π B t ). As long as the constraint L max is not aggregate income is determined by aggregate demand. Formally, we get:

Y t = AD t = C t + G t + π B t Total households income Y H t is total income minus the income of bankers, i.e. Y H t = Y -π B t .
Finally, let us discuss households' balance sheet dynamics. Let define β it = Z i /W it as the propensity to consume out of wealth of agent i at time t. It follows that, β it > 1 if the household is a borrower and β it <= 1 if she is a saver. We assume that consumption loans must be fully repaid at the end of each period. The same occurs for the remuneration of savings. It follows that the law of motion of agent's wealth is

W it+1 = (1 -τ )y it -(1 + r b )(β it -1)W it (4)
if the agent is a borrower, and

W it+1 = (1 -τ )y it + (1 + r s )(1 -β it )W it (5)
if the agent is a saver. In the above equations r b and r s (with r b ≥ r s ) are respectively the borrowers' and lenders' interest rate. Households who are unable to repay their debt go bankrupt. From 4 one gets that the condition for bankruptcy is simply that household's resources at the beginning of the period are lower than debt plus interests, i.e.

(1

-τ )y it < (1 + r b )(β it -1)W it
The above can also be re-expressed using consumption levels

(1 -τ )y it < (1 + r b )(C it -W it )
Once bankruptcy occurs the wealth of the household is reset to zero and the bank gets a credit loss equal to

BD it = (1 + r b )(C it -W it ) -(1 -τ )y it
Bankrupted households are denied access to the credit market for T def ault periods. The sequence of events unfolding in each period is the following

• Desired consumption of the period is determined. Households credit demand is determined

• Government consumption and the government balance is determined To find the steady state conditions of the model described in the previous section we focus on steady state conditions for households wealth. Indeed, household wealth plays a central role in the model, as it determines actual consumption, aggregate income, tax revenues (and thus government balance). In addition, by determining the overall amount of consumption loans it affects banks' revenues and net worth. First, notice that the assumption that each mill buys wheat from a household in proportion to the household's share of total wheat implies that each household is entitled a time-invariant share of total household income. Let us label this share by α i . If Y H t is total household income of the period, then disposable income of each household; y it is equal to

y it = α i (1 -τ )Y H t
As we said in the previous section, households whose marginal propensity to consume out of wealth is lower or equal to 1 always satisfy their consumption plans. In contrast, agents with β i > 1 need to borrow to achieve their desired consumption. We focus on steady states of the model where credit rationing is absent. Accordingly, all borrowers are able to get their desired consumption level. Using this fact, we get that the steady state level of wealth household i is

w * i = α i (1 -τ )Y H * [1 -(1 + r b )(1 -β * i )] (6)
if the agent is a borrower and

w * i = α i (1 -τ )Y H * [1 -(1 + r s )(1 -β * i )] (7)
if the agent is a saver. notice that the above steady state levels are stable if the following conditions are satisfied

| (1 + r b )(1 -β * i ) |< 1 for a borrower | (1 + r s )(1 -β * i ) |<
* i β * i α i [1 -(1 + r b )(1 -β * i )] = γ * i (9)
for a borrower and

β * i α i [1 -(1 + r s )(1 -β * i )] = γ * i (10)
for a saver. Solving for the above equations for β * i we get

β * i = γ * i r b [γ * i r b + (γ * i -α i )] (11)
and

β * i = γ * i r s [γ * i r s + (γ * i -α i )] (12)
for a saver. 1 . Furthermore, by exploiting Equations 11 and 12 we designed an algorithm that, for a given distributions of income shares {α i }, randomly assigns consumption weights to households and computes the values of β * i so that the fraction of borrowers in the population is 0 < η * < 1 (and the one of savers is 1 -η * ).

1 Incidentally, notice that the above functions are increasing in α * i and decreasing in γ * i . It follows that, if we impose that all agents have the same consumption share in steady state (γ i = 1/N, ∀i) then borrowers should be concentrated in the upper part of the income distribution.

Simulation Experiments

We now discuss the results of the simulation experiments using different fiscal rules in the framework of the model described in Section 2. In these simulation experiments we assumed that the economy was initially in the steady state described in Section 3. We then shocked the steady state by assuming that a given fraction of households goes bankrupt and tracking the dynamics of output under different types of fiscal rules. In particular, we considered three types of fiscal rules:

• deficit-spending rule. The government keeps the level government spending at the steady state level, and it is allowed to create a deficit if tax revenues fall below the steady state level

• balanced-budget rule I. Government spending is equal to tax revenues in every period

• balanced-budget rule II. The government keeps the level government spending at the steady state level. The tax rate is adjusted accordingly so that tax revenues are equal to government spending in every period. In each of the above experiment the level of government expenditure (and the tax rate when required) is set as a fraction of the steady state income level, according to a coefficient that captures the "intensity" of fiscal policy in the model. We then repeat each experiment for different intensities of fiscal policy.

Let us begin with the analysis the deficit-spending rule. Figure 1 shows the evolution of aggregate income as a fraction of steady state income in the aftermath of a bankruptcy shocks. The plot shows two results of the model. First, higher levels of government expenditures (corresponding to higher fiscal policies intensities in the figure) reduce the magnitude of the effects of the shocks. In other words the higher the level of government expenditure the lower is the fall in aggregate output in the aftermath of a shock. Moreover, higher levels of government expenditure lower the persistence of the shock and favors a quicker recovery to the steady state level of income.

In light of the above results we now proceed to calculate fiscal multipliers corresponding to different levels of fiscal policy intensity. More precisely, we calculate fiscal impact multipliers as follows

k ki (t) = Y k (t) -Y i (t) G k (t) -G i (t) (13)
In the above expression k ki (t) is the fiscal multiplier at time t corresponding to the experiment with govt. expenditure intensity k in relation to the experiment govt. expenditure intensity i.

Figure 2 shows the dynamics of fiscal impact multipliers for different intensities of government expenditure and following a bankruptcy shock. First, the values of fiscal multipliers follow an inverse U-shaped pattern, and peak values are reached after the through of the recession (compare Fig. 2 to Figure 1). Moreover, many fiscal multipliers have peak values that are significantly higher than one. In particular, the value of the peak multiplier is decreasing in the level of the government expenditure. This latter result can be understood by noticing that the impact multipliers are calculated as a ratio between the variation in aggregate income and the variation in government expenditure. As aftershock fall in output is at most of 10% (check the blue line in Figure 1, corresponding to the case where fiscal policy is absent), and given that aggregate output is bounded above by the capacity constraint L max it comes to no surprise that government intensities values equal or higher than 10% are associated with peak multipliers that are lower than one.

Let us now turn to analyse the evolution of fiscal multipliers in relation to the degree of credit rationing in the economy. Figure 3 compares the evolution of fiscal multipliers to the one of the ratio between actual and desired consumption. A value of the latter lower than one indicates the presence of credit rationing because it implies that some borrowers are unable to satisfy their desired consumption plans and their consumption is limited by their current liquid wealth. Figure 3 shows that the evolution of fiscal multipliers is inversely related to the evolution of actual-to-desired consumption ratio, and thus to credit rationing in the economy. However, the dynamics of fiscal multipliers appears to be lagging the one of credit rationing. This delivers an important results of the model, namely government expenditure is able to trigger variation in output insofar as it is able to let borrowers escape from the credit rationing trap and to get back to effective consumption levels compatible with their desired plans. Moreover, such a process of "re-leveraging" is faster the stronger is the level of government expenditure in the economy. This is showed by the fact that the peak multiplier is reached much faster for higher levels of government expenditure (see Figure 2). Figure 4: Deficit-spending vs. balanced-budget multipliers. In the plot "baseline" refers to the deficit-spending rule, "bal. budget 1" refers to the balanced-budget rule with endogenous government expenditure, "bal. budget 2" to the balanced budget rule with endogenous tax rate. Deficit-spending rule balanced-budget fiscal multipliers are always very low. On one hand, decreasing government expenditure to adapt it to falling tax revenues (as it is implied by the balanced-budget rule with endogenous government expenditure), does not help in triggering the re-leveraging process discussed before and restore pre-shock output levels. On the other hand rising the tax rate to fully balance government expenditure (as it is implied by the balanced budget rule with endogenous tax rate) further depresses household disposable incomes that are already falling after the shock and further reinforce credit rationing and de-leveraging in the economy.

To conclude, Figure 5 compares the evolution of fiscal multipliers (with the deficitspending rule) for different initial shares of borrowers in the households population. The latter is a measure of the degree of leverage in the economy as it captures the extent to which total consumption in the economy depends on consumption loans. The plot in the figure clearly shows that fiscal multipliers are growing in the degree of debt leverage in the economy. The same is observed for the magnitude and persistence of the effects of a bankruptcy shock (not shown). The result that multipliers are higher in more leveraged economy is explained by the fact that -to repeat -a greater share of consumption depends on bank loans in such economies and that bankruptcy shocks are going to generate higher degrees of credit rationing. It follows that the re-leveraging process triggered by deficitspending fiscal policy is going to be more pervasive in such economies and, thus, able to lead to higher increases in output in the aftermath of the shock.

Conclusions

We built a simple agent-based model to study how fiscal multipliers are affected by liquidity constraints in the economy. The model features two classes of heterogenous agents: savers, who can finance desired consumption with their liquid wealth, and borrowers who need bank's loans to satisfy their desired consumption plans. In the model bank's credit supply depends on bank's net worth that is affected by bad debt (if any). Moreover, the government taxes agents and purchases goods according to different fiscal rules. Finally, aggregate income is determined by aggregate demand. In such a framework we showed that small bankruptcy shocks can trigger wide and persistent falls in output. The magnitude and persistence of the shock is inversely related to the intensity of government expenditure in the economy. Moreover, we showed that fiscal multipliers can be significantly higher than one during recessions, and that the value of the multiplier is inversely related to the degree of credit rationing in the economy. Finally, we studied fiscal multipliers across different fiscal regimes and degrees of debt leverage in the economy. We showed that fiscal multipliers associated with balanced budget rules are very low compared to the ones associated with deficit-spending rules. Moreover, fiscal multipliers are higher in more leveraged economies. We plan to extend the model in several directions, by exploring further fiscal rules, and to the robustness of results under different consumption behaviors.

Figure 1 :

 1 Figure 1: Evolution of aggregate income as a fraction of steady state income. Deficitspending rule. The values in the legend correspond to different intensities of government expenditures and tax rates.

Figure 2 :

 2 Figure 2: Evolution of fiscal multipliers. Deficit-spending rule. All multipliers are calculated in relation to an intensity i = 0.02 of government expenditure (see also Eq. 13).

Figure 3 :

 3 Figure 3: Evolution of fiscal multipliers (left scale) and ratio between actual and desired consumption (right scale). Deficit-spending rule.

Figure 5 :

 5 Figure 5: Fiscal multipliers for different initial fractions of borrowers in the economy. Deficit-spending rule

  Debt is repaid. Household's wealth is determined. Household bankruptcy is determined3 Steady State ConditionsBefore carrying simulation experiments we identify the steady state conditions of the model. By steady state here, we mean a state of the economy where levels of all microeconomic variables (households wealth, households income, households consumption, debt, profits of the bank) are constant and where levels of all macroeconomic variables (aggregate consumption, aggregate government expenditure, tax revenues, aggregate income) are also constant.

• Total credit supply is determined. Credit is allocated to consumers. Actual private consumption is determined Aggregate income of the period is determined and distributed

• Taxes are collected •

  H * we get the following relations between individual consumption shares γ * i and individual marginal propensity to consume β

		1
	for a saver.
	Aggregate consumption must be constant in steady state. Thus, we can express steady
	state individual consumption, C * i as a fraction of steady state aggregate consumption C *
	N
	C * i = γ * i C * ,	γ * i = 1	(8)
	i=1
	Moreover, by exploiting the fact that, by definition, C * i = β * i w * i and that C

* = (1 -τ )Y