Approximation schemes for mixed optimal stopping and control problems with nonlinear expectations and jumps

We propose a class of numerical schemes for mixed optimal stopping and control of processes with infinite activity jumps and where the objective is evaluated by a nonlinear expectation. Exploiting an approximation by switching systems, piecewise constant policy timestepping reduces the problem to nonlocal semi-linear equations with different control parameters, uncoupled over individual time steps, which we solve by fully implicit monotone approximations to the controlled diffusion and the nonlocal term, and specifically the Lax-Friedrichs scheme for the nonlinearity in the gradient. We establish a comparison principle for the switching system and demonstrate the convergence of the schemes, which subsequently gives a constructive proof for the existence of a solution to the switching system. Numerical experiments are presented for a recursive utility maximization problem to demonstrate the effectiveness of the new schemes.


Introduction
Classical Markovian mixed optimal stopping and control problems, where the target is to maximise the (linear) expectation of a payoff on a finite time horizon T , are defined as u(t, x) = sup τ sup α E t,x τ t e −r(s−t) f (α s , X α,t,x s ) ds + e −r(τ −t) ξ(τ, X α,t,x τ ) , (1.1) Such nonlinear expectations arise naturally in financial mathematics, for instance as models for American options in a market with constrained portfolios [20], from recursive utility optimization problems [7], and for robust pricing and risk measures under probability model uncertainty [26]. It has been demonstrated in [15] that under suitable assumptions the value function u in (1.2) can be characterized by the viscosity solution to a more complicated HJBVI (1.3), which involves an extra nonlinearity resulting from the nonlinear expectation: with x = (t, x), nonlocal operators L α and B α , the driver f of the BSDE, and given functions ζ and g, which we will specify in Section 2. Particularly, in the case where the driver is additive in y and independent of z and k, i.e., f (α, x, y, z, k) ≡ f (α, x) − ry, the generalized control problem (1.2) reduces to the classical linear expectation case (1.1), and (1.3) reduces to an HJB obstacle problem. 1 As it is usually difficult to obtain analytic solutions of HJBVIs, it is necessary to design efficient and robust numerical methods for solving these fully nonlinear PIDEs.
We remark that, to the best of our knowledge, even for the case with linear expectations (i.e. f in the special form from above), there is no published numerical scheme covering the generality of (1.3). However, there is a vast literature on monotone approximations for local HJB equations (see, e.g., [10,2,11] and references therein) and a number of works covering specific extensions. For instance, monotone finite-difference quadrature schemes are proposed in [8,4,5,3] for nonlocal HJB equations. We refer the reader also to [12] for penalty approximations to nonlocal variational inequalities and to [28] for an application of policy iteration together with penalization to solve HJB obstacle problems. Probabilistic methods for solving HJB equations (without jumps and optimal stopping) can be found, for example, in [22].
However, this standard approach cannot be easily extended to nonlinear f which is only assumed to be Lipschitz (and generally is not semi-smooth [24]), which prevents a direct application of Newton-like solvers (see [29] for a special case of the discrete optimization problem with f (α, x, y, z, k) ≡ f (α, x, y) differentiable and concave in y, and A finite).
Moreover, at each step of policy iteration, one needs to identify the global optimal policy for each computational node. The nonlinear driver f and other PDE coefficients may have sufficiently complicate nonlinearities in the control variable such that the only way to construct a convergent algorithm is to discretize the admissible control set, and perform exhaustive search to determine the optimal policy at each node.
Another approach to solve (1.3) uses piecewise constant policy time stepping (PCPT) as in [25]. It is based explicitly on a discrete approximation of the admissible set by a finite set, say with J elements, and then defines a piecewise decoupled system of PDEs corresponding to these J (constant) controls. The information from the different solutions is assembled at the end of each timestep by taking the pointwise maximum.
As a specific scheme for these semi-linear PDEs, we propose an implicit Euler time discretization, monotone (semi-Lagrangian) approximations for local diffusions, a monotone quadraturebased scheme for the nonlocal terms, and the Lax-Friedrichs scheme for the nonlinearity in the gradient. The different solutions may be defined on different discretization grids by possibly high order monotonicity preserving interpolations. This approach not only avoids policy iteration, but also allows for an easier construction of convergent monotone schemes and an efficient parallel implementation of the individual semi-linear PDEs. Note that it is essential to obtain a monotone discretization, since it is well-known that non-monotone schemes may fail to converge or even converge to false "solutions" [11]. By Godunov's Theorem [18], in general, one can expect a monotone scheme to be at most first-order accurate.
The main contributions of our paper are: • We formulate our algorithm by approximating the solution of (1.3) by the solution to a switching system with small switching cost. We shall establish a comparison principle for the switching system and demonstrate that as the switching cost tends to zero, the solution of the switching system converges to the viscosity solution of (1.3), which extends the results in [4] to obstacle problems of switching systems and includes nonlinear drivers.
• We discretize the switching system piecewise in time by fully implicit monotone approximations. The convergence of the scheme is demonstrated, which subsequently gives a constructive proof for the existence of a viscosity solution to the switching system. Our results extend the one obtained in [25] from the case of standard control problems. In contrast to there, PCPT leads to coupled semi-linear PDEs rather than linear PDEs due to the nonlinear expectations. The optimal stopping right is treated as an additional control and included in the switching directly instead of the classical penalisation approach.
• By truncation of the singular jump measure, we obtain a stochastic control and optimal stopping problem whose value function is shown to converge to the value function of the initial problem and which satisfies a HJBVI equation. Our result extends earlier ones obtained only in the case of a linear expectation without control and optimal stopping (see e.g. [8]).
• For practical implementations, we propose a Picard-type iteration for the efficient numerical solution without the need to invert the dense matrices resulting from the nonlocal terms.
• Numerical examples for a recursive utility maximization problem are included to investigate the convergence order of the scheme with respect to different discretization parameters.
The remainder of this paper is organized as follows. In Section 2, we introduce the Markovian mixed optimal stopping and control problem with nonlinear expectations, and characterize its value function as the viscosity solution of a nonlocal HJBVI. We then derive numerical schemes in Section 3 by approximating the HJBVI with a switching system, PCPT, and ultimately fully discrete monotone schemes. Then we move on to the convergence analysis of our numerical schemes in Section 4. Numerical examples for a recursive utility maximization problem are presented in Section 5 to illustrate the effectiveness of our algorithms. In the Appendix, we include a rigorous proof of the comparison principle for the switching system and some complementary results that are used in this article.

Problem formulation and preliminaries
In this section, we formulate the mixed optimal stopping and control problem with nonlinear expectation and introduce the connection between such problems and HJBVIs, which is crucial for the subsequent developments. We start with some useful notation that is needed frequently in the rest of this work.
We write by T > 0 the terminal time, and by (Ω, F, P ) a complete probability space, in which two mutually independent processes, a d-dimensional Brownian motion W and a Poisson random measure N (dt, de) with compensator ν(de)dt, are defined. We assume ν is a σ-finite measure on E := R n \ {0} equipped with its Borel field B(E) and satisfies We denote by E the usual expectation operator with respect to the measure P . For any given t ∈ [0, T ], we define the t-translated Brownian motion W t := (W s − W t ) s≥t and the t-translated Poisson random measure N t := N (]t, s], ·) s≥t . We denote byÑ t (dt, de) = N t (dt, de) − ν(de)dt the compensated process of N t , and by F t = {F t s } s∈[t,T ] be the filtration generated by W t and N t augmented by the P -null sets.
Furthermore, we introduce several spaces: L 2 ν is the space of Borel functions l : We now proceed to introduce the control problem of interest. For each t ∈ [0, T ], let A t t be a set of admissible controls, which are F t -predictable processes (α s ) s∈[t,T ] valued in a compact set A, and T t t be the set of F t -stopping times which take values in [t, T ]. For any given initial state x ∈ R d , and control α ∈ A t t , we consider the controlled jump-diffusion process (X α,t,x s ) t≤s≤T satisfying the following SDE: for each s ∈ [t, T ], where b, η ∈ R d and σ ∈ R d×d are given measurable functions. We remark that although our analyses are performed only for jump-diffusion processes with time-homogenous coefficients, similar results are valid for controlled dynamics with time-dependent coefficients. The performance of the control problem, depending on α, is evaluated by a nonlinear expectation induced by a BSDE with a controlled driver f (α s , s, X α,t,x s , y, z, k). That is, for any given stopping time τ ∈ T t t and any bounded Borel function ξ, we define the nonlinear expectation where the process (Y α,τ,t,x s ) s≤τ is a solution in S 2 t of the following BSDE: for each s ∈ [t, τ ], and (Z α,t,x s,τ ), (K α,t,x s,τ ) are two associated processes, if they exist, lying in H t and H ν t , respectively. Now we are ready to state the generalized mixed optimal stopping and control problem. For each initial time t ∈ [0, T ] and initial state x ∈ R d , we consider the following value function: (2.4) subject to the controlled SDE (2.2), where ξ is the terminal position given by for some reward functions ζ and g. Note that the value function of our control problem is constant up to a P -null set. Throughout this work, we shall perform the analysis under the following standard assumptions on the coefficients: Assumption 1. The set of control values A is compact and the driver f is a measurable function of the form f (α, s, x, y, z, k) :=f (α, s, x, y, z, E k(e)γ(x, e) ν(de))1 s≥t for some functionsf and γ. Moreover, there exists a constant C > 0 such that for any α, is continuous in t and admits the properties: Assumption 1 is the same as that made in [16]. Note that under assumptions (1), (2), equation (2.2) admits a unique solution. Assumptions (3), (4.a), (4.c), (5) quarantee the existence and the uniqueness of the solution of equation (2.3). Consequently, the generalized mixed control problem is well-defined. For notational convenience, in the sequel, we will writef as f , and denote by ψ α a generic function ψ with control-dependence.
The rest of this section is devoted to the equivalence between the mixed control problem and a generalized nonlocal HJBVI. Specifically, we now consider a Hamilton-Jacobi-Bellman variational inequality of the following form: x) contains both the time t and the spatial coordinate x ∈ R d , and the nonlocal operators L α := A α + K α and B α satisfy, for φ ∈ C 1,2 (Q T ): where E = R n \ {0} is defined at the beginning of Section 2 and the nonlocal operators K α u and B α u are well-defined under Assumption 1.
We emphasize that since the matrix σ α (σ α ) T is only assumed to be nonnegative definite, both the diffusion coefficient σ α (σ α ) T and the jump intensity η of (2.5) are allowed to vanish at some points. Consequently, there is no Laplacian smoothing from the second-order differential operator nor fractional Laplacian smoothing from the nonlocal operator to this degenerate equation (2.5). Therefore, in general, this HJBVI will not admit classical solutions, and we shall interpret the equation in the following viscosity sense based on semi-continuous envelopes of the equation [1,25].
Definition 2.1 (Viscosity solution of HJBVI). An upper (resp. lower) semicontinuous function u is said to be a viscosity subsolution (resp. supersolution) of (2.5) if and only if for any point x 0 and for any φ ∈ C 1,2 (Q T ) such that φ(x 0 ) = u(x 0 ) and u − φ attains its global maximum (resp. minimum) at x 0 , one has A continuous function is a viscosity solution of the HJBVI (2.5) if it is both a a viscosity suband supersolution.
Under Assumption 1, the HJBVI (2.5) is well-posed in the class of bounded continuous functions (see [14,15]). The unique viscosity solution of (2.5) (after a change of time variable) can be further characterized as the optimal value function (2.4) of the mixed control problem. In other words, to obtain the optimal value function for all initial times t and initial states x, it is equivalent to design effective numerical schemes to solve (2.5).
Moreover, a strong comparison principle holds for the HJBVI (2.5), the proof of which is similar to that in [14] (without controls) and hence omitted. In particular, if U is a bounded viscosity subsolution and V is a bounded viscosity supersolution to (2.5)

Construction of numerical schemes
In this section, we will design numerical schemes for solving HJBVI (2.5). We carry out the following string of approximations to construct our numerical algorithm: • truncation of the singular jump measure (equation (3.3) and Appendix C); • approximation of the control set with a finite set (equation (3.4) and Theorem 4.2); • approximation of the discretized control problem with a switching system (equation (3.5) and Theorem 4.5); • discretization in time and space (equations (3.7, 3.8) and Theorem 4.6).
We start the derivations of our schemes by approximating the singular measure ν with a truncated non-singular measure and a modified diffusion coefficient as suggested in [8]. This can be done by introducing an approximative jump-diffusion dynamics and an approximative backward SDE (see Appendix C). More precisely, for any given ε > 0, let us define the truncated measure ν ε (de) = 1 |e|>ε ν(de) and the modified diffusion coefficientσ α (x) such thatσ α We further introduce the modified local operator A α ε as: (3.2) and truncated nonlocal operators K α ε and B α ε by replacing ν with ν ε in (2.7) and (2.8), respectively. With these operators in hand, we consider the following modified HJBVI: . These modified coefficients clearly satisfy Assumption 1, and hence (3.3) is well-posed in the viscosity sense. Remark 1. In Appendix C, we provide an alternative interpretation of the above approximation by identifying the viscosity solution of (3.3) as the value function of a mixed control problem in terms of modified SDE and BSDE. This characterization further enables us to establish the convergence of this approximation through a probabilistic argument.
We then approximate the admissible control set in (3.3) by a finite set. More precisely, for a finite subset A δ of the compact set A such that max α∈A miñ α∈A δ |α −α| < δ, we introduce the finite control HJBVI by where we denote for simplicity the nonlinear function f α ε [u] := f α (x, u δ , (σ α ) T Du δ , B α ε u δ ). Since (3.4) is a special case of (2.5) with a finite admissible set, it is clear that (3.4) admits a unique bounded viscosity solution.
Next, we approximate the finite control equation (3.4) by a switching system ( [2,4]). Suppose the finite control set is given by A δ = {α 1 , α 2 , . . . , α J }. We denote by U ε,δ,c j , j = 1, . . . , J the solution of the following system of HJB equations: where we define M j U := max k =j U k − c for any c > 0. The positive switching cost is needed for the well-posedness of the switching system (3.5).
We now proceed to introduce a discrete approximation to the switching system based on the idea of piecewise constant policy timestepping. Define a set of nodes {x j,i } and timesteps t n , with discretization parameters h and ∆t, i.e., By parameterizing the grid Ω j,h = {x j,i } i with the control index j, we are allowing the usage of different discretization grids for different controls. We denote by U n j,i the discrete approximation to U j at the point x n j,i = (t n , x j,i ), and extend it to the computational domain by interpolation. Let L α j ε,h and f α j ε,h be the discrete form of the operators L α j ε and f α j ε , respectively. We discretize (3.5) on the grid Ω j,h with the uniform time partition t n+1 − t n = ∆t by performing piecewise constant policy timestepping and applying the constraints at the beginning of a new timestep, whereŨ n k,i(j) is the value of the interpolant of {U n k,l } l∈Ω k,l at the i-th point of the grid Ω j,h . Now by rearranging the terms of (3.8), we obtain the following numerical scheme: for x n j,i ∈ Q T and j = 1, . . . , J, As seen from (3.9), performing switching at the beginning of a new timestep introduces two additional terms to both the switching part and the obstacle part of the equation, which will not appear in a straightforward discretization of the switching system (3.5). However, we will demonstrate in Section 4.4 that these terms vanish as ∆t, h → 0, and consequently our scheme (3.9) forms a consistent approximation to (3.5).
For notational simplicity, we label our approximations only by h and assume in the sequel that ∆t is a given function of h with ∆t → 0 as h → 0.

Remark 2.
To determine the optimal stopping strategy and the optimal controls, one can simply compare values of the obstacle and all components of the switching system at each grid point. As we will see in Section 4.2, each component of the switching system converges to the solution of the HJBVI (2.5) as the discretization parameters tend to 0; therefore, although we have no guarantee for the convergence of the control approximation, this numerically found control is close to optimal for the mixed control problem (2.4).
We now describe in detail how we perform spatial discretizations for L α ε and B α ε to construct a monotone discrete operator L α ε,h and f α ε,h , for a fixed control parameter α ∈ {α 1 , . . . , α J }. To simplify the presentation, we consider a piecewise linear or multilinear interpolation I h on a uniform spatial grid hZ d . That is, We start with the nonlocal terms. The definition of K α ε gives Then, by replacing the integrands by their monotone interpolants (c.f. [5]), we derive the following approximations for K α,1 ε and B α ε (where we have dropped the mesh index j in x for simplicity): with the coefficients which are well-defined and nonnegative, and consequently result in monotone approximations. These coefficients can be efficiently evaluated by using quadrature rules with positive weights, such as Gauss methods of appropriate order [5].
With these modified coefficients, we are ready to construct the following approximations of the local operators: for any k > 0, which, by using (3.10) and the fact that m ω m = 1, can be further written in the discrete monotone form with non-negative coefficients The approximation to the local operator A α ε falls into the class of semi-Lagrangian schemes (see e.g. [11]), and provides a consistent monotone approximation for possibly degenerate, non-diagonally dominant diffusion coefficients.
Before presenting our fully discrete scheme, we shall point out that by considering a truncated problem, one can without loss of generality assume that σ α is bounded, which consequently implies the Hamiltonianf is Lipschitz continuous with respect to p. Indeed, suppose σ α is unbounded, then for any given µ > 0, we define the cut-off function and consider a truncated HJBVI (3.4) by replacing σ α with the bounded diffusion coefficient . Using the fact that 1 − ξ µ → 0 uniformly on compact sets as µ → 0, one can easily prove this additional approximation is consistent with (3.4), and hence its viscosity solution converges to the solution to (3.4) uniformly on bounded sets.
The Lipschitz continuous Hamiltonian enables an approximation by the implicit Lax-Friedrichs numerical flux [10]. For each l = 1, . . . , d, we denote by ∆ (l) − U n j,i ) the one-step forward (resp. backward) difference operator along the l-th coordinate, and by ∆U n j,i = (∆ − U n j,i ) T the central difference operator at the grid point x n j,i . Then for any θ > 0, the Lax-Friedrichs numerical flux is given for any ( where we define λ = ∆t/h. A fully implicit time discretisation is finally given by Substituting (3.7) into (3.17), one can reformulate this implicit scheme in its equivalent form (3.9). We end this section with a remark about the implementation of the implicit scheme (3.17). To avoid solving linear systems with the dense matrices resulting from the discretization of the nonlocal operators, we write the solution to (3.17) as the fixed point of a sparse contraction mapping T , such that sufficient accuracy is achieved in practice by a few fixed point iterations.
Given bounded functions U n− 1 2 j and U n,(k) j , we define the following mapping T on ℓ ∞ (Z d ), i.e., the Banach space of bounded functions on hZ d employed with the sup-norm | · | 0 : It is clear that a fixed point U n j with T U n j = U n j is a solution to (3.17). Moreover, for any given functions U n,(k) j and V n,(k) j in ℓ ∞ (Z d ), we obtain from the Lipschitz continuity of f and the ℓ ∞ stability of the numerical fluxf (see Lemma 4.10) that Since we need h = o(ε) in general to achieve consistency of our scheme (see Lemma 4.8), it suffices to require ∆t ε 2 < 1 and 4dθ < 1 to ensure T is a contraction mapping on ℓ ∞ (Z d ). This establishes the well-posedness of (3.17) and enables us to solve the nonlinear equation (3.17) through Picard iterations by setting We emphasize that the criterion ∆t ε 2 < 1 is a sufficient condition in the worst case, but is often far from computationally optimal since we have used no information about the exact behavior of the singular measure ν around zero (see Remark 4 for details). For typical Lévy measures from finance [8], we only need ∆t = O(ε) (such as in the variance gamma case in our tests) or even ∆t independent of ε (for instance, for a Gaussian density) to guarantee T is a contraction mapping.
Moreover, since in practice we can evaluate the discrete nonlocal operators K α,1 ε,h U n j and B α ε,h U n j at all grid points in O(N log N ) operations using a FFT (see e.g. [12]), and in each iteration a sparse linear system is solved (in one dimension it is tridiagonal), the total complexity of the implicit scheme is still close to linear. It is also worth pointing out that even if we chose explicit approximations for the nonlocal operators, due to the nonlinearity of f , one still has to perform iterations to solve the resulting nonlinear equations. Because we only assume Lipschitz continuity but no higher regularity of f , we adopt the Picard iteration to solve for U n .

Convergence analysis
In this section, we establish the convergence of the numerical approximations in Section 3 to the viscosity solution of the HJBVI (2.5).
We start by outlining the convergence analysis of the truncation of singular measures. It is not difficult to see that (3.3) is a consistent approximation of (2.5) in the viscosity sense, such that the comparison principle of (2.5) enables us to conclude that the solution of (3.3) converges to that of (2.5) on compact sets as ε → 0. In Appendix C, we provide an alternative proof by identifying the viscosity solution of (3.3) as the value function of a modified control problem, and establish its convergence to the original value function (2.4) using a probabilistic argument.
The remainder of this section thus focuses on the convergence analysis for the control discretization, the approximation of HJBVIs with switching systems, and the discrete approximations of the switching systems.
To simplify the notation, we will occasionally drop the terms {K α u} α∈A and {B α u} α∈A in (3.3), (3.4) and (3.5), and simply denote hem by F (x, u, Du, D 2 u) = 0 in the sequel.

Approximation by finite control sets
In this section, we shall study approximations of HJBVI (2.5) with a finite control set. The following consistency result will be essential for our convergence analysis. The proof is similar to [25] (who consider the case without nonlocal term, obstacle, and nonlinear driver) and included in Appendix A for the convenience of the reader.
Now we are ready to conclude the convergence of our approximation with finite control sets.
The proof is a straightforward extension of the arguments in [25] and is hence omitted.

Approximation by switching systems
In this section, we study the approximation of (3.4) by switching systems. We adopt the following standard definition of a viscosity solution to switching systems of the form (3.5) (see [1,4,25] and references therein). Definition 4.3 (Viscosity solution of switching system). A R J -valued upper (resp. lower) semicontinuous function U is said to be a viscosity subsolution (resp. supersolution) of (3.5) if and only if for any point x 0 and for any φ ∈ C 1,2 (Q T ) such that U j − φ attains its global maximum (resp. minimum) at x 0 , one has A continuous function is a viscosity solution of the HJBVI (3.5) if it is both a a viscosity suband supersolution.
Note that in the definition of the viscosity solution of F j , the test function only replaces U j in the integrals and derivatives, while leaving the terms {U k } k =j unchanged. Now we present the comparison principle for bounded semicontinuous viscosity solutions of (3.5), which not only implies the uniqueness of bounded viscosity solutions of (3.5), but is also essential for our convergence analysis. The proof will be given in Appendix B. The following theorem demonstrates the convergence of the switching system to the finite control HJBVI (3.4) as the switching cost goes to 0. Convergence with order 1/3 is proved in [4] by a different technique, for nonlocal Bellman equations without obstacles and nonlinear source terms.
We momentarily assume the switching system (3.5) to admit a viscosity solution bounded independently of the (small enough) switching cost c. We give a constructive proof of existence through our numerical schemes in Section 4.3. ) and u ε,δ be the viscosity solution of (3.5) and (3.4), respectively. Then for fixed ε, δ > 0, we have for each j = 1, . . . , J that U ε,δ,c j → u ε,δ uniformly on compact sets as c → 0.
Proof. Since ε and δ are fixed for our analysis, we shall omit the dependence on ε and δ, and simply denote by U c the solution of (3.5). Consider a sequence of switching costs c m → 0 as m → ∞, and the corresponding viscosity solution U cm = (U cm 1 , . . . , U cm J ). We shall first prove by contradiction that U cm j (x) ≥ M j U cm , x ∈ Q T , j = 1, . . . , J.
Suppose the statement is false, then there would exist k = j and x 0 ∈ Q T such that U cm j (x 0 ) < U cm k (x 0 ) − c m . We then obtain from the continuity of U cm j and U cm k that there exists a nonempty open ball B around x 0 such that On the other hand, there exists a C 2 function φ such that U cm j − φ attains its minimum at some point in B, say x 1 . Hence we deduce from the fact that U cm j is a supersolution that which leads to a contradiction. We now introduce the following functions through a relaxed limit: for j = 1, . . . , J, Letting r → ∞ leads to the fact that U j ≥ U k for all j = k. The statement for {U j } can be shown similarly.
Since it is clear that U and U is bounded upper and lower semicontinuous, respectively, we now aim to show U and U is respectively a sub-and supersolution of (3.4). Then the strong comparision principle gives us U ≤ U , which implies U = U = U is the unique viscosity solution of (3.4). Uniform convergence on compact sets follows from a variation of Dini's theorem (See Remark 6.4 in [9]).
We start by showing U is a subsolution of (3.4). Let φ ∈ C 1,2 and U − φ have a strict global maximum atx 0 ∈Q T , then there will be a sequence c m → 0 such that for each j ∈ {1, . . . , J}, we havex j m →x 0 , U cm j (x j m ) → U (x 0 ), and U cm j − φ attains a global maximum atx j m . Since U cm j is a subsolution of (3.5) with c m , if we havex 0 ∈ {0} × R d , U cm j (x j m ) ≤ g(x j m ) for infinitly many m and a fixed j, then it is clear that U (x 0 ) ≤ g(x 0 ). Therefore, without loss of generality, we assume for all m and j that We have two cases. If there exists j ∈ {1, . . . , J} and a subsequence of c m such that U cm j (x j m )− ζ(x j m ) ≤ 0, then by passing to the limit m → ∞, we have U (x 0 ) − ζ(x 0 ) ≤ 0. Otherwise, by passing to subsequence, without loss of generality we can assume U cm j (x j m ) − ζ(x j m ) > 0 holds for all j and m. Then for each m ∈ N, we can choose j m ∈ {1, . . . , J} andx jm m such that and deduce from (4.4) that Our choice of j m implies (U cm jm −φ)(x jm m ) ≥ (U cm k −φ)(x jm m ) for all k = j m , and thus U cm jm (x jm m ) > M jm U cm (x jm m ). Consequently we obtain from (4.5) that Since we only have finite many choices of j m , by passing to a further subsequence if necessary, we can assume that j m → j 0 , then letting m → ∞ and using the continuity of the equation, we have Since α j 0 ∈ A δ is an admissible control, we obtain and conclude that U is a subsolution of (3.5).
We now proceed to show U is a supersolution. If φ ∈ C 1,2 and U − φ has a strict global mimimum atx 0 ∈Q T , then for any given j ∈ {1, . . . , J}, there will be sequences c m → 0, , and U cm j − φ attains a global mimimum atx m . Using the fact that U cm j is asupersolution to (3.5), we have (by ignoring the term U cm then passing m → ∞ enables us to conclude for any j ∈ {1, . . . , J}, which completes our proof.

General discrete approximation to the switching system
In this section, we establish the convergence of the piecewise constant policy approximation of (3.9) to the solution of the switching system (3.5). We will first summarize all the required conditions to guarantee the convergence, and perform the analysis under these assumptions. Then we will demonstrate in Section 4.4 that these conditions are in fact satisfied by the numerical scheme (3.17) proposed in Section 3 .
We assume the scheme (3.9) satisfies the following conditions introduced in [25]: (1) (Positive interpolation.) LetŨ n k,i(j) be the interpolant of the k-th grid onto the i-th point x n j,i of the j-th grid, and N k (j, i, n) be the neighbours 2 to the point x n j,i on the k-th grid Ω k,h . Then there exist weights {ω n k,i(j),a } a∈N k (j,i,n) satisfying ω n k,i(j),a ≥ 0 and a∈N k (j,i,n) ω n k,i(j),a = 1, such that we can writẽ U n k,i(j) = a∈N k (j,i,n) ω n k,i(j),a U n k,a . (4.6) (2) (Weak monotonicity.) The scheme (3.9) is monotone with respect to U n j,i andŨ n k,i(j) , i.e., if V n j,i ≥ U n j,i , ∀(i, j, n);Ṽ n k,i(j) ≥Ũ n k,i(j) , ∀(i, k, n), then we have (3) (ℓ ∞ stability.) The solution U n+1 j,i of the scheme (3.9) exists and is bounded uniformly in h and c.
Also note the contrast to the linear interpolant (3.10) used in (3.11) and (3.12) for the construction of a monotone approximation to the integral operators.
We now present the convergence of the discrete approximation to the switching system. The proof is essentially the same as that in [25] and is omitted. We remark that in the proof, we construct the solution of the switching system directly from the numerical solutions. Since the solution of the scheme (3.9) is uniformly bounded, Theorems 4.4 and 4.6 immediately give the existence and uniqueness of a bounded viscosity solution to the switching system (3.5).

A specific implicit scheme for the switching system
In this section, we analyze the implicit scheme (3.17) and demonstrate that it satisfies Condition 1, which subsequently implies its convergence to the switching system.
The following estimates are essential for our consistency and stability analysis.
Proof. We first derive the estimate for B α ε,h φ n+1 j,i . It follows from |η α | ≤ C and the definitions of B α ε,h φ and B α ε φ that We then infer from Taylor's theorem with an integral remainder that the truncation errors of the local terms can be bounded by for some function ω(x n+1 j,i , k) such that ω(·, k) → 0 as k → 0 uniformly on compact neighbourhoods of x n+1 j,i , which enables us to deduce that where κ α,n h,m,i , β α,n h,m,i , and d α,n h,k,m,i are defined in (3.13) and (3.15), respectively.
Proof. We shall only prove the estimate for κ α,n h,m,i , since the estimate for β α,n h,m,i follows from a similar argument, and the estimate for d α,n h,k,m,i follows directly from the fact that m ω m = 1. The definition of κ α,n h,m,i and the integrability property (2.1) of ν imply that Alternatively, it follows directly from the identity m∈Z d ω m (·; h) ≡ 1 that which leads us to the desired estimates.
Remark 4. Since we have not used any information on the exact behavior of the nonsingular measure ν around zero, the estimates for the nonlocal terms in Lemma 4.8 and 4.9 are not optimal for many specific cases. If one can estimate upper bounds of the density of the Lévy measure, or equivalently estimate the (pseudo-differential) orders of the nonlocal operators K α and B α , more precise results for the truncation error of the singular measure can be deduced ( [3]).
The next lemma presents some important properties of the Lax-Friedrichs numerical flux for Lipschitz continuous Hamiltonian, which are crucial for our subsequent analysis. We refer readers to [10] for a proof of these statements. Then the following hold:  .16) and (x n j,i , u, k) ∈ Ω j,h × R × R, and suppose Assumption 1 and the condition θ > Cλ hold, where C is the Lipschitz constant of the Hamiltonianf .

(1) (Consistency.) For any test functions
(2) (Monotonicity.) If V n j,i ≥ U n j,i , for all i, j, n, then we have (3) (Stability.) For any bounded functions U and V , we have Proposition 4.11. Suppose Assumption 1, the positive interpolation property in Condition 1 and the condition θ > Cλ hold. Then we have the following: (1) There exists a unique bounded solution U n of the scheme (3.17).
Proof. We start to establish the existence and uniqueness of a bounded solution of (3.17) in (1) (3.17), it suffices to establish that for small enough ρ, the operator P is a contraction on ℓ ∞ (Z d ), i.e., the Banach space of bounded functions on hZ d employed with the sup-norm, which along with the contraction mapping theorem leads to the desired results.
(Similar contraction operators have been introduced in [5,11] to demonstrate the well-posedness of their numerical schemes.) For any bounded functions U n j and V n j , the definitions of P, A α ε,h,k and K α,1 ε,h give that It remains to estimate (4.9), (4.10) and (4.11). Lemma 4.10 (3) enables us to bound (4.11) by −ρ2dθ(U n j,i − V n j,i ) + ρ2dθ|U n j − V n j | 0 . We then derive upper bounds for (4.9) and (4.10) depending the Lipschitz continuity of f in y enables us to bound (4.9) by ρ∆tC|U n j,i − V n j,i | = −ρ∆tC(U n j,i − V n j,i ). We then discuss the sign of B α ε,h U n j,i − B α ε,h V n j,i . Suppose B α ε,h U n j,i − B α ε,h V n j,i < 0, then we obtain from the monotonicity of f in k that (4.10) ≤ 0. Consequently we obtain that On the other hand, if B α ε,h U n j,i − B α ε,h V n j,i > 0, the Lipschitz continuity of f in k enables us to bound (4.10) by C(B α ε,h U n j,i − B α ε,h V n j,i ), which along with (3.12) implies again (4.12) provided that the the following condition is satisfied: (κ α,n h,m,i + β α,n h,m,i ) + C > 0, (4.13) which holds for small enough ρ. This completes the proof that P is a contraction operator.
We now proceed to establish the ℓ ∞ stability of the scheme. Let {U n−1 j } J j=1 be the solutions to (3.17). By expressing the discrete operators A α ε,h,k and K α,1 ε,h in the monotone form (3.14) and (3.11), and substituting them into (3.17), we obtain that Using similar arguments as those for the upper bound of (4.9), we deduce that (4.14) is bounded above by −∆tCU n j,i independent of the sign of U n j,i . Suppose now B α ε,h U n j,i < 0, then we obtain from the monotonicity of f in k that (4.15) is nonpositive. Then the ℓ ∞ stability of the numerical flux and the boundedness of f α (x, 0, 0, 0) yield that Here C is the constant from Assumption 1 and C 1 > 0 is a large enough constant that we will choose later. On the other hand, if B α ε,h U n j,i > 0, the Lipschitz continuity of f in k enables us to bound (4.15) by CB α ε,h U n j,i , which along with (3.11) implies again (4.16). With the estimate (4.16) in hand, we are ready to derive a uniform bound for the solutions {U n j }, which is independent of h and c. The proof follows from an inductive argument. Let us introduce the notation |U n | 0 = max 1≤j≤J |U n j | 0 for each n and define the term a 0 = max(|g| 0 , |ζ| 0 ), then it is clear that a 0 ≥ max(|U 0 | 0 , |ζ| 0 ). Suppose we have a n−1 such that a n−1 ≥ max(|U n−1 | 0 , |ζ| 0 ).
Then the definition of U n− 1 2 j,i implies that |U n− 1 2 j | 0 ≤ max(|ζ| 0 , |U n−1 | 0 ) ≤ a n−1 . Define the term a n := 1 1 + ∆tC a n−1 + ∆tC 1 , with the same constants as those in (4.16), then we have |U n | 0 ≤ a n . To proceed by induction, we further require a n ≥ |ζ| 0 . Since a n−1 ≥ |ζ| 0 and C is fixed, it suffices to require C 1 ≥ C|ζ| 0 . In this way, we can construct a sequence {a n }, such that |U n | 0 ≤ a n , but a n is uniformly bounded independent of c, h and ∆t, and hence this completes the proof of ℓ ∞ stability. We now study the weak monotonicity of the scheme. Let V n j,i ≥ U n j,i andṼ n k,i(j) ≥Ũ n k,i(j) for all i, j, k, n, then we have V Moreover the monotonicity of f in k and the weak monotonicity off imply that is nondecreasing with {U b+1 j,a } (a, b) = (i, n) , which gives the weak monotonicity of the scheme (3.17).
Finally we study the consistency of the scheme. By using the Lipschitz continuity of x → min(x, a), it is clear that it suffices to bound which can be estimated by using Lemma 4.8, Lemma 4.10, and the Lipschitz continuity of f .

Remark 5.
The contraction operator P is introduced to demonstrate our scheme admits a unique solution for any given discretization parameters ∆t, h, k and ε. However, due to its low convergence rate, it is not advisable to implement this contraction mapping directly to solve the nonlinear equation (3.17). In fact, Lemma 4.9 and the stability condition (4.12) restrict the contraction constant of P to admit a lower bound depending on the spatial discretization of the diffusion operator. This undesirable dependence of ∆t on k can be avoided by considering the mapping T defined by (3.18), which is implicit in the local terms. It has been shown that for small enough h, the contraction constant of T is proportional to θ, which can be chosen to achieve a rapid convergence.

Numerical experiments
In this section, we present several numerical experiments to analyse the effectiveness of the numerical scheme proposed in Section 3. We shall investigate the convergence of numerical solutions with respect to the switching cost, timestep, and mesh size, and show that a relatively coarse discretization of the admissible control set already leads to an accurate approximation.
We consider a portfolio optimization problem over a time interval [0, T ], in a framework of recursive utility. An investor can control his wealth process X t,x,α through a selection of the control process α ∈ A t t , say his or her portfolio strategy, and can also choose the duration of the investment via a stopping time τ . If the agent chooses a strategy pair (α, τ ), then the associated terminal reward is given by ξ t,x,α τ = ζ(τ, X t,x,α τ )1 t≤τ <T + g(X t,x,α τ )1 τ =T for some utilities ζ and g, and where τ ∈ T t t , the set of F t -stopping times valued in [t, T ]. The performance of this investment is evaluated under a particular nonlinear expectation, called the recursive utility process (see e.g. [7]), which is associated with a BSDE (with Lipschitz continuous drivers). It generalizes the standard additive utilities by including a dependence on the future utility (corresponding to the future wealth). Roughly speaking, the recursive utility depends on the future utility through the dependance of the driver f on y, and can also depend on the "variability" or "volatility" of future utility through the dependance of f on z and k.
Let x be the wealth at the initial time t, (α, τ ) be the chosen strategy, E α,t,x [·] be a recursive utility function associated with the BSDE with driver f α . The aim of the investor is to maximize the utility of the investment: over all admissible choices of (α, τ ). Under Assumption 1, it can be shown that the value function u of this mixed optimization problem coincides with the unique bounded viscosity solution of the (backward) HJBVI (2.5).
For the numerical tests, we consider a financial market with a risk-free asset with an interest rate r and a risky asset whose price follows where W is a Brownian motion andÑ (dt, de) = N (dt, de) − ν(de)dt is a compensated jump measure. If we denote by α t the percentage of the portfolio held in the risky asset at time t, then the dynamics of the portfolio is given by The performance will be evaluated by the recursive utility function induced by the BSDE with the following driver: for some instantaneous reward function ψ. Recall that any concave utility function admits a dual representation via a set of probability measures absolutely continuous with respect to the original probability measure P (see e.g. [21]). This result allows us to interpret κ ≥ 0 as an ambiguity-aversion coefficient relative to the Brownian motion as suggested in [7,Section 3.3].
The value function of this control problem satisfies the following HJBVI: We then specify the choice of data for our numerical experiments. We use the exponential utility function ζ(t, x) = g(x) = (1 − e −x ) + , which determines both the intermediate and terminal payoff, and acts as the initial condition and the obstacle to the HJBVI. Moreover, we consider the tempered stable Lévy measure ν(de) = e −µ|e| |e| de on R with intensity η(e) = 1 ∧ |e| for the jump component (which is a special case of the variance Gamma model in [8]). For simplicity, we choose a zero interest rate, i.e., r = 0.
We further choose the function ψ(t, x) = 0.8 exp(−(T − t)) exp(−x/2) as the instantaneous reward. As we will see later, this choice of ψ implies that the optimal control α varies in the state space and evolves in time, and there can be non-trivial stopping. The resulting HJBVI will be localized to the domain (0, 2) with u(t, x) = g(x) for (t, x) ∈ (0, T ) × R \ (0, 2). The numerical values for the parameters used in the experiments are given in Table 1.  Now we are ready to discuss the selection of the discretization parameters in detail. The density of the tempered stable measure ν enables us to improve the estimates in Lemma 4.9 to m =0 κ α,n h,m,i ≤ log(ε), and hence choosing ε = h and ∆t = O(h) leads us to a consistent approximation to the switching system (3.5). Moreover, choosing θ = 1 40 and ∆t = h 15 ensures the numerical flux is stable and the contraction constant of T in (3.18) is less than 1 10 . The coefficients of the nonlocal terms are evaluated by the midpoint quadrature formula, which is clearly monotone and consistent. We observe that for the control problem with the parameters as in Table 1, the optimal strategy α * will always be obtained at one of the endpoints of [0, 1]. In fact, using Taylor's theorem, we are able to approximate the nonlocal term K α u by at any given (t, x) for which the value function lies above the obstacle and is sufficiently smooth. Then we infer from the HJBVI (5.1) that the optimal control α * is the maximizer of a quadratic function on [0, 1], which is attained in the interior only if However, since we have b < σ, the above conditions can never hold for any x > 0. Consequently, we deduce that the admissible set is already finite, and replacing [0, 1] by A δ = {0, 1} in (5.1) will not introduce any discretiztion error. This has been confirmed with our numerical experiments. For the sake of simplicity, we discretise each component of the switching system on a single uniform mesh, thus Condition 1 (1) is trivially satisfied. Table 2 contains the numerical solutions to the last component of the switching system at the grid point (T, x 0 ) with different mesh size h and switching cost c. We examine the convergence of the numerical solutions, denoted as U h , in h for fixed c, as well as their convergence with respect to the cost c. For any fixed positive switching cost c, we infer from the lines (a) that the numerical solutions converge monotonically to the exact solution. Moreover, the lines (c) indicate the approximation error admits an asymptotic magnitude O(h) + O(∆t), which seems not to be affected by the size of the cost c. By considering the boldface values in Table 2 as an accurate approximation to the exact solution of the switching system with a given cost c, we can further conclude that the switching system is consistent to the HJBVI (5.1) with order 1. This follows from the approximate factor of four between the differences 0.00469, 0.00117, and 0.00029 between the last three pairs of values, proportional to the reduction in c. Therefore, by taking c = O(h) and ∆t = O(h), we can obtain a first-order scheme for the HJBVI.
We then proceed to analyze the effect of the control discretization. We pick the same parameters as those in Table 1, except that b = 0.25, which is chosen such that it is now possible that the optimal control is attained in the interior of (0, 1) (as seen from a similar argument as earlier). Computations are performed using Matlab R2016b on a 3.30GHz Intel Xeon ES-2667 16-Core processor with 256GB RAM to enable parallelization. Table 3 illustrates the numerical results for different control meshes (J = 1/δ + 1) with a fixed mesh size h = 0.005 and switching cost c = 1/2560, and also compares the runtime with or without parallelization.
We can clearly observe from line (a) second order convergence of the numerical solutions, and a relatively coarse control mesh has already yielded an accurate approximation with a negligible control discretization error.
Next, we discuss lines (b)-(f) which analyse the algorithm's parallel efficiency. Hereby, the implicit finite difference scheme for individual components of the switching system (i.e., (3.8), for different j) is solved independently on different processors, while the maximisation step (3.7) requires communication between processors.
The total execution time with and without parallelization are presented in line (b) and (c), respectively, which indicate a significant reduction of computational times. Moreover, by subtracting the communication time among clusters, as shown in line (d), from the total runtime, we can obtain the actual time spent on executing the numerical scheme (line (e)). The speed-up rate of the parallelization is shown in line (f), which grows with the number of controls, and remains stable at the number of cores. Therefore, together with parallelization, piecewise constant timestepping enables us to achieve a high accuracy in the control discretization without significantly increasing the computational time, which is an advantage over policy iterations, which do not parallelise naturally.
We finally examine the impact of the computational domain by performing computations on (0, 3) with h = 1/400, ∆t = h/20, c = 1/640 and the parameters as in Table 1. Compared to the results in Table 2, this larger domain leads to a relative difference of 7.53 · 10 −7 , which is negligible compared to the time and spatial discretization errors.  The numerical value function and the corresponded feedback control strategy with J = 21 are presented in Figure 1, in which the white area represents the region where the obstacle is active, and otherwise the colour indicates the value of the optimal control, as shown in the panel on the right. The approximation to the optimal control pair (τ, α) was found from the numerical solution as follows (see also (3.7) and Remark 2), noting that in our tests x j,i = x k,i for all j, k, and therefore no interpolation is needed: where α n i = α i * n is an approximation to the optimal policy and {(t n , x 1,i ) : θ n i = 1} is an approximation to the stopping region.

Conclusions
This paper provides a PDE approximation scheme for the value function of a mixed stochastic control/optimal stopping problem with nonlinear expectations and infinite activity jumps, which is the unique viscosity solution of a nonlocal HJB variational inequality. The approach that we have adopted is based on piecewise constant policy time stepping (PCPT), which reduces the problem to a system of semi-linear PDEs, and a monotone approximation scheme. We prove the convergence of the numerical scheme and illustrate the theoretical results with some numerical examples in the case of a recursive utility maximisation problem.
To the best of our knowledge, this is the first paper which proposes a numerical approximation for a control problem in such a generality. Natural next steps would be to derive theoretical results on the convergence rate and to extend this approach to the case of Hamilton-Jacobi-Bellman-Isaac equations obtained in [6]. and consequently we obtain from the Lipschitz continuity of f that for a suitable defined ω 1 (x, δ) with the properties of ω 0 (x, δ). Therefore, using (A.1), (A.2) and the fact that A δ ⊂ A, we have which completes the proof of our desired result.

B Comparison principle for switching systems
In this section, we establish the comparison principle for switching system (3.5), cf. Theorem 4.4. We consider a slightly more general switching system with no truncation of the singular measure in K ε and B ε , which includes as a special case the switching system (3.5). We first use a classical no-loop argument to reduce the problem into scalar cases, and then analyze the scalar HJBVI by extending the results for continuous solutions in [14] to semicontinuous viscosity solutions. For simplicity, we denote by σ the modified diffusion coefficientσ α defined as (3.1).

Proof of Theorem 4.4. Set
It suffices to show that M ≤ 0. For any given ε, ρ > 0, we introduce the functions for each t, s ∈ [0, T ] and x, y ∈ R d , and define the quantity M ε,ρ := sup j,t,s,x,y ψ ε,ρ j (t, s, x, y).
We now divide our analysis into three cases to establish M ≤ 0. If there exists a subsequence of {t ρ } such that t ρ = 0 for all ρ, we then deduce M ≤ 0 along this subsequence by adapting the arguments in [14] to semicontinuous solutions.
On the other hand, if t ρ is different from 0 for all ρ, then for any fixed ρ and small enough ε, using Lemma B.2, which can be proved similarly as Lemma 4.1 in [4], we know there exists j ε,ρ 0 ∈ {1, . . . , J}, which for simplicity is still denoted as j ε,ρ , such that U j ε,ρ (t ε,ρ , x ε,ρ ) > M j ε,ρ U (t ε,ρ , x ε,ρ ). In other words, at the point (t ε,ρ , s ε,ρ , x ε,ρ , y ε,ρ ), by considering the j ε,ρ component of the switching system, we can without loss of generality ignore the term U j ε,ρ − M j ε,ρ U in the definition of subsolutions and get back to the scalar HJBVI.
Then noticing the estimates derived in [14] for each term on the right-hand side of the above expression are uniform in the control α j ρ , and successively passing δ, ε and ρ to 0, we deduce that 0 < M ≤ 0, which leads to a contradiction. Thus we conclude M ≤ 0 and complete the proof.

C Truncation of singular measures
A possible way to work with a nonsingular jump measure is to introduce a Backward SDE with a modified driver and an approximative jump-diffusion dynamics where the small jumps part has been substituted by a rescaled diffusion coefficient of the Brownian motion W .
More precisely, we adopt the same probability space as introduced in Section 2, which supports the Brownian motion process W and the independent Poisson measure N (dt, de). For a given jump truncation size ε > 0, we define a modified diffusion coefficientσ α as in (3.1), and also introduce a modified driver f ε (α, t, x, y, z, k) :=f (α, t, x, y, z, |e|≥ε k(e)γ(α, x, e)ν(de)), where the function f is given in Assumption 1.
For any given initial state x ∈ R d , control α ∈ A t t and τ ∈ T t t , we consider the modified controlled jump-diffusion process (X ε,α,t,x s ) t≤s≤T satisfying the following SDE: for each s ∈ [t, T ], (e)Ñ t (ds, de), Y ε,α,t,x τ,τ = ξ(τ, X ε,α,t,x τ ). (C. 2) The coefficients of the above SDE and BSDE satisfy Assumption 1, and therefore the equations are well-posed. Now we are ready to state the modified mixed optimal stopping and control problem. For each initial time t ∈ [0, T ] and initial state x ∈ R d , we consider the following value function: with κ(ε) := |e|≤ε (1 ∧ |e| 2 )ν(de) and C a constant independent of α and ε.