Deep numerical schemes for systems of Ergodic BSDEs with applications to regime-switching forward utilities

Unknown

Conceptio › Archive › arXiv CS

arXiv CSopen access

Deep numerical schemes for systems of Ergodic BSDEs with applications to regime-switching forward utilities

Unknown · 2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

Open Source ↗Direct PDF ↓

machine learning, deep learning, neural networks

Deep numerical schemes for systems of Ergodic BSDEs with applications to regime-switching forward utilities

arXiv:2606.24271v1 [math.NA] 23 Jun 2026

Guillaume Broux-Quemerais1

Sarah Kaakai2

Anis Matoussi1

Wissal Sabbagh1

Abstract In this paper, we introduce two neural-network-based numerical schemes for solving systems of coupled ergodic Backward Stochastic Differential Equations (eBSDEs), motivated by the approximation of optimal strategies within the framework of forward utilities in a regime-switching stochastic factor model. Our approach builds on the representation of such models through systems of eBSDEs introduced in [HLT20]. We first establish a link between the solution of the system of ergodic BSDEs and that of an associated multidimensional BSDE with random terminal time, given by the hitting time of the positive recurrent stochastic factor. Building on this representation, we introduce a locally additive deep learning scheme obtained by minimizing aggregated local error terms. We then present a new Deep Galerkin Method (DGM) inspired algorithm that minimizes the residual of the associated ergodic PDE system, relying on a representation of the ergodic cost. Finally, we apply this framework to regime-switching forward utilities in a stochastic factor model. We first derive a general consistency SPDE that characterizes regime-switching forward utilities and retrieve their representation with systems of ergodic BSDEs in the homothetic case. Numerical experiments demonstrate the performance of the proposed methods, with a particular focus on the impact on forward preferences of taking into account regime switches.

∗

Acknowledgements: The authors research is part of the ANR project DREAMeS (ANR-21-CE46-0002) and benefited

from the support of respectively the "Chair Risques Emergents en Assurance" and "Chair Impact de la Transition Climatique en Assurance" under the aegis of Fondation du Risque, a joint initiative by Risk and Insurance Institute of Le Mans, and MMA-Covéa and Groupama respectively. 1

Laboratoire Manceau de Mathématiques & FR CNRS No 2962, Institut du Risque et de l’Assurance, Le

Mans Université 2

Laboratoire Analyse, Géométrie et Applications, UMR CNRS 7539, Université Sorbonne Paris Nord

Introduction In this paper, we are interested in the numerical approximation of certain classes of forward performance/utility processes in a regime-switching financial market, and their associated optimal decision criterion. Introduced by [MZ06], forward utilities offer an interesting alternative to the classical setting of expected utility maximization at a terminal time. This forward-looking approach enables the dynamic adjustment of the decision criteria, starting from preferences which are known at an initial time, rather than imposing a potentially distant and arbitrary time horizon. In a continuous setting, the agent’s preferences are represented by a utility random field U (t, ·). This preference criterion maintains time consistency within the given investment or decision-making context, in the sens that if Xtπ is the observable process (typically the wealth), resulting from the admissible strategy π, then the preference process U (t, Xtπ ) is a supermartingale, and there exists an optimal strategy such that the preference process is a martingale. Since their introduction, there have been significant theoretical advancements in the field of forward utilities. In a general framework, [MZ10] established a sufficient condition for time-consistency when the utility random field follows dynamics of Itô type. This condition takes the form of a nonlinear SPDE of Hamilton-Jacobi-Bellman (HJB) type satisfied by U . By establishing a correspondence between the SPDE and the compound of two SDEs, the authors in [EKM13] provide a method to construct forward utilities. This work has been extended to consistent utility of investment and consumption in [EKHM18], and has been applied, for instance, to derive consistent utilities in stochastic factor market models (see e.g. [NZ14], [ASS20], [LZ17]). To the best of our knowledge, the study of forward utilities in regime-switching markets has only been tackled in [HLT20]. Forward utilities have found diverse applications over recent years, including but not limited to option valuation, insurance, mean field games ([LZ19], [DRP21]), long term interest rate modeling ([EKHM22]), risk measures ([Cho19] or more recently pension design ([HKM24],[NC24]). Surprisingly, the subject of numerical methods for forward preferences remains largely unexplored, despite its critical importance for practical applications. In [GM18], a general approach is proposed using strong approximations of compounds of random maps. In this paper, we focus on a different point of view, introduced in [HLT20], to develop new numerical methods for forward utilities in regime-switching markets, based on a representation with means of systems of ergodic BSDEs. Regime-switching models, also called Markov-modulated dynamics, are designed to capture financial markets that exhibit multiple states. In addition to the standard Brownian motion, which drives stock price dynamics, these models incorporate an additional source of randomness represented by a finite-state continuous-time Markov chain, where each state corresponds to a distinct market regime. For instance, a two-regime model may describe a "bull market," characterized by rising asset prices, and a "bear market," marked by declining prices. In this framework, the Brownian motion captures persistent microeconomic effects, while the Markov chain accounts for occasional, macroeconomic shifts in market conditions. Regime-switching models have been extensively studied in the context of option pricing; see, for example, [BE02], [GZ04], [YZZ06], [JR06]. Portfolio optimization problems within regime-switching frameworks have also been considered, including mean-variance portfolio selection in [ZY03] and utility maximization across various models, see [Zar92], [BR04], [SC09], and [FWY14]. Our motivation is the representation of an agent’s preferences in an incomplete financial market 2

presenting multiple states. Formally, the regime transitions are governed by a continuous-time Markov chain α, with a finite state space I = {1, . . . , I}. In this context, the decision criteria are built upon a family of utility random fields U i i∈I , each representing the agent’s preferences in regime i. The regime-switching utility is the càdlàg process defined as U (t, x) = U αt (t, x) =

U i (t, x)1{αt =i} .

i∈I

We focus on a stochastic factor model in which the dynamics of stock prices are influenced by a stochastic factor V . This leads to the study of regime-switching preferences based on families of homothetic utilities, which are expressed as separable functions of the stochastic factor. For each i

i ∈ I, these are given by U i (t, x) = u(x)ef (t,Vt ) , where u is a standard exponential or power utility function (with an additive expression in the logarithmic case). While the main focus of [HLT20] is on existence and uniqueness of the Markovian solution to systems of eBSDEs, they also characterize the functions f i as solutions of such systems for homothetic utilities in power form. Systems of ergodic BSDEs generalize the notion of ergodic BSDEs introduced in [FHT09] (see also [BP99] with a different formulation with stationary processes ), [DHT11], [HL19]), to the case of multidimensional and coupled components Y . Informally, a system of ergodic BSDEs with generators (F i )1≤i≤I and coupling terms (Gi )1≤i≤I is a system of backward stochastic differential equation on an infinite horizon, whose solution is a triplet ((Y i , Z i )1≤i≤I , λ) where λ ∈ R and for all i = 1, ..., I, processes Y i , Z i are adapted and satisfy for any T > 0, P-a.s for any 0 ≤ t ≤ T : Yti

YTi +

Z T

F (Vs , Zsi ) + Gi (Ys ) ds − λ(T − t) −

Z T

(Zsi )⊤ dWs .

(0.1)

Existence and uniqueness of a Markovian solution to systems of eBSDEs of the form (0.1) are established in [HLT20], under assumptions analogous to those for classical ergodic BSDEs: locally Lipschitz conditions on F i and strong dissipativity of the stochastic factor. Such systems are genuinely multidimensional quadratic BSDEs, for which existence and uniqueness have been extensively studied, see [FDR11] [Tev08, Fre14, HT16, Xv18, HR19]. In [HLT20], the quadratic growth is handled by a truncation argument exploiting a specific property of the solution, namely the boundedness of the Z component. Numerical methods for ergodic BSDEs remain, by contrast, largely unexplored. The authors of [GRS24] propose a fully implementable scheme with an approximation-error analysis based on nonlinear Feynman–Kac formulas. In a different direction, machine-learning methods for high-dimensional nonlinear PDEs and BSDEs have been investigated extensively in recent years [HJ+ 17, HPW20, GPW+ 23, KT24], and were first extended to the ergodic setting in [BQKMS24]. We pursue this line for systems of ergodic BSDEs arising from regime-switching forward utilities. Objective and contributions. The aim of this paper is to develop numerical methods for the approximation of Markovian solutions to systems of ergodic BSDEs, with a particular focus on those arising from homothetic regime-switching forward utilities. Our first contribution is theoretical. We state rigorously the equivalence between the system of

ergoidc BSDEs, the ergodic PDE system and the ergodic BSDE with jumps. These correspondences underpin the foundations of our algorithms. Computing the dynamics of U (t, Xtπ ) by means of a generalized Itô–Ventzell formula with jumps [ØZ07, MM22], we derive a consistency SPDE with jumps (Theorem 3.1) that characterizes regime-switching forward utilities in our framework, providing a general consistency condition in regime-switching models. Applied to homothetic preferences, it identifies the systems of ergodic BSDEs associated with power, exponential and logarithmic utilities, recovering in particular the power-utility system of [HLT20]. Our second contribution is numerical. We introduce two neural-network schemes for the simultaneous approximation of the ergodic cost λ and the processes Y and Z: a locally additive deep BSDE scheme (LAeBSDE) and an extension of the Deep Galerkin Method (DGM). The probabilistic foundation of the first scheme rests on a random-horizon reformulation of the ergodic system. Since the initial agent’s utility is known, an initial condition Y0α0 ∈ R is naturally available. Leveraging the recurrence of the stochastic factor V , we introduce the first return time τ of V to its initial point v0 , which is almost surely finite under the dissipativity assumption and satisfies Yτ = Y0 . We show that the solution of the ergodic system (0.1) then coincides with that of an associated system of BSDEs with random terminal time τ , and, under suitable integrability of τ , that this solution is unique, building on [Par98] and on the fact that λ is pinned down by the bridge condition Yτ = Y0 . The argument is more involved than in [BQKMS24] since the normalization only fixes one point v0 in one state i0 of the solution. Building on this reformulation, the LAeBSDE scheme discretizes V by an Euler scheme, estimates the horizon τ along each trajectory, and treats λ as a trainable parameter, optimized jointly with two networks approximating Y and Z through an aggregation of local losses. The DGM scheme follows a distinct route by approximating the Markovian solution directly through the minimization of the residual of the associated ergodic PDE system at randomly sampled points, the ergodic cost being recovered from its invariant-measure representation. Finally, we assess the performance of the two algorithms on several examples. We first construct a class of examples admitting an explicit solution, by starting from a sublinear ansatz and choosing the market price of risk accordingly. We then study a two-state market model with regime-switching power utility and an affine market price of risk, where a comparison with the non-switching case highlights the economic implications of our framework for forward preferences. The paper is organized as follows. In Section 1, we recall the framework and the existence and uniqueness results associated with the system of ergodic BSDEs. We also state the ergodic PDE and BSDE with jumps associated with the system of ergodic BSDE. In Section 2, we establish a correspondence between systems of eBSDEs and BSDEs with random terminal time, which is needed for the definition of the locally additive scheme. We then extend the Deep-Galerkin method to the approximation of ergodic PDEs. Section 3 introduce the regime-switching stochastic factor market model and the notion of forward utilities. We then derive the consistency SPDE that allows the identification of systems of ergodic BSDEs associated with power, exponential and logarithmic homothetic utilities. Finally, in Section 4 we present numerical results on one benchmark example and interpret financial implications of regime-switches through an example of power utility with affine market price of risk.

Notations: All stochastic processes in the sequel are defined on a standard probability space (Ω, F, F, P), on which we consider a d-dimensional Brownian motion W , and denote F = (Ft )t≥0 its augmented filtration. 1

For x ∈ Rd , we denote x⊤ the transpose of vector x, ∥.∥ the usual norm ∥x∥ = Tr(xx⊤ ) 2 , dist(x, Π) the distance function of x to a closed convex subset Π ⊂ Rd and diag(x) = (diag(x)ij )1≤i,j≤d ∈ Rd×d , the matrix such that diag(x)ii = xi and diag(x)ij = 0 for i ̸= j. For X = (Xt )t≥0 a càdlàg stochastic process, we denote Xt− the left-limit of this process at time t. For any integer k ≥ 1 and positive constant B ∈ R+∗ , let φB denote the projection on the centered ball of Rk of radius B. We denote L2 the space of square integrable random variables and also introduce the usual spaces of solution for γ ∈ R and τ a F stopping time: S 2 (γ, τ )

2 (φt )t≥0 , Rk -valued progressively measurable process s.t. E sup eγs ∥φs ∥ < ∞.

Z τ 2 eγs ∥φs ∥ < ∞. (φt )t≥0 , Rk -valued progressively measurable process s.t. E

0≤s≤τ

M(γ, τ )

System of coupled ergodic BSDEs

We are interested in the numerical study of coupled systems of ergodic BSDEs, as introduced in [HLT20]. Let I > 1, and denote I = {1, . . . , I}. In the following, we call systems of ergodic BSDEs with generators (F i )i∈I and coupling terms (Gi )i∈I , a multidimensional BSDE on an infinite horizon, whose solution is a triplet ((Y i , Z i )i∈I , λ) where λ ∈ R and (Y i , Z i )1≤i≤I is a family of F-adapted processes taking values in R × Rd , satisfying for any T > 0 and P-a.s for any 0 ≤ t ≤ T : Yti = YTi +

Z T

F i (Vsv0 , Zsi ) + Gi (Ys ) ds − λ(T − t) −

Z T

(Zsi )⊤ dWs ,

(1.1)

where V v0 is a stochastic factor solution of ′

dVtv0 = µ(Vtv0 )dt + κ⊤ dWt ,

V0v0 = v0 ∈ Rd ,

′

κ ∈ Rd×d ,

(1.2)

such that ′

Assumption 1.1. µ is Lipschitz, and there exists a constant Cµ > 0 such that for any v, v̄ ∈ Rd (µ(v) − µ(v̄))⊤ (v − v̄) ≤ −Cµ ∥v − v̄∥2 .

(1.3)

and κκ⊤ is a positive definite matrix. In particular, the diffusion V v0 is exponentially ergodic. The next lemma provides an estimate on the running supremum of the stochastic factor V v0 , which will be used in Section 2.1 to establish integrability properties of the Markovian solution to the system of ergodic BSDEs. The proof is postponed to Appendix B. Lemma 1.1. Suppose that Assumption 1.1 holds and let V v0 be the solution of (1.2). Then for every q ≥ 1, there exists a constant Cq > 0, depending only on q, Cµ and κ, such that for every stopping

time τ , v0 2q E sup |Vt | ≤ Cq 1 + |v0 |2q + E[τ q ] .

(1.4)

0≤t≤τ

1.1

Existence and uniqueness of Markovian solution

We consider an extension of the setting studied in [HLT20], allowing for a more general coupling term (Gi )1≤i≤I , defined in (1.7), instead of the specific coupling function g(y) = ey − 1 considered therein. Under Assumption 1.2 (iii), the proof follows essentially the same arguments as in Theorem 3.2 of [HLT20], and the corresponding existence and uniqueness result for bounded Markovian solutions to (1.1) remains valid. We recall this result below under the additional assumptions imposed in our framework. Assumption 1.2. The following properties hold: ′

i) There exists a positive constant K such that ∀v ∈ Rd , |F (v, 0)| ≤ K. ′

ii) There exists positive constants Cv and Cz such that for i ∈ I, ∀v, v̄ ∈ Rd , ∀z, z̄ ∈ Rd×I F i (v, z) − F i (v̄, z) i

F (v, z) − F (v, z̄)

≤ Cv (1 + ∥z∥)∥v − v̄∥,

(1.5)

≤ Cz (1 + ∥z∥ + ∥z̄∥)∥z − z̄∥.

(1.6)

Moreover, we require that Cv < Cµ with Cµ as introduced in (1.3). iii) There exists an increasing C 1 function g : R → R, such that, for all i ∈ I and y ∈ RI Gi (y) =

I P q ij g y j − y i ,

(1.7)

j=1

g(0) = 0 and for all y ∈ RI , g(y) − g(−y) ≥ y. Note that, under Assumption 1.2 (iii), the function g is Locally Lipschitz on R. Hence, setting P Lg (M ) = sup g ′ (x) and q ∗ = max q ij , the map G : RI → RI is Lipschitz on |x|≤M

i∈I j̸=i

DM := y ∈ RI , y j − y i ≤ M, ∀ i, j ∈ I . Indeed, for any y, ȳ ∈ DM and any i ∈ I, Gi (y) − Gi (ȳ)

q ij g(y j − y i ) − g(ȳ j − ȳ i )

j̸=i

≤ Lg (M )

q i (y j − y i ) − (ȳ j − ȳ i )

j̸=i

q ij y j − ȳ j + y i − ȳ i

j̸=i

≤ Lg (M )2∥y − ȳ∥

X j̸=i

≤ Lg (M )2∥y − ȳ∥q ∗ . 6

q ij

and hence, for all i ∈ I, Gi is also locally Lipschitz on R. In addition to the usual ergodicity of the stochastic factor V v0 , ensured by Assumption 1.1, the study of Markovian solutions to systems of ergodic BSDEs (1.1) also requires the irreducibility of the transition rate matrix associated with the coefficients q ij . Assumption 1.3. Q = (qij )1≤i,j≤I is a transition rate matrix, i.e. , q ii = −

P ij q , for all i ∈ I.

j̸=i

Furthermore, we assume that for all i ̸= j q ij > 0 and we denote qmin = min q ij . i̸=j

′

Theorem 1.2. Let i0 ∈ I, v0 ∈ Rd and y0 ∈ R. Under Assumptions 1.1, 1.2, 1.3, there exists a unique Markovian solution (y i (Vtv0 ), z i (Vtv0 ))i∈I , λ t≥0 to the system of ergodic BSDEs (1.1), such ′

that y i0 (v0 ) = y0 ∈ R and for all i ∈ I, and v ∈ Rd , y i (v) z i (v) y i (v) − y j (v) Remark 1.1.

≤ C(1 + ∥v∥) Cv ≤ ∥κ∥ Cµ − Cv Cv Cµ Cz 1 ≤ K+ := CY . qmin (Cµ − Cv )2

(1.8) (1.9) (1.10)

1. Notably, the constant λ in the solution of (1.1), called the ergodic cost, is shared

by all equations for all i ∈ I, and do not depend on the choice of i0 , v0 or y0 . 2. A further important point is that uniqueness of Markovian solutions is ensured by fixing only one component, y i0 (v0 ), of the vector-valued solution y = (y 1 , . . . , y I ) at i0 , v0 . Hence, one cannot fix the whole vector (y 1 (v0 ), . . . , y I (v0 )) simultaneously. This constitutes a nontrivial obstacle for the numerical approximation of the solution. 3. The assumption of sublinearity of the solution components y i is required for uniqueness in Theorem 1.2. In fact, [HLT20] provide examples of solutions to ergodic BSDEs that do not satisfy this growth property. Link with semilinear PDE The system of Makovian ergodic BSDEs (1.1), represented by (y i (Vtv0 ), z i (Vtv0 ))i∈I , λ t≥0 , gives a probabilistic representation of the following system of elliptic PDE, for all i ∈ I Ly i (v) + F i (v, ∇y i (v)κ) +

q ij g(y j (v) − y i (v)) = λ

(1.11)

j∈I

where L denotes the generator of the semi-group associated to the SDE (1.2) Proposition 1.3. Under Assumptions 1.1, 1.2, 1.3, (y i )i∈I is a viscosity solution of (1.11). Proof. We consider, for all ρ > 0, the following infinite horizon markovian BSDE system: for t ≥ 0

and i ∈ I Yti,ρ,v = YTi,ρ,v +

Z T





F i (Vsv , Zsi,ρ,v ) +

q ij g(Ysj,ρ,v − Ysi,ρ,v ) − ρYsi,ρ,v ds −

Z T

(Zsi,ρ,v )⊤ dWs .

j∈I

(1.12) ′

and define for a fixed reference point, say v0 ∈ Rd , y i,ρ (v) := Y0i,ρ,v , ȳ i,ρ (v) := Y0i,ρ,v − Y0α0 ,ρ,v0 , for ′

all v ∈ Rd . By standard arguments, see e.g. the proof of Theorem 5.74 in [PR14], ȳ i,ρ is a viscosity solution of the elliptic PDE X ij Ly i (v) + F i v, ∇y i (v)κ + q g(y j (v) − y i (v)) = ρy i (v) + ρy α0 ,ρ (v0 )

′

v ∈ Rd .

(1.13)

j∈I

Following the arguments of [HLT20], we can prove that there exists a sequence (ρn )n∈N such that ′

ρn ↘ 0, ρn y α0 ,ρn (v0 ) → λ, and ȳ i,ρn (v) → u uniformly on Rd as n → +∞. Then, Remark 6.3 in [CIL92] implies that y i is a viscosity solution of (1.11).

Link with eBSDE with jumps A useful representation of the solution of (1.1) introduced in [HLT20] relies on the transformation of the system of ergodic BSDEs (1.1) into an ergodic BSDEs with jumps, where jumps occur at jump times of a continuous time Markov chain (CTMC) on I, of transition rate matrix Q = (qij )1≤i,j≤I . Assume that the probability space supports a random Poisson measure N (dt, di, dj) on R+ × I 2 , of intensity qij dtγ(di)γ(dj), with γ the counting measure on I. Let α be the CTMC of transition rate matrix Q, solution of Z tZ αt = α0 + I2

(j − i)1{αs− =i} N (ds, di, dj),

In the following, we will denote F αt (.) =

∀t ≥ 0.

i αt i=1 F (.)1{αt =i} and G (.) =

(1.14)

i i=1 G (.)1{αt =i} .

Lemma 1.4. Let (y i (Vtv0 ), z i (Vtv0 ))i∈I , λ t≥0 be the Markovian solution of the system of ergodic BSDEs (1.1), as defined in Theorem 1.2. Then, (y αt (Vtv0 ), z αt (Vtv0 ), ψt (·, ·, Vtv0 ), λ)t≥0 is solution of the following ergodic BSDE with jumps:

Z T Yt = YT +

 F αs (Vsv0 , Zs ) + Gαs (Ys ) −

I X

 q αs− ,j ψs (αs− , j, Vtv0 )ds

j=1

Z T − λ(T − t) −

Z TZ

⊤

(Zs ) dWs − t

ψs (i, j, Vtv0 )1{αs− =i} Ñ (ds, di, dj), (1.15)

where ψt (i, j, Vtv0 ) := y j (Vtv0 ) − y i (Vtv0 ), and Ñ the compensated Poisson measure N (dt, di, dj) − qij dtγ(di)γ(dj). 8

(1.16)

1.2

Representation of the ergodic cost

Finally, the ergodic cost λ admits a representation of a weighted space-time average with respect to the invariant measures of the recurrent state Markov chains Vtv0 . A similar identity was used in the fix -point algorithm of [GRS24]. It may be used to obtain an approximation of the ergodic cost provided that we know or approximate correctly the invariant distribution ν. Proposition 1.5. Let ν be the unique invariant probability measure of the diffusion V v0 with infinites imal generator L. Let y i (.), z i (.), λ denote the unique Markovian solution to the system of ergodic BSDEs (1.1) such that y α0 (v0 ) = y0 . Then the ergodic cost λ satisfies, for all i ∈ I !

Z λ=

Rd′

F i (v, z i (v)) +

(1.17)

q ij g(y j (v) − y i (v)) ν(dv).

i∈I

Proof. The Markovian solution to the system of ergodic BSDEs (1.1) verifies the following equation  y i (v0 ) = E y i (VTv0 ) +

Z T

 F i (Vsv0 , z i (Vsv0 )) +

I X

  q ij g(y j (Vsv0 ) − y i (Vsv0 )) − λds . (1.18)

j=1

Integrating with respect to the measure ν yields Z

Rd′

y (v0 )ν(dv0 ) =

′ Rd

Z +

Rd′

E y i (VTv0 ) ν(dv0 )     Z T I X F i (Vsv0 , z i (Vsv0 )) + E q ij g(y j (Vsv0 ) − y i (Vsv0 )) − λds ν(dv0 ). 0

j=1

Then, we can conclude (1.17) by applying Fubini theorem and the definition of the invariant probability measure ν.

Deep learning algorithms for the simulation of systems of eBSDEs

In this section, we focus on the numerical simulation of the Markovian solution to the system of ergodic BSDEs introduced above. We present two deep learning-based approaches. The first is a deep BSDE scheme of locally additive type, in which the solution and its gradient are represented by neural networks trained against a sum of local losses, each accumulating the dynamics up to the terminal condition, following the methodology of LaDBSDE in [KT24]. The second approach relies on a Galerkin-type approximation, see [SS18]. Deep learning methods for BSDEs have attracted considerable attention in recent years, mainly because of their ability to handle high-dimensional nonlinear PDEs through probabilistic BSDE representations; see, for instance, [HJ+ 17], [CWNMW19], [HPW20], [GPW+ 23], and [KT24]. In the present framework, we adapt these ideas to systems of ergodic BSDEs, in the spirit of the approach developed in [BQKMS24]. We first establish a connection between the regime-switching system of eBSDEs and a system of eBSDEs with random terminal time. This representation will 9

serve as the probabilistic foundation for the deep BSDE scheme introduced below. The Galerkin-type method follows a different numerical strategy and is presented independently afterwards.

2.1

Connection with a system of eBSDEs with random terminal time

In the spirit of [BQKMS24], we exploit the recurrence property of the stochastic factor V v0 to establish a correspondence between the system of ergodic BSDEs (1.1) and a multidimensional BSDE with random terminal time. Unless stated otherwise, we assume throughout this section that the stochastic factor is one-dimensional, namely d′ = 1. We fix a minimal deterministic horizon T0 and define the random horizon τ as the first return time after T0 of the diffusion V to its initial value v0 : τ = inf {t > T0 , , Vtv0 = v0 }.

(2.1)

Let i0 ∈ I and v0 , y0 ∈ R. Under Assumptions 1.1 and 1.2, there exists a unique solution i (y (Vtv0 ), z i (Vtv0 ))i∈I , λ t≥0 to the system of ergodic BSDEs (1.1) such that y i0 (v0 ) = y0 . Moreover, for all i ∈ I, the function y i has sublinear growth with respect to v, the function z i is bounded v , and the differences y i − y j are uniformly bounded. by Zmax := ∥κ∥ CµC−C v

By construction, (y i (Vtv0 ), z i (Vtv0 ))i∈I , λ t≥0 also solves the following “ergodic” BSDE with random terminal time and the same initial and terminal condition: for each component i = 1, . . . , I, Ytr,i = Yτr,i + Yτr,i = Y0r,i ,

Z τ

F (Vsv0 , Zsr,i ) + Gi (Ysr ) ds − λ(τ − t) −

Z τ

(Zsr,i )⊤ dWs ,

(2.2)

Yτr,i0 = y0 .

Here, τ is the return time defined in (2.1). The term “ergodic” BSDE with random terminal time is a slight abuse of terminology, since the equation is no longer formulated on an infinite time horizon. We use it to emphasize that the ergodic cost λ remains part of the unknowns in (2.2). The infinite horizon in (1.1) is replaced by the random terminal time τ , together with the additional constraint that the initial and terminal values coincide, namely Y0r,i = Yτr,i , i ∈ I. Conversely, solutions (Y r,i , Z r,i )i∈I , λ t∈[0,τ ] of (2.2) can be studied in their own right. The main difficulty, compared with the one-dimensional case, is that the normalization cannot be imposed componentwise. Indeed, only the value of the solution at one point v0 and in one state i0 , say Y0r,i0 = Yτr,i0 = y0 , can be fixed. On the other hand, the remaining initial values Y0r,i , for i ̸= i0 , cannot be chosen due the coupling term G = (Gi )i∈I and the shared parameter λ. Hence, the random-time formulation does not reduce to a collection of independent scalar BSDEs, the coupling must be handled simultaneously with the ergodic constant λ and the periodic-type constraint Y0r,i = Yτr,i for all i ∈ I. The following theorem states the existence and uniqueness of a Markovian solution to (2.2). It ensures that this solution indeed coincides on [0, τ ] with the unique solution of the system of ergodic BSDEs (1.1). Recall that CY is defined in (1.10). Let Kz := Cz (1 + 2Zmax ) be the Lipschitz constant with respect to z of the truncated driver F ◦ φZmax , and let KG be the Lipschitz constant of maxGi on the i∈I set x ∈ RI ; |xi − xj | ≤ CY , ∀ i, j ∈ I . 10

Theorem 2.1. Suppose that Assumption 1.1 and 1.2 and 1.3 hold. Assume that the random horizon τ admits an exponential moment of order Γ > 2KG + Kz2 . Then, the “ergodic” BSDE with random horizon (2.2) with Yτr,i = Y0r,i for all i ∈ I, and Yτr,i0 = y0 , admits a unique Markovian solution ((y i (.), z i (.))i∈I , λ), such that for any γ < Γ, y(Vtv0 ) ∈ S 2 (γ, τ ), z(Vtv0 ) ∈ M(γ, τ ). In addition, if for all i, j ∈ I, y i is sublinear, differences y i − y j are uniformly bounded by CY , and z i is bounded by Zmax , then the solution ((y i (.), z i (.))i∈I , λ) coincides on [0, τ ] with the unique Markovian solution of the system of eBSDEs (1.1), as defined in Theorem 1.2, verifying y i0 (v0 ) = y0 . Remark 2.1. Sufficient conditions for the exponential integrability of the hitting time τ can be obtained as a straightforward extension of Lemma 3.1 in [BQKMS24]. Proof. Existence - Let γ < Γ. By construction, the unique solution (y i (Vtv0 ), z i (Vtv0 )))i∈I , λ t≥0 of the system of ergodic BSDEs (1.1) is also solution of the BSDE with random terminal time (2.2), with z = (z i )i∈I bounded by Zmax so that Z ∈ M(γ, τ ). For all i ∈ I, y i has sublinear growth. Hence, there exists a constant C > 0 such that for all i ∈ I v0 2 2 γτ 2 γτ i 2 γτ ≤ 2C E[e ] + 2C E e sup |Vt | . E e sup Yt 0≤t≤τ

0≤t≤τ

Under the assumption that τ admits exponential moments of order Γ > γ, the first expectation on the p be conjugate exponents, such that pγ < Γ. By Hôlder’s right side is finite. Let p > 1 and q = p−1

inequality, 1 q 1 v0 2 v0 2q γτ pγτ p E e sup |Vt | ≤ E([e ]) E sup |Vt | . 0≤t≤τ

0≤t≤τ

The first factor is finite by assumption, while the second one is finite by Lemma 1.1. Y ∈

In turn

S 2 (γ, τ ).

r,i r,i Uniqueness of λ - Let (Y r,i , Z r,i )i∈I , λ t∈[0,τ ] and (Y , Z )i∈I , λ

t∈[0,τ ]

be two Markovian solu-

tions of the ergodic system of BSDEs with random terminal time (2.2). Let α be the CTMC defined in (1.14). This will serve only as a representation tool for uniqueness and without loss of generality, we can assume it starts from any state l ∈ I. Define for all 0 ≤ t ≤ τ , r

r,αt

), and

(Ytr , Ztr ) := (Y r,αt , Z r,αt ),

(Y t , Z t ) := (Y

ψtr (i, j, Vtv0 ) := Ytr,j − Ytr,i ,

ψ t (i, j, Vtv0 ) := Y t − Y t .

,Z r,j

r,i

Then, using the same notation as in Lemma 1.4, (Y r , Z r , ψ r , λ) and (Y , Z , ψ , λ) are solutions on [0, τ ] of dYt = −(F

αt

(Vtv0 , Zt ) + Gαt (Yt ))dt + λdt + (Zt )⊤ dWt + r

Z I2

ψt (i, j, Vtv0 )1{αt− =i} N (dt, di, dj),

∀0 ≤ t ≤ τ.

Let ∆Ytr = Ytr − Y t , ∆Ztr = Ztr − Z t , ∆λ = λ − λ̄, and ∆ψtr (i, j) = ψtr (i, j) − ψ t (i, j). The equation verified by ∆Y r can be linearized by using a change of measure as in the proof of Theorem 3.2 in [HLT20]. Indeed, using the previous equation and by definition of (Gi )i∈I in Assumption 1.2 11

(iii), we have ∀ 0 ≤ t ≤ τ r

d∆Ytr = −(F αt (Vtv0 , Ztr ) − F αt (Vtv0 , Z t ))dt + ∆λdt −

I X

(2.3)

r q αt− ,j g(ψtr (αt− , j)) − g(ψ t (αt− , j)) − ∆ψtr (αt− , j) dt

j=1

+ (∆Ztr )⊤ dWt +

∆ψsr (i, j)1{αt− =i} Ñ (dt, di, dj) (2.4) Z = ∆λdt + (∆Ztr )⊤ (dWt − γtF dt) + ∆ψtr (i, j)1{αt− =i} Ñ (dt, di, dj) − γtg (i, j)dtγ(di)γ(dj) , I2

with αt (V v0 , Z r,αt− ) − F αt (V v0 , Z r,αt− ) t t t t ∆Ztr 1{∆Ztr ̸=0} , and ∥∆Ztr ∥2 r g(ψtr (i, j)) − g(ψ t (i, j)) − ∆ψtr (i, j) γtg (i, j) = q i,j 1{∆ψtr (i,j)̸=0} . ∆ψtr (i, j)

F γtF =

By Assumption 1.2 and the boundedness of Z r , Z , (Y r,i −Y r,j ) and (Y

r,i

−Y

r,j

), the processes γtF and

(γtg (i, j))i,j∈I 2 are uniformly bounded. Hence, by Girsanov’s theorem, for each initial state α0 = l ∈ I r − ∆λ(t ∧ τ )) and all T ≥ 0, there exists an equivalent probability measure P̃l under which (∆Yt∧τ 0≤t≤T

is a martingale (for ease of notations, we omit the v0 in P̃l ). Consequently

∆Y0r,l = Ẽl [∆YTr∧τ ] − ∆λẼl [T ∧ τ ]. = Ẽl [∆Yτr 1T ≥τ ] + Ẽl [∆YTr 1τ >T ] − ∆λẼl [T ∧ τ ]. Since each solution satisfies Yτr,i = Y0r,i , for all i ∈ I, the differences also satisfy ∆Yτr,i = ∆Y0r,i . Therefore ∆Y0r,l =

P̃l (τ < T, ατ = k)∆Y0r,k + Ẽl [∆YTr 1τ >T ] − ∆λẼl [T ∧ τ ].

(2.5)

k∈I

Using the sublinear growth property of ∆Y , there exists a constant C such that 1

Ẽl [∆YTr 1τ >T ] ≤ Ẽl [|∆YTr |2 ] 2 P̃l (τ > T ) 2 2

≤ C(1 + Ẽl [ VTv0 ]) 2 P̃l (τ > T ) 2 .

(2.6)

Under P̃l , the process (V v0 , α) satisfies dVtv0 = µ(Vtv0 ) + γtF dt + κ⊤ (dWt − γtF dt) Z dαt = (j − i)1{α− =i} N (ds, di, dj), I2

with W −

R·

F 0 γs ds is a Brownian motion under P̃l

γsg (i, j))s≥0,(i,j)∈I 2 under P̃l . 12

and N is a point process with intensity (q ij +

Moreover, γsF can be written as a uniformly bounded function of Vsv0 . By Theorem 2.6 and Lemma 3.4 in [DHT11], this yields that V v0 is recurrent under P̃l . Hence P̃l (τ > T ) −−−−→ 0. T →∞

Furthertmore, by Proposition 5 in [HL19], sup Ẽl (VTv0 )2 < ∞.

T ≥0

Thus, taking the limit as T goes to infinity in (2.6) leads Ẽl [∆YTr 1τ >T ] −−−−→ 0. T →∞

By the equivalence of Pl and P̃l on FT , P̃l (τ < T, ατ = k) > 0 ⇐⇒ Pl (τ < T, ατ = k) > 0 By recurrence of (V v0 , α) on Pl , Pl (τ < T, αT ∧τ = k) > 0, for T large enough. By monotone convergence, this yields that P̃l (τ < T, ατ = k) ↑ P̃l (τ < ∞, ατ = k) > 0. Let P̃ = (P̃lk )1≤l,k≤I := (P̃l (τ < ∞, ατ = k))k,l∈I . Then, P̃ is a strictly positive stochastic matrix on the finite state space I. It therefore admits a unique invariant measure π̃ = (π̃ l )l∈I . Combining this with (2.5) and letting T → ∞, we obtain ! X

π̃l ∆Y0r,l =

l∈I

X X k∈I

π̃l P̃lk

∆Y0r,k − ∆λ

π̃l Ẽl [τ ].

l∈I

Since π̃ P̃ = π̃, this becomes X l∈I

Hence, ∆λ

π̃l ∆Y0r,l =

π̃k ∆Y0r,k − ∆λ

k∈I

π̃l Ẽl [τ ].

l∈I

π̃l Ẽl [τ ] = 0. Finally, since Ẽl [τ ] > T0 for every l ∈ I, we conclude that ∆λ = 0.

l∈I

Uniqueness of terminal values - Since ∆λ = 0, identity (2.5) rewrites in vector form ∆Y0r = P̃ ∆Y0r , where ∆Y0r = ∆Y0r,l

1≤l≤I

. Thus ∆Y0r is a P̃ -harmonic function on the finite state space I. Since P

is irreducible, Lemma 1.16 in [LP17] implies that ∆Y0r is constant on I. By (2.2), Y r,i0 = Ȳ r,i0 = y0 . Hence, for every l ∈ I, ∆Y0r,l = ∆Y0r,i0 = 0. Finally, since ∆Yτr,l = ∆Y0r,l by (2.2), we also get that for all l ∈ I, ∆Yτr,l = 0. 13

(2.7)

Uniqueness of Markovian solution (y i (.), z i (.))i∈I - Define for all i ∈ I, y ∈ RI , a truncated version of the coupling term in the generator, namely GiCY (y) =

q ij g ◦ φCY (y j − y i ).

j̸=i

The rest of the proof follows by applying Theorem 3.1 in [Par98] to the BSDE with random terminal time τ , terminal condition Yτr (unknown but unique from the last step) and generator F ◦ φZmax (v, z) + GCY (y) − λ, where λ is fixed. Under Assumption 1.2, the generator F ◦ φZmax + GCY is Lipschitz in z and satisfy the usual monotonicity condition in y, with monotonicity constant KG . By assumption, the random horizon τ admits exponential moments of order Γ > 2KG +Kz2 . Theorem 3.1 in [Par98] applies and ensures existence and uniqueness of a Markovian solution ((Y r,i , Z r,i )i∈I ) ∈ S 2 (γ, τ ) × M(γ, τ ). From the first step of the proof, this solution coincides necessarily with the unique Markovian solution (y i (Vt ), z i (Vt ))i=1,...,I , λ t≥0 of the system of ergodic BSDEs (1.1), which concludes. The solution of the system of ergodic BSDE (1.1) thus coincides with the Markovian solution (Y r,i , Z r,i )i∈I , λ of the ergodic BSDE with random time horizon and fixed initial condition (2.2) on [0, τ ]. We will omit the subscript r in the sequel. This point of view provides our simulation problem with a random horizon τ allowing to adapt numerical schemes introduced in [BQKMS24].

2.2

Locally additive deep solver

We can now introduce the locally additive deep solver for approximating the solution (Yti , Zti )i∈I , λ t≥0 of the random time horizon ergodic BSDE (2.2), associated with the forward stochastic factor given by (1.2). The main idea is to minimize a global loss obtained from the aggregation of local residuals of a forward discretization of the BSDE. In contrast with the backward formulation from [KT24], the forward discretization is more convenient for simulating systems of ergodic BDSEs. In our multidimensional setting, the terminal conditions of the random time horizon ergodic BSDE (2.2) are unknown for components i ̸= i0 . However, we address this issue by incorporating the constraint Yτi = Y0i in the loss function, while the normalization constraint Y i0 = y0 is included as a penalization term. The component Y is represented by a neural network, common to all time steps, while Z can either be computed with automatic differentiation or with another neural network. In the context of systems of ergodic BSDE, the ergodic cost λ will be estimated as a trainable parameter of the model, denoted by λ̄. Stochastic factor and random horizon approximations - The forward stochastic factor V v0 is approximated by a Euler discretization on a time grid π = (ti )i≥0 with t0 = 0 and time step h. Denoting for all i ≥ 0, ∆Wti = Wti+1 − Wti the Brownian increment at time ti : V ti+1

= V ti + µ(V ti )h + κ∆Wti ,

V 0 = v0 .

We denote by τ the first hitting time in the time grid of V to v0 after T0 . More precisely, assuming T0 ∈ π, we set (2.8)

τ = inf ti > T0 , ti ∈ π ; (V T0 − v0 )(V ti − v0 ) ≤ 0 . We denote Nτ = hτ the corresponding number of time steps and set Vτ = v0 .

Forward discretization of the eBSDE. Let Y θ1 : R → RI and Z θ2 : R → Rd×I be two neural networks, shared across all time steps, intended to approximate the Markovian components y and z of the system of ergodic BSDEs, associated with (2.2). The initial value of the forward discretization is taken as Y θ1 (v0 ). Without loss of generality, we can assume that i0 = 1. Given a realization of the Euler scheme V , we consider a forward discretization of Equation (2.2) on the time grid π, θ ,λ̄

θ ,λ̄

1 Y ti+1 = Y ti1 − F (V ti , Z θ2 (V ti ))h − G(Y ti1 )h + λ̄h + Z θ2 (V ti )⊤ ∆Wti ,

∀ 0 ≤ i ≤ Nτ .

(2.9)

The construction of local residuals relies on the iteration of the time discretization (2.9). In fact, for 1 ≤ i ≤ Nτ , θ ,λ̄

Y ti1

= Y θ1 (v0 ) + λ̄ti −

i−1 h X

i θ ,λ̄ F (V tk , Z θ2 (V tk )) + G(Y tk1 ) h − Z θ2 (V tk )⊤ ∆Wtk .

(2.10)

k=0

Loss function. The error at each time step is evaluated by computing the squared difference between the network prediction Y θ1 (V ti ) with the value (2.10) obtained by iteration of the forward discretization with the initial condition. For 1 ≤ i ≤ Nτ , define ϕti (V , θ1 , θ2 , λ) =

i−1 h X

i F (V tk , Z θ2 (V tk )) + G(Y θ1 (V tk )) h − Z θ2 (V tk )⊤ ∆Wtk − λ̄ti .

(2.11)

k=0

The quantity ϕti represents the cumulative increment between time 0 and time ti induced by the forward discretizations. The residual at time ti is then defined by 2

Errti = Y θ1 (V ti ) + ϕti (V , θ1 , θ2 , λ) − Y θ1 (v0 ) ,

(2.12)

At the initial time, the normalization condition is enforced through 2

(2.13)

Err0 = Y 1,θ1 (v0 ) − y0 .

At the terminal time tNτ = τ , the condition Vτ = v0 implies that the initial and terminal value cancels out in (2.12), leading to the residual 2

ErrNτ̄ = ϕtNτ (V , θ1 , θ2 , λ̄) .

These errors are then aggregated into a global loss function, which takes the form " Lloc (θ1 , θ2 , λ̄) = E

1,θ1

(v0 ) − y0

Nτ X

θ1

Y (V ti ) + ϕti (V , θ1 , θ2 , λ̄) − Y (v0 )

# (2.14)

i=1

In practice, this expectation is approximated by Monte Carlo over a set of trajectories V̄ j j=1,...B , up to their return time (τ j )j=1,...B , yielding the empirical loss

LB loc (θ1 , θ2 , λ̄) =

B 1 X





 Y 1,θ1 (v0 ) − y0

j=1

τ̄j X

j j Y θ1 (V ti ) + ϕti (V , θ1 , θ2 , λ̄) − Y θ1 (v0 )  (2.15)

i=1

Automatic differentiation variant. A variant consists of computing Z using automatic differentiation, rather than through an additional neural network. In that case, one sets Ztθi1 = κ

∂Y θ1 (v) ∈ Rd×I , ∂v v=V t

(2.16)

and optimizes only over (θ1 , λ). We refer to this modified algorithm as the ADLAeBSDE solver. Implementation. The complete procedure is summarized in Algorithm 1. At each gradient step, one first simulates a batch of size B of Euler trajectories then evaluates the empirical loss LB loc , and updates the parameters by stochastic gradient descent. Remark 2.2.

1. We also tested a global solver, in the spirit of [BQKMS24], based on the minimiza-

tion of a terminal loss at the random horizon. In the switching setting, this approach requires to introduce trainable values for the unknown starting points of forward discretization, for i ̸= i0 . Numerically, the method turned out to be unstable; the forward iteration generated large values in the loss, most likely due to the coupling term G. Clipping the differences inside g at CY did not improve stability, since it also suppressed the corresponding gradients. For this reason, we focused on the locally additive formulation. 2. Algorithm 1 could be extended to work for stochastic factor of dimension d′ > 1, by extending the definition of the random horizon τ to the return time in a centered ball around v0 . However, we expect the random horizon to have a greater mean, thus leading to high computational cost and possible time discretization error propagation.

2.3

Deep Galerkin method for ergodic PDE

Finally, we close this section by introducing an alternative residual-based neural-network algorithm for simulating systems of ergodic BSDEs, in the spirit of the Deep Galerkin Method (DGM) introduced in [SS18]. Since the DGM was introduced for high-dimensional parabolic PDE, a broad family of neuralnetwork solvers for PDEs has developed around the same idea, of representing the solution by a neural network and training it to satisfy the equation. The DGM and the closely related physics-informed neural networks of [RPK19] both minimize the residual of the PDE, and differ in the sampling used. Extensions of the DGM to the standard HJB equations is presented in [AACJ+ 22]. 16

Algorithm 1: Locally additive eBSDE Algorithm - (LAeBSDE) Let Y θ1 be a neural network defined on R, valued in RI with parameters θ1 and Z θ2 be a neural network defined on R, valued in Rd×I , with parameters θ2 . Let θ0 = (θ10 , θ20 , λ̄0 ) ∈ R3 be the initialization of the neural network and ergodic cost parameters. Define NT0 = ⌊ Th0 ⌋. for m = 0, ..., M do for j = 1, ..., B do j for k ∈ {0, ..., NT0 }, starting from V 0 = v0 do Sample ∆Wtjk ∼ N (0, hId ) and compute j

V tk+1 = V tk + µ(V tk )h + κ⊤ ∆Wtjk , Let Nj = NT0 + 1. j j while (V tN − v0 )(V tN − v0 ) > 0 do j

Sample ∆WtjN ∼ N (0, hId ) and compute j

V tN +1 = V tN + µ(V tN )h + κ⊤ ∆WtjN , j j j j Nj = Nj + 1 Set, hNj = τj . for j = 1, ..., B do Set, ϕt0 = 0. m j for k ∈ {1, ..., Nj }, starting from Y 0 = Y θ1 (v0 ) do m m m m j j j j ψtθk ,j = hF (V tk−1 , Z θ2 (V tk−1 )) + hG(Y θ1 (V tk−1 )) − λ̄m h − Z θ2 (V tk−1 )⊤ ∆Wtk−1 , m

ϕθtk ,j = ϕθtk−1,j + ψtθk ,j m

Compute Errθj (tk ) = Y θ1 (V tk ) + ϕθtk ,j − Y θ1 (v0 ) .   Nj B X X 2 m 1 m m m  Y 1,θ1m (v0 ) − y0 + Compute LB Errθj (tk ). loc (θ1 , θ2 , λ̄ ) = B j=1

k=1

m Update θm+1 = θm − ρm ∇θ LB loc (θ ).

In the following, we extend the application of the Deep-Galerkin method to the system of ergodic PDE (1.11), which we recall below Ly i (v) + F i (v, ∇y i (v)κ) +

q ij g(y j (v) − y i (v)) = λ,

′

∀ v ∈ Rd , ∀ i ∈ I.

j∈I

First, the equation is stationary and present an ergodic constant unknown λ, common to all equations, that is determined jointly with the solution (y i (.))i∈I . Second, the equation is posed on the whole space ′

Rd , so that there is no boundary condition and the choice of the sampling distribution for the residual becomes a modeling decision. While uniform sampling is common on bounded domain, we will prefer sampling from the invariant law of the underlying stochastic factor, in line with the unbounded ergodic setting. Instead, uniqueness is enforced through the normalization y i0 (v0 ) = y0 , in accordance with Theorem 1.2. A key feature of the DGM is that the dimension d′ of the state variable V v0 does not affect the formulation of the method: no spatial grid is required, and the PDE residual is evaluated ′

directly at randomly sampled points in Rd . The method approximates the unknown solution by a neural network and minimizes a loss function built from the PDE residual. This loss typically contains squared residuals evaluated at interior sampling points, together with a normalization term enforcing 17

y i0 (v0 ) = y0 . ′

Let Yθ be a neural network defined on Rd and valued in RI , with parameter θ. For each i ∈ I, define the residuals Rθi (v) = LYθi (v) + F i (v, ∇Yθi (v)κ) +

q ij g(Yθj (v) − Yθi (v)) − λθ .

(2.17)

j∈I

Ergodic cost approximation. The approximation of the ergodic cost λθ is based on its invariantmeasure representation (2.18). Replacing the exact solution by the neural network approximation, we define, for each regime i in I,  λiθ =

Z Rd′



F i (v, ∇Yθi (v)κ) +

q ij g(Yθj (v) − Yθi (v)) ν(dv).

(2.18)

j∈I

where ν denotes the invariant distribution of the stochastic factor V v0 . For the exact solution, the ergodic unknown λ is common to all states and thus does not depend on the regime i. During the numerical approximation, however, the quantities λiθ obtained from the network need not coincide. We therefore form a common approximation λθ by averaging over the regimes, λθ =

1X i λθ . I

(2.19)

An alternative approach would consist in treating λθ as an additional trainable scalar parameter, optimized jointly with the parameters of the neural network, as presented in the LAeBSDE and in [BQKMS24]. In this article, we will use the invariant measure representation 2.19 to approximate the ergodic cost in the DGM solver. In this setting, direct sampling from the invariant measure is possible, as considered in [GRS24]. In the general ergodic case, however, the invariant measure is typically not known explicitly and has to be approximated numerically. This can be done, for instance, using the recursive algorithm introduced in [LP02]. Loss function. We now define the loss minimized by the DGM algorithm. Since Theorem 1.2 requires the uniform bound (1.10), we add a penalization term in order to encourage the neural approximation to remain in this admissible region, namely Pθi (v) =

max( Yθi (v) − Yθj (v) − CY , 0).

i̸=j

The loss function is then defined by " L(θ) = Eν

Rθi (v)

# + Yθi0 (v0 ) − y0

i∈I

+ Eν

Pθi (v)

# ,

(2.20)

i∈I

where Eν denotes the expectation with respect to the invariant distribution. In practice, the expectations entering the loss are approximated by empirical averages based on a set of batchsize B. This 18

builds the DGM adaptation for ergodic BSDE, described in Algorithm 2. Remark 2.3. Although we choose to sample the collocation points from the invariant distribution ν, this choice is not mandatory for the residual term. One may alternatively sample points from another distribution, for instance uniformly on a compact region of interest, while still using invariant-measure sample for the approximation of λθ . Algorithm 2: Deep Galerkin method for eBSDE - (DGM eBSDE) ′

Let Y θ be a neural network defined on Rd , valued in RI with parameters θ. Let i0 ∈ I, ′ v0 ∈ Rd and y0 ∈ R be fixed, ensuring uniqueness of the Markovian solution to (1.1) such that y i0 (v0 ) = y0 . for m = 0, ..., M do R ⊤ ∞ Sample (V̄k )1≤k≤B from the invariant distribution N m, 0 e−µs κκ⊤ e−µ s ds for i = 1, ..., B do P P j i i ij i Compute λiθm = B k=1 (F (V̄k , ∇Yθm (V̄k )κ) + j∈I q g(Yθm (V̄k ) − Yθm (V̄k )). P Compute λθm = I1 Ii=1 λiθm for k = 1, ..., B do Set, ϕt−1 = 0. for i = 1, ..., I do P Rθi m (V¯k ) = LYθim (V¯k ) + F i (V¯k , ∇Yθim (V¯k )κ) + j∈I q ij g(Yθjm (V¯k ) − Yθim (V¯k )) − λθm , P Pθim (V̄k ) = i̸=j max( Yθi (v) − Yθj (v) − CY , 0). 1 Compute LB DGM (θm ) = B

B I X X

Rθi m (V¯k ) +

k=1 i=1 m Update θm+1 = θm − ρm ∇θ LB DGM (θ ).

I X

Pθim (V̄k ) + Yθi0m (v0 ) − y0

! .

i=1

Regime switching forward utilities in a stochastic factor model

In this section, we develop a general framework for decision-making in a regime-switching market using forward utilities, generalizing the results of [HLT20], where forward utilities are introduced in such a setting. We start by recalling the financial market model introduced in [HLT20]. Within this setting, we first derive a consistency stochastic partial differential equation that provides a broad characterization of regime-switching forward utilities. As a particular case, we recover the characterization of homothetic (power-type) forward utilities via systems of ergodic BSDEs obtained in [HLT20], and we derive systems of ergodic BSDEs associated with exponential and logarithmic forward utilities in regime switching markets. These ergodic BSDE systems enable the numerical approximation of the corresponding forward utilities and optimal strategies.

3.1

Regime switching stochastic factor model

We consider an agent who can invest her wealth in one riskless asset and n risky assets S = (S 1 , . . . , S n ), in a regime switching incomplete financial market. The regime switches are modeled by the continuous19

time Markov chain (CTMC) α on the state space I and transition rate matrix Q = (qij )i,j∈I , as defined in (1.14). In each regime market i ∈ I, the risky assets are characterized by a market price of risk (θi (Vtv0 ))t≥0 and a volatility matrix (σ i (Vtv0 ))t≥0 driven by the d′ -dimensional stochastic factor V v0 , as defined in (1.2). The bond is assumed to be the numeraire, and hence the stock price dynamics discounted by the interest rate is given by dSt = diag(St )σ αt− (Vtv0 )(θαt− (Vtv0 )dt + dWt ), Assumption 3.1.

(3.1)

′

1. For all v ∈ Rd and i ∈ I, the (n, d) matrix σ i (v) has full row rank n. ′

2. For all market regime i ∈ I, the market price of risk θi : Rd → Rn is uniformly bounded and Lipschitz continuous. The agent invests a proportion π̄ = π̄ 1 , ..., π̄ n

⊤

of her wealth X π in the n risky assets. For an

initial value X0π = x0 ∈ R+ , assuming the self-financing condition holds and rescaling the strategy vector by the volatility, the wealth process X evolves as dXtπ = Xtπ πt · (θαt (Vtv0 )dt + dWt ),

πt = σ αt− (Vtv0 )⊤ π̄t ∈ Rd .

(3.2)

Let Πi i∈I be a family of closed convex subsets of Rd , representing the constraint sets associated with each market regime. For any i ∈ I, we assume that π i = σ i (Vtv0 )⊤ π̄t i ∈ Ri , with Ri a vector space i

of Rd . In the following, for any a ∈ Rd , aR denotes the orthogonal projection on Ri and ai,⊥ its orthogonal projection onto Ri . For any t ≥ 0, the set of admissible strategies in [0, t] is defined as ( At =

πs =

πsi 1{αs− =i} , for s ∈ [0, t], ; ∀ i ∈ I, πsi ∈ Ri ,

i∈I

π is F prog measurable and i

Z t

2 πsi ds < ∞, P − a.s.

. (3.3)

For each t ≥ 0, the predictable strategy (πt )t≥0 is assumed to lie in Rαt− , which depends on the current S market regime. Finally, the set of admissible strategies is denoted by A = t≥0 At .

3.2

Forward utilities in a regime-switching environment

In this section, we give a general characterization of forward utilities in a regime-switching market. The introduction of regular random field spaces for the study of differentiability of Itô random fields is recalled in Annex A, see also [EKM13]. The decision criteria are built upon a family of random utilities U i i∈I , each modeling the agent’s preferences in regime i. For all i ∈ I, the utility random field U i : R+ × R+ × Ω → R associated with market regime i is an F-adapted random field such that: • The functions x ∈ R+ 7→ U i (t, x, ω) are nonnegative, strictly concave increasing functions of class C 2 on ]0, ∞[, (ω, t) a.s. • ui0 := U i (0, ·) is a standard (deterministic) utility function. 20

• Inada conditions: limx→0+ U i (t, x) = 0,

limx→0+ Uxi (t, x) = +∞,

limx→+∞ Uxi (t, x) = 0 a.s.

More specifically, we consider regime specific utilities U i (t, x) i∈I that are Itô random fields with 3,ϵ 3,ϵ local characteristics β i , γ i ∈ Kloc × Kloc with ϵ > 0, namely for i ∈ I (3.4)

dU i (t, x) = β i (t, x)dt + γ i (t, x) · dWt . ′

3,ϵ By Theorem 2.2 in [EKM13], U i is a Kloc semimartingale, for any ϵ′ < ϵ.

Regime switching utility The regime-switching utility U is the G-adapted càdlàg process defined as follows: X

U (t, x) = U αt (t, x) =

(3.5)

U i (t, x)1{αt =i} .

i∈I

At each regime-switching time Tk (i.e jump times of α), the random field U jumps from the utility U αTk−1 in market regime αTk−1 to U αTk in market regime αTk . The regime-switching utility U (3.5) have thus the following dynamics:

U (t, x) = U α0 (0, x) +

XZ t i∈I

Z t + 0





1{αs− =i} β i (s, x) +

q ij (U j (s, x) − U i (s, x)) ds

j∈I

1{αs− =i} γ i (s, x) · dWs +

Z tZ 0

U j (s, x) − U i (s, x) 1{αs− =i} Ñ (ds, di, dj), (3.6)

with Ñ the compensated Poisson measure as defined in Lemma 1.4. The decision criterion maintains time consistency within the given investment context, in the sense of the following definition, as introduced in [HLT20]. The optimal strategy provides maximal satisfaction to the agent, which is preserved at all times in the future. Definition 3.1 (Regime switching forward utility). A family of utility random fields U i (t, x) i∈I generates a regime switching forward utility if the G-adapted càdlàg process U defined by (3.5) satisfies the time consistency property: • For any admissible strategy π ∈ A, U (t, Xtπ ), t ≥ 0, is a locale supermartingale, i.e. U (t, Xtπ ) ≥ E[U (s, Xsπ )|Gt ]. ∗

• There exists an admissible optimal strategy π ∗ ∈ A such that U (t, Xtπ ), t ≥ 0, is a local martingale, i.e h i ∗ ∗ U (t, Xtπ ) = E U (s, Xsπ )|Gt . Time consistency In Theorem 3.1 we give a sufficient characterization of time consistency for regime-switching forward utilities. A general framework for forward utilities with jumps and the related nonlinear consistency SPDE has been introduced [MM22], extending the results of [MZ10], 21

[EKM13] in the continuous case, where the authors rely on a generalization of Itô-Ventzel’s formula with jumps introduced in [ØZ07]. Here, the regime-switching setting can be dealt with directly. Indeed, despite the jumps of the market price of risk θi , the wealth process X π is a continuous-time diffusion process. Hence, the dynamics of the compound process U (t, Xtπ ) is first obtained by applying the Itô-Ventzel’s formula with jumps. For U to satisfy the martingale consistency property from Definition 3.1, a necessary condition is that its drift β U (t, Xtπ ) ≤ 0, with equality holding along some optimal strategy. A natural candidate optimal strategy πt∗ (x) is thus the one that maximizes β U (t, ·), when seen as function of π. It is actually sufficient to impose that this maximum equals 0 to ensure that β U (t, Xtπ ) ≤ 0 for all admissible strategies π, thus providing time consistency of the regime-switching forward utility U . Finally, the following assumption ensures the admissibility of the optimal strategy. Assumption 3.2. Let U be a regime-switching utility as defined in (3.6). There exists a process K ∈ L2 (dt) such that a.s., for all i ∈ I, t ≥ 0, x > 0, i i,⊥ (t, x) , γxx (t, x) ≤ Kt Uxx

γxi,⊥ (t, x) ≤ Kt Uxi (t, x),

(3.7)

where γ i,⊥ = Proj(Ri )⊥ γ i . Theorem 3.1. Let U be a regime-switching utility as defined in (3.6), verifying Assumption 3.2. Assume that for each i ∈ I, the utility random field U i of local characteristics (β i , γ i ) verifies: X 1 i β i (t, x) = − Uxx (t, x) inf Qi (t, Xtπ , π) − U j (t, x) − U i (t, x) q ij , πt ∈Π 2

(3.8)

where

(3.9)

j∈I i U (t, x)θi (Vt ) + γxi (t, x) Qi (t, x, π) = ∥xπ∥2 + 2xπ · x . i (t, x) Uxx

Then, the regime switching utility U defined in (3.5) is a regime switching forward utility in the sense of Definition 3.1. The optimal investment strategy πt∗ is given by πt∗ =

πti,∗ 1{αt =i} ,

(3.10)

i∈I ∗

where πti,∗ := πti,∗ (Xtπ ) is the optimal strategy in regime i, with πti,∗ (x) =

−1 i (t, x) xUxx

ProjRi Uxi (t, x)θi (Vt ) + γxi (t, x) .

In particular, Qi (t, x, πti,∗ ) = − xπti,∗

(3.11)

Remark 3.1. Note that the expression for the optimal portfolio in regime i coincides with the nonswitching case, see e.g. [EKM13]. However, the dynamics of the regime-switching optimal strategy differ from those in the single regime setting since the consistency condition for each regime (3.8) induce a coupling of the dynamics of regime-specifc utilities . P Proof. Let (πt )t≥0 = ( i∈I πti 1{αt− =i} )t≥0 ∈ A be an admissible strategy, associated with the wealth X π . Under regularity assumptions on the local characteristics, we can apply the Itô-Ventzel’s formula 22

with jumps to U (t, Xtπ ). By (3.6), this yields that:  dU (t, Xtπ ) =



1{αt− =i} β i (t, Xtπ ) +

i∈I

X j∈I

1{αt− =i} γ (t, Xtπ ) · dWt +

i∈I

q ij (U j (t, Xtπ ) − U i (t, Xtπ ))dt

1{αt− =i}

Z I2

Uxi (t, Xtπ ) Xtπ πti ·

U j (t, Xtπ ) − U i (t, Xtπ ) 1{αt− =i} Ñ (dt, di, dj) i

θ (Vt )dt + dWt

i∈I

1 i π i 2 π i π i + Uxx (t, Xt ) Xt πt dt + Xt γx (t, Xt ) · πt dt . 2

Note that there is no jump quadratic variation because the wealth process (Xt )t≥0 is continuous. In integral form, the Brownian terms and compensated terms of the jump part are local martingales. Hence, a sufficient condition for consistency in this framework is that the drift of U (t, Xtπ ) is negative for any admissible strategy, and that there exists an optimal strategy π ∗ along which it vanishes. Accounting for the jumps’ contribution, the drift of U (t, Xtπ ) takes the form 

 β U (t, Xtπ ) =

X i∈I

1 i (t, Xtπ )Qi (t, Xtπ , πti ) + 1{αt =i} β i (t, Xtπ ) + Uxx 2

U j (t, x) − U i (t, x) q ij  ,

j∈I

(3.12) where Qi is given by (3.9). For i ∈ I, assume that the solution U i (t, .) obtained with (3.8) is concave. Then, its second derivative is negative, and thus each term of the sum in (3.12) is maximal when ∗

Qi is minimal. The candidate optimal policy in regime i is then defined by πti,∗ = πti,∗ (Xtπ ), with πti,∗ (x) := argminQi (t, x, π), and from the first order condition, π∈Ri

πti,∗ (x) = ProjRi

−1 i i i Ux (t, x)θ (Vt ) + γx (t, x) . i (t, x) xUxx ∗

∗

For all i ∈ I the quadratic form Qi attains its minimal value Q(t, Xtπ , π i,∗ ) = − Xtπ πti,∗

for this

strategy. The drift of U (t, Xtπ ) thus satisfies for any admissible strategy π ∈ A 



1 i 1{αt =i} β i (t, Xtπ ) + Uxx (t, Xtπ )Qi (t, Xtπ , πti ) + U j (t, x) − U i (t, x) q ij  2 i∈I j∈I   h i X X 2 1 i ∗ ∗ ∗ ≤ 1{αt =i} β i (t, Xtπ ) − Uxx (t, Xtπ ) Xtπ πti,∗ + U j (t, Xtπ ) − U i (t, Xtπ ) q ij  . 2

β U (t, Xtπ ) =

i∈I

j∈I

By the assumption (3.8), the right hand side of the inequality is equal to 0 and thus the process (U (t, Xtπ ))t≥0 is a supermartingale. It remains to show that for all i ∈ I, U i is an increasing and concave positive random field and the π ∗ ∈ A. This can be obtained as our framework falls under the framework of forward utilities with jumps introduced in [MM22], with the Lévy measure replaced by the jump measure of the CTMC α evolving on the finite state space I. Under the Inada condition and by Assumption 3.2, Theorem 3.8 in [MM22] is verified, which allows us to conclude. 23

Regime switching power forward utilities A typical choice for the regime specific utility random fields (U i )i∈I , are homothetic power type utilities, of separable form U i (t, x) = Pti

x δi , δi

δi ∈ (−∞, 0) ∪ (0, 1),

(3.13)

and with P i an Itô diffusion process with the following dynamics: dPti = Pti (bit dt + νti · dWt ). The following proposition is a corollary of Theorem 3.1, and shows that regime-switching power forward utilities have necessarily the same risk aversion in each regime. Proposition 3.2. Let (U i )i∈I be a family of power type utilities as above, defining a regime switching forward utility U by (3.5). Then, under the assumptions of Theorem 3.1, the regime specific utilities have necessarily the same risk aversion coefficient: δi = δ,

∀ i ∈ I,

and the consistency assumption (3.8) of Theorem can be written as: δi (δi − 1) ProjRi bit = 2

3.3

i 2 X j νt + θi (Vt ) Pt ij − q . 1 − δi Pti j∈I

Link between homothetic switching forward utilities and systems of ergodic BSDEs

A straightforward application of Proposition 3.2 and Itô’s formula allows us to establish a correspondence between regime-switching forward utility in power form and a system of ergodic BSDEs, thus recovering the characterization of [HLT20]. Corollary 3.3 ([HLT20]). Consider a family of utility random fields U i (t, x) i∈I of homothetic power type U i (t, x) =

xδ yi (Vt )−λt e , δ

δ ∈ (−∞, 0) ∪ (0, 1),

(3.14)

where (y i (.), z i (.))i∈I , λ is a Markovian solution of the system of ergodic BSDE (1.1), with generator F i and coupling term Gi given by 1 z + θi (v) δ ∥z∥2 2 δ(δ − 1) dist2 Ri , + z + θi (v) + , 2 1−δ 2(1 − δ) 2 X j i Gi (y) = ey −y − 1 q ij .

F i (v, z) =

j∈I

(3.15) (3.16)

This family generates a Markovian regime switching forward performance process U , and the associated optimal strategy in regime i takes the form πti,∗ = ProjRi

i z (Vt ) + θi (Vt ) . 1−δ

(3.17)

Similary, Theorem 3.1 also allows us to identify systems of ergodic BSDEs associated with regimeswitching forward utilities of exponential and logarithmic types. For switching exponential forward utilities, it is more convenient to use the discounted amount of wealth invested in the stock αt = Xtπ πt as a control variable, leading to the following wealth process dynamics: dXtα = αt⊤ (θαt− (Vt )dt + dWt ).

(3.18)

Proposition 3.4 (Exponential forward utility in regime-switching market). Consider a family of utility random fields U i (t, x) i∈I of homothetic exponential type i

U i (t, x) = −e−γx ey (Vt )−λt ,

γ∈R

(3.19)

where (y i (.), z i (.))i∈I , λ is a Markovian solution of the system of ergodic BSDE (1.1), with generator F i and coupling term Gi given by i ∥z∥2 γ2 1 2 2 i z + θ (v) F (v, z) = dist R , − z + θi (v) + , 2 γ 2 2 X j i ey −y − 1 q ij Gi (y) = i

(3.20) (3.21)

j∈I

This family generates a Markovian regime switching forward performance process U , and the associated optimal strategy in regime i takes the form αti,∗ = ProjRi

i z (Vt ) + θi (Vt ) . γ

(3.22)

Proof. The proof is a straightforward application of Itô’s formula and Theorem 3.1. The generator (3.20) coincides with the one associated to homothetic exponential forward utility from [LZ17]. We also recover the same exponential coupling term (3.21) as for regime switching forward performance process in power form. Finally, we state a similar result for utilities of logarithmic type. Proposition 3.5 (Logarithmic forward utility in regime-switching market). Consider a family of utility random fields U i (t, x) i∈I of homothetic logarithmic type U i (t, x) = ln(x) + y i (Vt ) − λt,

(3.23)

where (y i (.), z i (.))i∈I , λ is a Markovian solution of the system of ergodic BSDE (1.1), with generator

F i and coupling term Gi given by 1 1 2 F i (v) = − dist2 Ri , θi (v) + θi (v) , 2 2 X Gi (y) = y j − y i q ij .

(3.24) (3.25)

j∈I

This family generates a Markovian regime switching forward performance process U , and the associated optimal strategy in regime i takes the form αti,∗ = ProjRi θi (Vt ) .

(3.26)

The generator (3.24) also coincides with the one from [LZ17] in the single state setting. The coupling term (3.25) however is not exponential, because of the additive form of the logarithmic utility (3.23).

Numerical results

In this section, we present numerical experiments illustrating the performance of Algorithms 1 (LAeB SDE) and 2 (DGM) for the approximation of the Markovian solution (y i (.), z i (.))i∈I , λ of the system of ergodic BSDEs (2.2), with a particular emphasis on systems arising from regime-switching utilities of power type, see Corollary 3.3. Common numerical setting. Numerical experiments are conducted using Intel(R) Xeon(R) CPU @ 2.20GHz with 25GB of RAM. Unless otherwise specified, the stochastic factor V is an OrnsteinUhlenbeck (OU) process dVt = µ(m − Vt )dt + κ⊤ dWt , ′

′

V0 = v0 ,

(4.1)

′

with m ∈ Rd , µ ∈ Rd ×d and κ ∈ Rd ×d , which satisfies the dissipativity Assumption 1.1. The invariant measure is Gaussian, namely Z ∞ ⊤ N m, e−µs κκ⊤ e−µ s ds .

(4.2)

Points for the DGM and for validation error are sampled from this invariant distribution. Both solvers use feedforward neural networks with two hidden layers of 20 + Id neurons and hyperbolic tangent activations. We checked that increasing the depth and width does not improve the accuracy in our tests. The networks are trained by stochastic gradient descent using Adam optimize with initial learning rate 7 × 10−4 , batch size B = 100, and 10000 gradient descent steps. The algorithms are implemented in Python with Tensorflow library. The code of both solvers is available on github : https://github.com/gubrx/DeepSolvers-Ergodic-BSDE-Systems.

4.1

Some explicitly solvable examples

We first construct a family of systems of ergodic BSDEs admitting explicit Markovian solutions, which will serve as benchmarks for both algorithms. The construction starts from a family of ansatz (y i )i∈I satisfying the condition of uniqueness of Theorem 1.2, and a prescribed ergodic cost λ. We then define the market price of risk so that the ergodic PDE system (1.11) is satisfied. We consider the unconstrained power-type generator of Proposition 3.3, namely Π = Rd and F i (v, z) =

∥z∥2 δ 2 z + θi (v) + , 2(1 − δ) 2

Gi (y) =

j i ey −y − 1 q ij ,

(4.3)

j∈I

where δ ∈ (0, 1). Proposition 4.1. Let (y i (.))i∈I be a family of C 2 functions such that, for all i ∈ I: • y i has a bounded first derivative, • y i − y j is bounded uniformly in i, • Ly i is bounded. Set z i (v) = κ∇y i (v) ∈ Rd and let λ ∈ R be such that

1 i 2 λ > sup Ly (v) + G (y(v)) + z (v) . 2 i∈I, v i

(4.4)

Let 1 ∈ Rd denote the unit vector and define s i

θ (v) = −z (v) + 1

1 i 2(1 − δ) 2 i i λ − Ly (v) − G (y(v)) − ∥z (v)∥ . δ 2

(4.5)

Then (y i (.), z i (.))i∈I , λ is the unique Markovian solution of the systems of ergodic BSDEs (1.1) with generator and coupling term given by (4.3) and market price of risk given by (4.5). In the following numerical experiments, the factor process is the Ornstein-Uhlenbeck process (4.1) ′

with V valued in Rd and W is a d-dimensional Brownian motion. We consider benchmark solutions of the form y i (v) = ψ(v) + ci + ai ϕ(v),

ci , ai ∈ R

(4.6)

where ψ is sublinear and ϕ is bounded and smooth with bounded first and second derivatives. The differences y i − y j are then explicit bounded functions of ϕ and the coupling term is given by Gi (y(v)) =

(exp (cj − ci + (aj − ai )ϕ(v)) − 1)q ij .

j∈I ′

We fix a direction ℓ ∈ Rd and consider functions of the projection ℓ · v. ′

Example 4.2. Let w > 0, ℓ ∈ Rd . We consider the two profiles 27

(4.7)

(T) Hyperbolic tangent ϕT (v) = tanh(w ℓ · v). (R) Rational

ℓ·v . ϕR (v) = p 1 + w2 (ℓ · v)2

with w > 0 and ψ ≡ 0 in the sequel, so that Proposition 4.1 can be applied. Since the training losses of the two methods correspond to different objectives, they are not directly comparable. We therefore consider L2 validation errors computed under the invariant measure ν of the stochastic factor. Given approximations ȳ, z̄, and (v1 , . . . , vM ) an i.i.d sample from ν, we define M

EL2 (ν) (ȳ, M ) =

1 XX i 2 y (vm ) − ȳ i (vm ) , MI

EL2 (ν) (z̄, M ) =

m=1 i=1

1 XX i 2 z (vm ) − z̄ i (vm ) ,(4.8) MI m=1 i=1

with M = 1000 in all reported experiments. Hyperbolic tangent example (T) - We first consider a two-regime example (I = 2) with a onedimensional OU process (d′ = 1) initialized at v0 = 0, with mean-reversion coefficient µ = 2, asymptotic mean m = 0 and volatility coefficient κ = 0.65. The transition rates are q12 = 0.4, q21 = 0.8, and the reference solution is y 1 (v) = 1 − 0.3 tanh(0.8v),

y 2 (v) = 1 + 0.3 tanh(0.8v).

normalized by y 1 (0) = 1. The ergodic cost is set to λ = 0.811, which satisfies condition (4.4). For the LAeBSDE algorithm, the discretization parameters are h = 0.01 and the minimal horizon T0 = 1.

(b) Convergence of |λ̄ − λ| as a function of the number of epochs.

(a) Convergence of EL2 (ν) (ȳ, 1000) and EL2 (ν) (z̄, 1000) as a function of the number of epochs.

Figure 1. Error convergence for the hyperbolic tangent example (H).

The training behavior of both solvers is reported in Figure 1a. Both methods show a fast initial decrease of the L2 (ν)-errors, which stabilize around 10−2 after 4000 epochs. The LAeBSDE approximation provides the smallest final errors for both y and z. Figure 1b illustrates the error convergence 28

of the ergodic cost during training. For the DGM, the ergodic constant error decreases sharply and stabilizes around 10−3 once the error on y and z have themselves stabilized. This is consistent with the structure of the DGM solver, where λ̄ is approximated by the invariant-measure representation (2.18), whose accuracy is driven by that of (ȳ, z̄). By contrast, in the LAeBSDE, λ̄ is a trainable scalar optimized jointly with the networks; its trajectory displays spikes, inherited from the random horizon and from the time-discretization error that enters the loss at each gradient step. Ultimately, it reaches a higher accuracy, of order 10−4 . Figure 2 compares the approximated solutions ȳ i and z i produced by the two solvers against the exact solution, over the main support of the invariant law. Both approaches recover the qualitative structure of the solution, the two regimes crossing near v = 0 where the normalization y 1 (0) = 1 is enforced, while the gradients z 1 and z 2 are of opposite sign, small, and nearly flat over the relevant range. Overall the two methods provide approximations close to the exact solution on the major part of the invariant distribution and degrades only mildly towards the tails, where samples are sparce.

(a) Approximated solutions ȳ 1 , ȳ 2 against the exact solution.

(b) Approximated gradient component solution z̄ 1 , z̄ 2 against the exact solution.

Figure 2. Regime-wise solutions y i , z i as functions of the factor v, for the LAeBSDE and DGM solvers compared with the exact solution.

Time discretization parameter for LAeBSDE The sensitivity of the LAeBSDE scheme with respect to the minimal horizon T0 and to the time step h is reported in Tables 1 and 2. Increasing T0 from 0.1 to 2.0 improves all errors, particularly the error on the ergodic constant, which decreases from order 10−2 to 10−7 over the tested range. This behavior is coherent with the ergodic nature of the problem: increasing T0 lengthens the simulated trajectories, hence enriches the loss with information on the long-run behavior of the solution. From a numerical point of view, the term λ̄τ̄ appearing in the aggregation of error terms (2.11) acts as a long-time average over the whole trajectory, which makes the loss sensitive to a misspecification of λ̄. However, this effect saturates for large horizons, see the case T0 = 5 in Table 1, when the gain in long-run information is balanced by the accumulation of discretization errors. As shown in Table 2, refining the time step h improves the discretization and yields smaller errors on ȳ, z̄, and λ, but the gain is moderate. This comes at a significant computational cost: decreasing h increases the number of time steps before the random horizon, and therefore substantially increases the training time.

EL2 (ν) (ȳ, M )

EL2 (ν) (z̄, M )

|λ̄ − λ|

time (s)

0.10 0.25 0.50 1.00 2.00 5.00

2.25 × 10−2 1.74 × 10−2 1.10 × 10−2 2.96 × 10−3 4.68 × 10−3 5.02 × 10−2

1.33 × 10−2 1.03 × 10−2 9.97 × 10−3 8.91 × 10−3 8.42 × 10−3 4.22 × 10−3

1.25 × 10−3 7.33 × 10−3 3.11 × 10−3 1.19 × 10−8 1.19 × 10−8 1.19 × 10−8

453.9 460.7 502.5 542.7 648.6 873.5

Table 1. LAeBSDE errors and training times for different values of the minimal horizon T0 , with h = 0.01, in example (T).

EL2 (ν) (ȳ, M )

EL2 (ν) (z̄, M )

|λ̄ − λ|

time (s)

0.005 0.010 0.020 0.050 0.100

5.95 × 10−3 7.38 × 10−3 4.10 × 10−3 5.25 × 10−3 1.08 × 10−2

8.28 × 10−3 7.83 × 10−3 7.70 × 10−3 9.03 × 10−3 1.07 × 10−2

1.38 × 10−7 1.38 × 10−7 4.37 × 10−4 1.43 × 10−4 1.51 × 10−4

722.65 405.87 298.19 192.79 169.16

Table 2. LAeBSDE errors with respect to the time step h in example (T), with T0 = 1.

Switching rate sensitivity We next evaluate the behavior of the two algorithms as the number of regimes I increases. For each value of I, the transition matrix has diagonal entries −0.8 and constant off-diagonal entries q ij = 0.8/(I − 1), i ̸= j. The results are gathered in Table 3. Both methods remain accurate, the LAeBSDE yields the smallest errors on ȳ and z̄ for all values of I, and these errors stay essentially constant as the number of regimes grows, the error on ȳ remaining of order 10−3 from I = 2 to I = 20. The DGM, in contrast, presents a mild increase of its error on ȳ with the number of regimes, from 8.62 × 10−3 at I = 2 to 1.14 × 10−2 at I = 20. This is coherent with the deterioration of approximation quality of deep learning PDE with the dimension of the problem, as observed in [SS18]. The errors on z̄ remain of order 10−3 for both methods and do not display this trend. The ergodic constant is recovered up to numerical precision by both methods, with an error at the level of round-off in almost all cases. The main difference lies in the computational cost. The cost of the LAeBSDE grows rapidly with I since each gradient step requires the simulation of a batch of trajectories up to their random horizon and the evaluation of the I-dimensional coupling term along each trajectory. The DGM solves the stationary PDE system directly and avoids any time discretization, thus remaining 10 times faster compared to the LAeBSDE.

4.2

Regime switching power forward utility

We now focus on the systems of ergodic BSDEs associated with regime-switching forward utilities in power form. We consider a market alternating between two economic regimes: a growth regime (i = 1) with a high risk premium, and a more conservative regime (i = 2) with a lower one. The single risky asset is driven by a one-dimensional Brownian motion, as in (3.1), with regime-dependent volatility σ i > 0. When there is no constraint on the portfolio, that is Ri = Rd for all i ∈ I, the generator (3.15)

Method

EL2 (ν) (ȳ, M )

EL2 (ν) (z̄, M )

|λ̄ − λ⋆ |

time (s)

8.11 × 10−1

LAeBSDE DGM

4.44 × 10−3 8.62 × 10−3

7.98 × 10−3 1.03 × 10−2

5.53 × 10−5 1.19 × 10−8

420.9 30.1

6.03 × 10−1

LAeBSDE DGM

2.85 × 10−3 1.10 × 10−2

6.05 × 10−3 9.08 × 10−3

1.19 × 10−8 1.19 × 10−8

522.9 49.4

5.70 × 10−1

LAeBSDE DGM

3.82 × 10−3 1.11 × 10−2

5.95 × 10−3 8.10 × 10−3

1.19 × 10−8 1.19 × 10−8

1199.3 110.9

5.57 × 10−1

LAeBSDE DGM

4.52 × 10−3 1.14 × 10−2

5.65 × 10−3 8.18 × 10−3

1.19 × 10−8 1.19 × 10−8

3421.9 261.8

Table 3. Approximation errors across the number of regimes I (M = 105 Monte-Carlo samples).

reduces to F i (v, z) =

∥z∥2 δ 2 z + θi (v) + . 2(1 − δ) 2

(4.9)

We consider a truncated affine risk premium θi (v) = φb (ai + θi v) for i ∈ I, with ai ∈ R, θi > 0 and b > 0, where φb denotes the projection on the centered ball of radius b. Denoting θmax = max θi , i∈I

|δ| Assumption 1.2 holds with Cv = 1−δ max(1, b)θmax .

We assign a higher market price of risk and a lower stock’s volatility to the growth regime, namely a1 = 0.4,

a2 = −0.1

θ1 = 0.2,

θ2 = 0.05,

σ 1 = 0.15,

σ 2 = 0.3.

The remaining parameters are set to κ = 0.8, µ = 1.5, T0 = 1, v0 = m = 0, δ = 0.25, and b = 1, and the normalization is fixed to y 1 (0) = 1. The stochastic factor V then admits the Gaussian invariant law N (0, 0.213), whose 95% interval is approximately D = [−0.91, 0.91]. Over this range, the growth premium θ1 (v) ∈ [0.22, 0.58] remains positive, with mean 0.4, and the conservative premium θ2 (v) ∈ [−0.15, −0.06] remains negative with mean −0.1. The regime process is the continuous-time Markov chain α with asymmetric transition rate matrix Q= Its unique invariant distribution is ᾱ =

−0.3

0.3

1.0

−1.0

10 3 13 , 13

! (4.10)

, the mean holding times are approximately 3.33 in the

growth regime and 1 in the conservative one. This model that the market spends, on average, more than three times as long in expansion as in stagnation. Since no closed-form solution is available in this example, we assess the quality of the learned solutions through the residuals of the ergodic PDE (2.17). For each regime i, define Ri (ȳ, z̄, v) = Lȳ i (v) − F i (v, z̄ i (v)) +

I X

q ij g(ȳ j (v) − ȳ i (v)) − λ̄.

(4.11)

i=1

and the mean squared residual under the invariant measure, estimated by Monte Carlo over a sample 31

(v1 , . . . , vM ) from ν, M

EP DE (ȳ, z̄, M ) =

1 1 XX i 2 R (ȳ, z̄, vm ) . MI

(4.12)

m=1 i=1

We also report the normalization error Enorm = ȳ i0 (v0 ) − y0 . Figure 3 displays the training losses of both solvers, which decrease and stabilize around 10−2 for the LAeBSDE and 10−5 for the DGM, as well as the trajectories of the estimators λ̄, which converge to consistent values across the two methods.

(b) Convergence of λ̄ as a function of the number of epochs.

(a) Empirical loss functions for the LAeBSDE and DGM during training.

Figure 3. Convergence of empirical loss and λ̄ estimation for both the LAeBSDE and DGM.

This offers a good framework for the study of the associated regime-switching forward utility (3.5) generated by a family of utility random fields in power form as (3.14), which takes the form U (t, x) =

xδ yαt (v)−λt e , δ

y αt (v) =

y i (v)1{αt =i} .

(4.13)

i∈I

Simulating the jump times of the two Poisson processes N 1,2 and N 2,1 , and the stochastic factor with its Euler scheme, we are able to simulate the regime-switching forward utility for all time t ≥ 0. Figure 4a displays a realization of the regime-switching forward utility together with its non-switching counterparts, the vertical lines marking the jump times of the regime chain. The switching utility is higher than the regime frozen utility in regime 1, and lower in regime 2. This provide a trajectory of the regime switching forward utility that sits in between the trajectories of regime frozen preferences, this because of the coupling and the shared ergodic cost solution. Figure 4b displays the solution ȳ in each state as a function of v. The approximations produced by the two methods coincide on the high density region of the invariant law’s support, with discrepancies growing as v approaches 1, where samples are sparce. Impact of regime-switching frequency We investigate how the intensity of the regime switches affects the solution and its numerical approximation. We scale the transition rate matrix Q by a factor q ∈ {0.01, 0.1, 1, 10, 100}, that is Q(q) = qQ, which leaves the stationary distribution α unchanged 32

(b) Numerical approximation ȳ for LAeBSDE and DGM.

(a) Trajectory of the regime switching forward utility (4.13)

Figure 4. Trajectory of the regime-switching forward utility and solution

while accelerating the chain. The corresponding PDE residuals and normalization errors are reported in Table 4. We observe that the maximal distance between the two regime solutions max|y 1 − y 2 | v∈D

q (Q = q Q0 )

Method

EL2 (ν) (ȳ, M )

Enorm

λ̄

max|y 1 − y 2 |

0.01

LAeBSDE DGM

7.79 × 10−4 8.10 × 10−4

4.33 × 10−4 6.64 × 10−4

−1.15 × 10−2 −1.13 × 10−2

2.29 × 100 2.25 × 100

1.22 × 102

0.1

LAeBSDE DGM

7.92 × 10−4 8.15 × 10−4

2.34 × 10−4 4.69 × 10−4

−2.72 × 10−2 −2.66 × 10−2

3.13 × 10−1 3.09 × 10−1

1.22 × 101

LAeBSDE DGM

8.13 × 10−4 8.47 × 10−4

7.69 × 10−4 4.67 × 10−5

−2.94 × 10−2 −2.86 × 10−2

4.76 × 10−2 4.80 × 10−2

1.22 × 100

LAeBSDE DGM

1.24 × 10−3 8.59 × 10−4

3.99 × 10−2 8.71 × 10−4

−3.09 × 10−2 −2.91 × 10−2

7.63 × 10−3 5.55 × 10−3

1.22 × 10−1

100

LAeBSDE DGM

3.65 × 10−3 1.70 × 10−3

9.79 × 10−1 2.41 × 10−3

−2.88 × 10−2 −5.66 × 10−3

1.64 × 10−3 1.02 × 10−3

1.22 × 10−2

v∈D

Table 4. PDE residuals and normalization errors for different switching-rate scales q, and δ = −1 (M = 105 Monte-Carlo samples).

decreases by a factor 10 with q, which is coherent with the theoretical coupling bound CY , which scales as q −1 . In other words, when the chain switches rapidly, the regime specific solutions y 1 and y 2 collapse onto one common corrector, whereas for slow switching (q = 0.01), the two regimes are nearly decoupled and the solutions differ. The numerical error grows as q increases for both method, and the normalization error of the LAeBSDE deteriorates at high switching rates, since the switching term has a high weight. Risk aversion sensitivity Table 5 gathers the errors and estimated ergodic costs for several values of the risk-aversion parameter δ, from the risk-tolerant case δ = 0.5 to the strongly risk-averse case δ = −5. Several observations stand out. First, the residuals deteriorate as |δ/(1 − δ)| grows: the EPDE of both solvers increases by a factor 10 between δ = 0.25 and δ = 0.5. This is expected, as the factor δ 2(1−δ) governs both the size of the quadratic nonlinearity in (4.9) and the growth constant Cv , and

hence the magnitude of z. The two methods are comparable on the PDE residual, and both keep the normalization error below 10−3 . Additionally, the estimated ergodic cost λ̄ reproduces across all runs δ

Method

EPDE

Enorm

λ̄

time (s)

0.5

LAeBSDE DGM

5.32 × 10−3 5.61 × 10−3

1.30 × 10−4 3.08 × 10−5

8.02 × 10−2 7.89 × 10−2

331.8 8.5

0.25

LAeBSDE DGM

4.95 × 10−4 4.91 × 10−4

2.28 × 10−4 1.44 × 10−4

2.35 × 10−2 2.23 × 10−2

329.3 8.5

−1.0

LAeBSDE DGM

8.13 × 10−4 8.39 × 10−4

4.10 × 10−5 1.81 × 10−4

−2.93 × 10−2 −2.84 × 10−2

331.4 8.7

−2.0

LAeBSDE DGM

1.42 × 10−3 1.41 × 10−3

3.12 × 10−3 1.21 × 10−4

−3.84 × 10−2 −3.67 × 10−2

331.0 8.7

−5.0

LAeBSDE DGM

2.11 × 10−3 2.09 × 10−3

4.54 × 10−4 6.25 × 10−4

−4.59 × 10−2 −4.48 × 10−2

330.4 8.5

Table 5. PDE residuals EPDE , normalization errors and estimated ergodic costs for the LAeBSDE and DGM solvers, for several values of the risk-aversion parameter δ.

the sign relation sign λ̄ = sign δ. The two solvers agree on λ̄ to within 5% in every case, indicating a coherent approximation form both methods. The relation follows by evaluating the criterion at the admissible constant strategy π ≡ 0 between 0 and t, and exploiting the supermartingale property of the preference criterion U , which reads xδ0 −λt h yαt (Vt ) i xδ0 yα0 (V0 ) ≤ e E e e . δ δ Recall that y αt is bounded by sub-linearity of y i and Proposition 5 in [HL19]. Dividing by

(4.14) xδ0 δ

and

taking the logarithm in (4.14) yields α • if δ > 0, 1t ln E ey t (Vt ) − y α0 (V0 ) ≤ λ, so that letting t → ∞ leads 0 ≤ λ. α • if δ < 0, 1t ln E ey t (Vt ) − y α0 (V0 ) ≥ λ, so that letting t → ∞ leads 0 ≥ λ. Optimal strategies and comparison with the single-regime case. Recall that π i,∗ (v) = z̄ i (v)+θi (v) denotes the optimal strategy rescaled by the stock volatility, see (3.2). Figure 5 displays the 1−δ effective proportion of wealth invested in stock i, that is π̄ i,∗ (v) = σ1i π i,∗ (v), obtained from the DGM.

We report it both in the coupled (switching) case and in the decoupled case, where each regime is treated as a stand-alone market and the corresponding scalar ergodic BSDE of [LZ17] is solved with the same parameters. For δ = −1, the switching strategy is very close to the decoupled one and leads to realistic allocation levels. For δ = 0.5, the allocations are heavily leveraged and the discrepancy between the coupled and decoupled strategies becomes visible. This behavior can be explained from the structure of the solution. The regime coupling affects the δ in the generator strategy only through the z̄ i component. The size of z i depends on the weight 2(1−δ)

(4.9), which tends to 0 as δ → 0− and remains bounded by 1/2 as δ → −∞, whereas it blows up as δ → 1− . This weight also rescales the allocation size, leading to leveraged position for δ close to 1. 34

Consequently, for extreme risk averse agents (δ negative), the optimal allocation is dominated by the myopic component, and the anticipation of regime switches encoded in z i is minimal. For risk-tolerant agents (δ positive, close to 1), the component z is amplified, and the switching effect becomes visible, at the cost of extreme leverage. Figure 5 illustrate this observation.

(a) Optimal allocation for δ = −1.

(b) Optimal allocation for δ = 0.5

Figure 5. Optimal portfolio allocation for different risk aversion parameters.

Interpretation of the ergodic cost. For δ = 0.5, the ergodic cost is estimated at λ̄ = 8.02 × 10−2 . δ

In the representation U i (t, x) = xδ ey (Vt )−λt , the constant λ plays the role of an endogenous timepreference rate: a positive λ exponentially discounts future utility and thus reflects a preference for the present, whereas a negative λ up-weights future utility. The sign of this rate is tied to the agent’s risk preference: risk tolerant agents (δ > 0) discount future utility as a positive rate, while risk averse agents (δ < 0) exhibit a negative rate that favors future utility. Across the whole range of risk aversions considered, the ergodic cost remains of order 10−2 in absolute value, which is consistent with the magnitudes considered for long-term social discount rates in the policy literature. The magnitude of such rates is critical in long-term decision problems, and has been widely debated, in particular in environmental policy and sustainable development models, where excessively high discount rates are criticized for undervaluing the welfare of future generations; see [LH04], [Ste07]. A noticeable feature of forward utilities derived from systems of ergodic BSDEs is that this rate is not an exogenous input of the model but emerges as part of the solution (y, z, λ), jointly determined by the market dynamics, namely the risk premia and the factor’s ergodicity, and by the switching mechanism. This opens the possibility of calibrating the model parameters so as to reflect a desired time-preference structure.

Spaces for regularity of semimartingales

Let Ustd denote the set of standard deterministic utility functions. We next introduce spaces for studying the regularity of semimartingales in terms of their local characteristics, following [Kun97], [IW14]. Definition

seminorms

- Let

Rk -valued

random

field

class

C m,δ (]0, +∞[), with m a nonnegative integer and δ a number in (0, 1], i.e. β is m times differentiable in x and its mth derivative is δ-Hölder, for any t, almost surely. We introduce the following family of random Hölder K-seminorms aiming to control asymptotic behavior of β and the regularity of its Hölder derivatives, for any K ⊂]0, +∞[ X ∥β(t, x, ω)∥ sup ∂xj β(t, x, ω) + x x∈K x∈K

∥β∥m,K (t, ω) = sup ∥β∥m,δ,K (t, ω) =

1≤j≤m ∥∂xm β(t, x, ω) − ∂xm β(t, y, ω)∥ ∥β∥m,K (t, ω) + sup . |x − y|δ x,y∈K

As mentioned in [EKM13], the first term of these random semi-norms is divided by x, in order to control behavior in the neighborhood of x = 0, by imposing at most linear vanishing at the boundary. This normalization also preserves traditional asymptotic results from [Kun97]. Associated function spaces - The previous norms are related to the space parameter. We add the temporal dimension by requiring these seminorms (or their square) to be integrable in time with respect to Lebesgue measure on [0, T ]. We then define the following sets: m

m (resp. K ) denotes the set of C m -random fields β such that β and ∂ k β for k ≤ m are L1 1. Kloc loc x x

(resp. L2 )-locally bounded, that is for any compact K ⊂]0, +∞[ and any T , Z T 0

∥β∥m,K (t, ω)dt < ∞,

Z T 2 resp. ∥β∥m,K (t, ω)dt < ∞ . 0

m,δ

m,δ (resp. Kloc ) denotes the set of C m,δ -random fields such that for any compact K ⊂]0, +∞[ 2. Kloc

and any T , Z T 0

∥β∥m,δ,K (t, ω)dt < ∞,

Z T 2 resp. ∥β∥m,δ,K (t, ω)dt < ∞ . 0

3. When these norms are defined on the whole space ]0, +∞[, the derivatives up to a certain order are bounded in the spatial parameter, with an integrable (resp. square integrable) random bound, m

m,δ

so that we use the notation Kbm , Kb or Kbδ,m , Kb .

Proof of Lemma 1.1

For the sake of completeness, we provide the proof of the estimate of the running supremum in Lemma 1.1. Proof. Applying Itô’s formula to |Vt |2 , and using the dissipative Assumption 1.1 yields

d|Vt |2 = ≤

2Vt µ(Vt ) + |κ|2 dt + 2Vt κ⊤ dWt ! 2 |µ(0)| + |κ|2 dt + 2Vt κ⊤ dWt . −Cµ |Vt |2 + Cµ

Dropping the negative term and taking the supremum over s ∈ [0, τ ], Z t

sup |Vt |2 ≤ |v0 |2 + C0 τ + 2 sup 0≤t≤τ

0≤t≤τ

Vs κ⊤ dWs .

We now aim to obtain a bound of the above quantity in Lq . Elevating to the power q and taking expectation, we get that there exists a constant Cq such that Z t q 2q 2q q ⊤ ≤ Cq |v0 | + E[τ ] + E sup E sup |Vt | Vs κ dWs . 0≤t≤τ

0≤t≤τ

An application of the Burkholder-Davis-Gundy inequality in Lq , then allow to bound the local martingale, leading "Z q/2 #! τ 2q 2q q 2 q E sup |Vt | ≤ Cq |v0 | + E[τ ] + Cq |κ| E |Vs | ds 0≤t≤τ

≤ Cq

q 2q q q q |v0 | + E[τ ] + Cq |κ| E τ 2 sup |Vt | . 0≤t≤τ

(B.1)

Next applying Young’s inequality, for any ϵ > 0, there exists Cϵ such that q

τ 2 sup |Vt |q ≤ ϵ sup |Vt |2q + Cϵ τ q . 0≤t≤τ

0≤t≤τ

Injecting this back into (B.1) yields 2q 2q q 2q q q E sup |Vt | ≤ Cq |v0 | + E[τ ] + Cq |κ| (ϵE sup |Vt | + Cϵ E[τ ] . 0≤t≤τ

0≤t≤τ

2q For ϵ small enough, the term E sup |Vt | can be absorbed on the left hand side, leading 0≤t≤τ

2q E sup |Vt | ≤ Cq 1 + |v0 |2q + E[τ q ] . 0≤t≤τ

(B.2)

Proof of Proposition 3.2

Proof. Let i ∈ I and U i utility random field of power type verifying (3.13), so that U i have the local characteristics β i (t, x) = bit U i (t, x) and γ i (t, x) = νti U i (t, x). We start this proof by recalling useful relations for power utilities defined by (3.13): xUxi (t, x) = δ i U i (t, x), i (t, x) = δ i (δ i − 1)U i (t, x), x2 Uxx i (t, x) = (δ i − 1)U i (t, x). xUxx

Combining this with (3.11) yields that the optimal policy in regime i is given by: πti,∗ = ProjRi

i νt + θi (Vt ) . 1 − δi

Then, by Theorem 3.1, we have Q (t, x, πti,∗ ) = − i

xπti,∗

i i i νt + θ (Vt ) = −x dist R , , 1 − δi 2

and the time consistency assumption (3.8) becomes X ν i + θi (Vt ) 1 i (t, x)x2 Proj2 Ri , t bit U i (t, x) = Uxx − U j (t, x) − U i (t, x) q ij 2 1 − δi j∈I   X j i + θ i (V ) ν U (t, x) δ (δ − 1) t i i Proj2 Ri , t − − 1 q ij  . = U i (t, x)  2 1 − δi U i (t, x) j∈I

The time consistency assumption can thus be rewritten as: bit =

X j δi (δi − 1) ν i + θi (Vt ) Pt δi δi −δj ij Proj2 Ri , t q . − x 2 1 − δi Pti δj j∈I

By assumption 3.1, for j ̸= i q ij > 0. Thus, we have necessarily δi = δj for all i, j ∈ I, which concludes the proof.

References [AACJ+ 22]

Ali Al-Aradi, Adolfo Correia, Gabriel Jardim, Danilo de Freitas Naiff et Yuri Saporito : Extensions of the deep galerkin method. Applied Mathematics and Computation, 430:127287, 2022.

[ASS20]

Levon Avanesyan, Mykhaylo Shkolnikov et Ronnie Sircar :

Construction of a

class of forward performance processes in stochastic factor models, and an extension of Widder’s theorem. Finance and Stochastics, 24(4):981–1011, 2020. [BE02]

John Buffington et Robert J Elliott : American options with regime switching. International Journal of Theoretical and Applied Finance, 5(05):497–514, 2002.

[BP99]

Rainer Buckdahn et Shige Peng : Ergodic backward sde and associated pde. In Seminar on Stochastic Analysis, Random Fields and Applications: Centro Stefano Franscini, Ascona, September 1996, pages 73–85. Springer, 1999.

[BQKMS24]

Guillaume Broux-Quemerais, Sarah Kaakaï, Anis Matoussi et Wissal Sabbagh : Deep learning scheme for forward utilities using ergodic bsdes. Probability, Uncertainty and Quantitative Risk, 9(2):149–180, 2024.

[BR04]

Nicole Bauerle et Ulrich Rieder : Portfolio optimization with markov-modulated stock prices and interest rates. IEEE Transactions on Automatic Control, 49(3):442– 447, 2004.

[Cho19]

Wing Fung Chong : Pricing and hedging equity-linked life insurance contracts beyond the classical paradigm: The principle of equivalent forward preferences. Insurance: Mathematics and Economics, 88:93–107, 2019.

[CIL92]

Michael G. Crandall, Hitoshi Ishii et Pierre-Louis Lions : Users guide to viscosity solutions of second order partial differential equations. Bulletin of the American Mathematical Society, 27:1–67, 1992.

[CWNMW19] Quentin Chan-Wai-Nam, Joseph Mikael et Xavier Warin : Machine learning for semi linear PDEs. Journal of scientific computing, 79(3):1667–1712, 2019. [DHT11]

Arnaud Debussche, Ying Hu et Gianmario Tessitore : Ergodic BSDEs under weak dissipative assumptions. Stochastic Processes and their applications, 121(3):407–426, 2011.

[DRP21]

Goncalo Dos Reis et Vadim Platonov :

Forward utilities and mean-field games

under relative performance concerns. In From Particle Systems to Partial Differential Equations: International Conference, Particle Systems and PDEs VI, VII and VIII, 2017-2019 VIII, pages 227–251. Springer, 2021. [EKHM18]

Nicole El Karoui, Caroline Hillairet et Mohamed Mrad : Consistent utility of investment and consumption: a forward/backward spde viewpoint. Stochastics, 90(6): 927–954, 2018. 39

[EKHM22]

Nicole El Karoui, Caroline Hillairet et Mohamed Mrad :

Ramsey rule with

forward/backward utility for long-term yield curves modeling. Decisions in Economics and Finance, 45(1):375–414, 2022. [EKM13]

Nicole El Karoui et Mohamed Mrad : An exact connection between two solvable SDEs and a nonlinear utility stochastic PDE. SIAM Journal on Financial Mathematics, 4(1):697–736, 2013.

[FDR11]

Christoph Frei et Gonçalo Dos Reis : A financial market with interacting investors: does an equilibrium exist? Mathematics and financial economics, 4:161–182, 2011.

[FHT09]

Marco Fuhrman, Ying Hu et Gianmario Tessitore : Ergodic BSDEs and optimal ergodic control in Banach spaces. SIAM journal on control and optimization, 48(3): 1542–1566, 2009.

[Fre14]

Christoph Frei : Splitting multidimensional bsdes and finding local equilibria. Stochastic Processes and their Applications, 124(8):2654–2671, 2014.

[FWY14]

Jun Fu, Jiaqin Wei et Hailiang Yang : Portfolio optimization in a regime-switching market with derivatives. European Journal of Operational Research, 233(1):184–192, 2014.

[GM18]

Emmanuel Gobet et Mohamed Mrad : Convergence rate of strong approximations of compound random maps, application to SPDEs. Discrete & Continuous Dynamical Systems-Series B, 23(10), 2018.

[GPW+ 23]

Maximilien Germain, Huyên Pham, Xavier Warin et al. : Neural networks-based algorithms for stochastic control and PDEs in finance. Machine Learning and Data Sciences for Financial Markets, pages 426–452, 2023.

[GRS24]

Emmanuel Gobet, Adrien Richou et Lukasz Szpruch : Numerical approximation of ergodic bsdes using non linear feynman-kac formulas. arXiv preprint arXiv:2407.09034, 2024.

[GZ04]

X. Guo et Q. Zhang : Closed-form solutions for perpetual american put options with regime switching. SIAM Journal on Applied Mathematics, 64(6):2034–2049, 2004.

[HJ+ 17]

Jiequn Han, Arnulf Jentzen et al. : Deep learning-based numerical methods for highdimensional parabolic partial differential equations and backward stochastic differential equations. Communications in mathematics and statistics, 5(4):349–380, 2017.

[HKM24]

Caroline Hillairet, Sarah Kaakai et Mohamed Mrad : Time-consistent pension policy with minimum guarantee and sustainability constraint. Probability, Uncertainty and Quantitative Risk, pages 1–30, 2024.

[HL19]

Ying Hu et Florian Lemonnier : Ergodic BSDE with unbounded and multiplicative underlying diffusion and application to large time behaviour of viscosity solution of HJB equation. Stochastic Processes and their Applications, 129(10):4009–4050, 2019. 40

[HLT20]

Ying Hu, Gechun Liang et Shanjian Tang : Systems of ergodic BSDEs arising in regime switching forward performance processes. SIAM Journal on Control and Optimization, 58(4):2503–2534, 2020.

[HPW20]

Côme Huré, Huyên Pham et Xavier Warin :

Deep backward schemes for high-

dimensional nonlinear PDEs. Mathematics of Computation, 89(324):1547–1579, 2020. [HR19]

Jonathan Harter et Adrien Richou : A stability approach for solving multidimensional quadratic bsdes. Electronic Journal of Probability, 24:1 – 51, 2019.

[HT16]

Ying Hu et Shanjian Tang : Multi-dimensional backward stochastic differential equations of diagonally quadratic generators. Stochastic Processes and their Applications, 126(4):1066–1086, 2016.

[IW14]

Nobuyuki Ikeda et Shinzo Watanabe : Stochastic differential equations and diffusion processes. Elsevier, 2014.

[JR06]

Arnaud Jobert et Leonard CG Rogers :

Option pricing with markov-modulated

dynamics. SIAM Journal on Control and Optimization, 44(6):2063–2078, 2006. [KT24]

Lorenc Kapllani et Long Teng :

Deep learning algorithms for solving high-

dimensional nonlinear backward stochastic differential equations. Discrete and Continuous Dynamical Systems - B, 29(4):1695–1729, 2024. [Kun97]

Hiroshi Kunita :

Stochastic flows and stochastic differential equations, volume 24.

Cambridge university press, 1997. [LH04]

Franck Lecocq et Jean-Charles Hourcade : Le taux d’actualisation contre le principe de précaution? L’Actualité Economique, 80:41–65, 01 2004.

[LP02]

Damien Lamberton et Gilles Pagès : Recursive computation of the invariant distribution of a diffusion. Bernoulli, 8(3):367–405, avril 2002.

[LP17]

David A Levin et Yuval Peres :

Markov chains and mixing times, volume 107.

American Mathematical Soc., 2017. [LZ17]

Gechun Liang et Thaleia Zariphopoulou : Representation of homothetic forward performance processes in stochastic factor models via ergodic and infinite horizon bsde. SIAM Journal on Financial Mathematics, 8(1):344–372, 2017.

[LZ19]

Daniel Lacker et Thaleia Zariphopoulou : Mean field and n-agent games for optimal investment under relative performance criteria. Mathematical Finance, 29(4):1003–1038, 2019.

[MM22]

Anis Matoussi et Mohamed Mrad : Dynamic utility and related nonlinear spdes driven by levy noise. International Journal of Theoretical and Applied Finance, 25(01): 2250004, 2022.

[MZ06]

Marek Musiela et Thaleia Zariphopoulou : Investments and forward utilities. Technical report, 2006.

[MZ10]

Marek Musiela et Thaleia Zariphopoulou : Stochastic partial differential equations and portfolio choice. In Contemporary Quantitative Finance: Essays in Honour of Eckhard Platen, pages 195–216. Springer, 2010.

[NC24]

Kenneth Tsz Hin Ng et Wing Fung Chong : Optimal investment in defined contribution pension schemes with forward utility preferences. Insurance: Mathematics and Economics, 114:192–211, 2024.

[NZ14]

Sergey Nadtochiy et Thaleia Zariphopoulou : A class of homothetic forward investment performance processes with non-zero volatility. Inspired by Finance: The Musiela Festschrift, pages 475–504, 2014.

[ØZ07]

Bernt Øksendal et Tusheng Zhang : The itô-ventzell formula and forward stochastic differential equations driven by poisson random measures. Osaka Journal of Mathematics, 44:207–230, 2007.

[Par98]

Étienne Pardoux : Backward stochastic differential equations and viscosity solutions of systems of semilinear parabolic and elliptic PDEs of second order. In Stochastic Analysis and Related Topics VI: Proceedings of the Sixth Oslo—Silivri Workshop Geilo 1996, pages 79–127. Springer, 1998.

[PR14]

Etienne Pardoux et Aurel Rascanu : Stochastic Differential Equations, Backward SDEs, Partial Differential Equations. Springer, 2014.

[RPK19]

M. Raissi, P. Perdikaris et G. E. Karniadakis : Physics-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 378:686–707, 2019.

[SC09]

Luz Rocío Sotomayor et Abel Cadenillas :

Explicit solutions of consumption-

investment problems in financial markets with regime switching. Mathematical Finance: An International Journal of Mathematics, Statistics and Financial Economics, 19(2): 251–279, 2009. [SS18]

Justin Sirignano et Konstantinos Spiliopoulos : Dgm: A deep learning algorithm for solving partial differential equations. Journal of computational physics, 375:1339–1364, 2018.

[Ste07]

Nicholas Stern : The economics of climate change: the stern review. HM Treasury, 2007.

[Tev08]

Revaz Tevzadze :

Solvability of backward stochastic differential equations with

quadratic growth. Stochastic Processes and their Applications, 118(3):503–515, 2008. [Xv18]

Hao Xing et Gordan Žitković : A class of globally solvable markovian quadratic bsde systems and applications. The Annals of Probability, 46(1):491 – 550, 2018. 42

[YZZ06]

David D. Yao, Qing Zhang et Xun Yu Zhou : A Regime-Switching Model for European Options, pages 281–300. Springer US, Boston, MA, 2006.

[Zar92]

Thaleia Zariphopoulou : Investment-consumption models with transaction fees and markov-chain parameters. SIAM Journal on Control and Optimization, 30(3):613–636, 1992.

[ZY03]

Xun Yu Zhou et George Yin : Markowitz’s mean-variance portfolio selection with regime switching: A continuous-time model. SIAM Journal on Control and Optimization, 42(4):1466–1482, 2003.

Record · ID 303217 · SHA-256 99bd55e8efb44b76
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.