A deep neural network algorithm for semilinear elliptic PDEs with applications in insurance mathematics

In insurance mathematics optimal control problems over an infinite time horizon arise when computing risk measures. Their solutions correspond to solutions of deterministic semilinear (degenerate) elliptic partial differential equations. In this paper we propose a deep neural network algorithm for solving such partial differential equations in high dimensions. The algorithm is based on the correspondence of elliptic partial differential equations to backward stochastic differential equations with random terminal time.


Introduction
Classical optimal control problems in insurance mathematics include finding risk measures like the probability of ruin or the expected discounted future dividend payments. While in mathematical finance control problems (e.g., investment problems) are studied over relatively short time horizons, in insurance mathematics they are often considered over the whole lifetime of the insurance company. A standard method for solving control problems is to derive the associated Hamilton-Jacobi-Bellman (HJB) equation -a semilinear (sometimes integro) partial differential equation (PDE) the solution of which solves the control problem. In the case of infinite time horizon problems, these HJB equations are (degenerate) elliptic. In this paper we propose a deep neural network algorithm for semilinear (degenerate) elliptic PDEs in high dimensions.
We need this method for solving the dividend maximization problem in insurance mathematics. This problem originates in the seminal work by De Finetti [20], who introduced expected discounted future dividend payments as a valuation principle for a homogeneous insurance portfolio. This constitutes an alternative risk measure to the (even more) classical probability of ruin. Classical results on the dividend maximization problem are [61,42,56,2]. Overviews can be found in [1,3], for an introduction to optimization problems in insurance we refer to [60,4]. The problem becomes high-dimensional, if, for example, we allow for changes in the underlying economy. Such models have been studied, e.g., in [44,63,66,47,64,65]. In this paper we solve the problem studied in [65] in high-dimensions, where classical numerical methods fail.
Classical algorithms for solving semilinear (degenerate) elliptic PDEs like finite difference or finite element methods suffer from the so-called curse of dimensionality -the computational complexity for solving the discretized equation grows exponentially in the dimension. In highdimensions (say > 10) one has to resort to costly quadrature methods such as multilevel-Monte Carlo or the quasi-Monte Carlo-based method presented in [45]. In recent years, deep neural network (DNN) algorithms for high-dimensional PDEs have been studied extensively. Prominent examples are [33,22], where semilinear parabolic PDEs are associated with backward stochastic differential equations (BSDEs) through the (non-linear) Feynman-Kac formula and a DNN algorithm is proposed that solves these PDEs by solving the associated BSDEs. In the literature there exists a variety of DNN approaches for solving PDEs, in particular (degenerate) parabolic ones. Great literature overviews are given, e.g., in [32,8], out of which we list some contributions here: [5,6,7,10,11,12,16,17,21,23,25,26,28,34,35,36,37,41,48,49,50,51,55,57,62]. For elliptic PDEs, a multi-level Picard iteration algorithm is studied in [8], a derivative-free method using Brownian walkers without explicit calculation of the derivatives of the neural network is studied in [35], and a walk-on-the-sphere algorithm is introduced in [29] for the Poisson equation, where the existence of DNNs that are able to approximate the solution to certain elliptic PDEs is shown.
In the present article, we extend the DNN algorithm from [33] to a large class of semilinear (degenerate) elliptic PDEs by using a correspondence between these PDEs and BSDEs with random terminal time. This correspondence was first presented in [53], and elaborated afterwards, e.g., in [19,15,54,59,18]. As these results are not as standard as the BSDE correspondence to parabolic PDEs, we summarize the theory in Section 2 for the convenience of the reader.

BSDEs associated with elliptic PDEs
This section contains a short survey on scalar backward stochastic differential equations with random terminal times and on how they are related to a certain type of semilinear elliptic partial differential equations.

BSDEs with random terminal times
Let (Ω, F, P, (F t ) t∈[0,∞) ) be a filtered probability space satisfying the usual conditions and let W = (W t ) t∈[0,∞) be a d-dimensional standard Brownian motion on it. We assume that (F t ) t∈[0,∞) is equal to the augmented natural filtration generated by W . For all real valued row or column vectors x, let |x| denote their Euclidean norm. We need the following notations and definitions for BSDEs.
Optimal control problems which can be treated using a BSDE setting have for example been studied in [18,Section 6]. For this they consider generators of the forward-backward form where X is a forward diffusion (see also the notation in the following subsection), U is a Banach space, r is a Hilbert space-valued function (in their setting z takes values in the according dual space) with linear growth, g a real valued function with quadratic growth in u, and λ ∈ (0, ∞).
In the sequel we focus on generators of forward-backward form.

Semilinear elliptic PDEs and BSDEs with random terminal time
In this subsection we recall the connection between semilinear elliptic PDEs and BSDEs with random (and possibly infinite) terminal time. The relationship between the theories is based on a nonlinear extension of the Feynman-Kac formula, see [53,Section 4].
We define the forward process X as the stochastic process satisfying a.s., where x ∈ R d and µ : R d → R d and σ : R d → R d×d are globally Lipschitz functions.
In this paper we consider the following class of PDEs.
where the differential operator L acting on C 2 (R d ) is given by and F : R d × R × R 1×d → R is such that the process (t, y, z) → F (X t , y, z) is a generator of a BSDE in the sense of Definition 2.1.
• We say that a function u satisfies equation (4) with Dirichlet boundary conditions on the open, bounded domain G ⊆ R d , if where g : R d → R is a bounded, continuous function. (3), and the solution satisfies a.s. for all T ∈ (0, ∞) that 2. A BSDE associated to the PDE (6) with Dirichlet boundary conditions is given by the (3), and the solution satisfies a.s. for all T ∈ (0, ∞) that In order to keep the notation simple, we do not highlight the dependence of X, Y, Z on x.
For later use we also introduce the following notion of solutions of PDEs, which we will use later.
• A function u ∈ C(R d ) is called viscosity subsolution of (4), if for all ϕ ∈ C 2 (R d ) and all points x ∈ R d where u − ϕ has a local maximum, (4), if for all ϕ ∈ C 2 (R d ) and all points x ∈ R d where u − ϕ has a local minimum, (4), if it is a viscosity sub-and supersolution.
A similar definition of viscosity solutions can be given for the case of Dirichlet boundary conditions (6), see [53].
For later use, note that (8) can be rewritten in forward form as The following theorems link the semilinear elliptic PDEs (4) and (6) to the associated BSDEs.
solve the BSDE (7). An equivalent statement can be established for the system with boundary conditions (6) and equation (8), see [53].
Note that for all x ∈ R d , Y and Z are stochastic processes adapted to . Assume that for some K, K , p ∈ (0, ∞), γ ∈ (−∞, 0) the function F satisfies for all x, y, y , z, z , To conclude, the correspondence between PDE (6) and BSDE (8) is given by . For tackling elliptic PDEs which are degenerate (as it is the case for our insurance mathematics example) we need to take the relationship a little further in order to escape the not so convenient structure of the Z-process. We factor Zσ(X) = Z for cases where this equation is solvable for Z (σ needs not necessarily be invertible) and define f (x, y, ζ) := F (x, y, ζσ(x)) 1 , giving the correspondence Y t = u(X t ), Z t = ∇u(X t ), ξ = g(X τ ). This relationship motivates us to solve semilinear degenerate elliptic PDEs by solving the corresponding BSDEs forward in time (cf. (9)) for Y 0 by approximating the paths of Z = ∇u(X) by a DNN, see Section 3. Doing so, we obtain an estimate of a solution value u(x) for a given x ∈ R d .

Algorithm
The idea of the proposed algorithm is inspired by [33], where the authors use BSDEs in order to solve semilinear parabolic PDEs. In the same spirit, we construct a DNN algorithm to solve BSDEs related to semilinear elliptic PDEs. The details of the algorithm are described in three steps of increasing specificity. Here we give a rough idea. Algorithm 1 below provides a pseudocode. Our program code builds upon the code from [33] and is published on Github 2 . The goal is to calculate u(x) for some given x ∈ G. Let 0 = t 0 < t 1 < · · · < t N = T , ∆t n = t n+1 − t n , and ∆W n = W t n+1 − W tn . First we simulate M paths ω 1 , . . . , ω M of the Brownian motion W to approximate the forward process using the Euler-Maruyama scheme, that is X 0 = x and Using these paths we compute the solution to the PDE via the BSDE forward in time until the path of X hits the stopping time τ by For all t n , Z tn = ∇u(X tn ) are approximated by DNNs, each mapping G to R d . Therefore, also u(X t n+1 ) is approximated by a DNN as a combination of DNNs. At the cutoff time T (which we choose sufficiently large) we compare u(X τ ∧T ) with the terminal value ξ. This defines the loss function, which is used for the training of our DNNs: Though the last point is already studied in the literature (see, e.g., [8,9,13 We close this section with some comments on the implementation.

Remark 3.2.
• All DNNs are initialized with random numbers.
• For each value of x we average u(x) over 5 independent runs. The estimator for u(x) is calculated as the mean value of u(x) in the last 3 network training epochs of each run, sampled according to the validation size.  • For efficient implementation in a multi-core environment we create a stopping matrix that contains the values 0 and 1 for each sample ω and each time step t n , in dependence of whether the path X(ω) was stopped or not.
• We choose a non-equidistant time grid in order to get a higher resolution for earlier (and hence probably closer to the stopping time) time points.
• We use tanh as activation function.
• We compute u(x) simultaneously for 8 values of x by using parallel computing.

Examples
In this section we apply the proposed algorithm to three examples. The first one serves as a validity check, the second one as an academic example with a non-linearity. Finally, we apply the algorithm to solve the dividend maximization problem under incomplete information.

The Poisson equation
The first example we study is the Poisson equation -a linear PDE.
To obtain a reference solution for this linear BSDE we use an analytic formula for the expectation of τ , see [52,Example 7.4.2,p. 121]. This yields

Numerical results
We compute u(x) on the R 2 and the R 100 for 15 different values of x. Figure 1 shows Note that as the expected value of τ decreases linearly in d, we adapt the cut off time T for d = 100 accordingly.

Quadratic gradient
The second example is a semilinear PDE with a quadratic gradient term. Let r ∈ (0, ∞), G = x ∈ R d : |x| < r , and ∂G = x ∈ R d : |x| = r . We consider the PDE corresponding to the BSDE f (x, y, ζ) = |ζ| 2 − 2e −y , ξ = log Also for this example we have an analytic reference solution given by u(x) = log |x| 2 + 1 d .

Numerical results
As in the previous example we compute u(x) for 15 different values of x on the R 2 and the R 100 . Figure 2 shows

Dividend maximization
The goal of this paper is to show how to use the proposed DNN algorithm to solve highdimensional control problems that arise in insurance mathematics. We finally arrived at the point where we are ready to do so.
Our example comes from [65], where the author studies De Finetti's dividend maximization problem in a setup with incomplete information about the current state of the market. The hidden market-state process determines the trend of the surplus process of the insurance company and is modelled as a d-state Markov chain. Using stochastic filtering, in [65] they achieve to transform the one-dimensional problem under incomplete information to a d-dimensional problem under complete information. The cost is (d − 1) additional dimensions in the state space. We state the problem under complete information using different notation than in [65] in order to avoid ambiguities. The probability that the Markov chain modelling the market-state is in state i ∈ {1, . . . , d − 1} is given by where x i ∈ (0, 1), B is a one-dimensional Brownian motion, a 1 , . . . , a d ∈ R are the values of the surplus trend in the respective market-states of the hidden Markov chain, and (q i,j ) i,j∈{1,...,d} ∈ R d×d denotes the intensity matrix of the chain. Let ( t ) t∈[0,∞) be the dividend rate process. The surplus of the insurance company is given by where x d , ρ ∈ (0, ∞). For later use we define also the dividend free surplus The processes (16) and (18) form the d-dimensional state process underlying the optimal control problem we aim to solve.
The goal of the insurance company is to determine its value by maximizing the discounted dividends payments until the time of ruin η = inf{t ∈ (0, ∞] : X d t < 0}, that is it seeks to find where δ ∈ (0, ∞) is a discount rate, A is the set of admissible controls, and E x 1 ,...,x d [·] denotes the expectation under the initial conditions π i (0) = x i for i ∈ {1, . . . , d − 1} and X d 0 = x d . Admissible controls are (F X d t ) t≥0 -progressively measurable, [0, K]-valued for K ∈ (0, ∞), and fulfill t ≡ 0 for t > η, cf. [65].
In order to tackle the problem, we solve the corresponding Hamilton-Jacobi-Bellmann (HJB) equation 4 from [65], where L is the second order degenerate elliptic operator The supremum in (21) is attained at Plugging this into (21) we end up with a d-dimensional semilinear degenerate elliptic PDE: The boundary conditions in x d direction are given by No boundary conditions are required for the other variables, cf. [65].
In [65,Corollary 3.6] it is proven that the unique viscosity solution to (22) solves the optimal control problem (20). Hence, we can indeed solve the control problem by solving the HJB equation.
For the numerical approximation we cut off x d at r ∈ (0, ∞). Hence, For the convenience of the reader we derive the BSDE corresponding to (22). The forward equation is given by where W = (B, W 2 , . . . , W d ) , x = (x 1 , . . . , x d ) , and We claim that the BSDE associated to (22) is given in forward form by Applying Itô's formula to u(X) yields Combining (23) and (24) gives Cancelling terms verifies (in a heuristic manner) (22). Hence, the corresponding BSDE has the parameters and ξ = K/δ, X d τ = r, 0, X d τ = 0, if τ < ∞.

Numerical results
As for this example we have no analytic reference solution at hand, we use the solution from [65] for the case d = 2, which was obtained by a finite difference method and policy iteration. Then we show that the DNN algorithm also provides an approximation in high dimensions in reasonable computation time. As in the previous examples we compute u(x) on the R 2 and on the R 100 for 15 different values of x. Figure 3 shows the approximate solution of the HJB equation and hence the value of the insurance company obtained by the DNN algorithm (in blue) and the reference solution from [65] (in green) for the case d = 2. For d = 100 we have no reference solution at hand. Figure 4 shows the loss for a fixed value of x in the case d = 100. The parameters we use are Figure 3: Approximate solution (blue) and reference solution (green) for the dividend problem on the R 2 for fixed π 1 = π 2 = 0.5 (left) and on the R 100 for fixed π 1 = · · · = π 100 = 0.01 (right, without reference solution).

Conclusion
The goal of this paper was to extend the DNN algorithm proposed by Han, Jetzen, and E [33] to the case of (degenerate) elliptic semilinear PDEs as they appear frequently in insurance mathematics, where often problems over an infinite time horizon are considered. We attacked the problem inspired by a series of results by Pardoux [53]. Of course, in low dimensions one would not use the proposed DNN algorithm -classical methods are much more efficient. However, if a solution in high dimensions is needed, a DNN algorithm like the one presented here can still be applied in order to get an idea of how the solution looks like, while classical methods suffer from the curse of dimensionality.