A deep neural network algorithm for semilinear elliptic PDEs with applications in insurance mathematics
AA deep neural network algorithm for semilinear elliptic PDEswith applications in insurance mathematics
Stefan Kremsner Alexander Steinicke Michaela Sz¨olgyenyi Preprint, December 2020
Abstract
In insurance mathematics, optimal control problems over an infinite time horizon arisewhen computing risk measures. An example of such a risk measure is the expected dis-counted future dividend payments. In models which take multiple economic factors into ac-count, this problem is high-dimensional. The solutions to such control problems correspondto solutions of deterministic semilinear (degenerate) elliptic partial differential equations. Inthe present paper we propose a novel deep neural network algorithm for solving such partialdifferential equations in high dimensions in order to be able to compute the proposed riskmeasure in a complex high-dimensional economic environment. The method is based on thecorrespondence of elliptic partial differential equations to backward stochastic differentialequations with unbounded random terminal time. In particular, backward stochastic differ-ential equations which can be identified with solutions of elliptic partial differential equationsare approximated by means of deep neural networks.
Keywords:
Backward stochastic differential equations, semilinear elliptic partial differentialequations, stochastic optimal control, unbounded random terminal time, machine learning, deepneural networks.
Mathematics Subject Classification (MSC 2020):
Classical optimal control problems in insurance mathematics include finding risk measures likethe probability of ruin or the expected discounted future dividend payments. Mathematically,these are problems over a potentially infinite time horizon, ending at an unbounded randomterminal time – the time of ruin of the insurance company. In recent models which take mul-tiple economic factors into account, the problems are high dimensional. For computing theserisk measures, optimal control problems need to be solved numerically. A standard method forsolving control problems is to derive the associated Hamilton-Jacobi-Bellman (HJB) equation– a semilinear (sometimes integro) partial differential equation (PDE) and show that its (nu-merical) solution also solves the original control problem. In the case of infinite time horizonproblems, these HJB equations are (degenerate) elliptic.
In this paper we propose a novel deepneural network algorithm for semilinear (degenerate) elliptic PDEs associated to infinite timehorizon control problems in high dimensions. Department of Mathematics, University of Graz, Heinrichstraße 36, 8010 Graz, Austria. stefan.kremsner @ uni-graz.at (cid:12) Department of Mathematics and Information Technology, Montanuniversitaet Leoben, Peter Tunner-Straße25/I, 8700 Leoben, Austria. alexander.steinicke @ unileoben.ac.at Department of Statistics, University of Klagenfurt, Universit¨atsstraße 65-67, 9020 Klagenfurt, Austria. michaela.szoelgyenyi @ aau.at a r X i v : . [ q -f i n . M F ] D ec e apply this method to solve the dividend maximization problem in insurance mathemat-ics. This problem originates in the seminal work by De Finetti [20], who introduced expecteddiscounted future dividend payments as a valuation principle for a homogeneous insurance port-folio. This constitutes an alternative risk measure to the (even more) classical probability ofruin. Classical results on the dividend maximization problem are [64, 44, 58, 2]. Overviews canbe found in [1, 3], for an introduction to optimization problems in insurance we refer to [63, 4].Recent models for the surplus of an insurance company allow for changes in the underlyingeconomy. Such models have been studied, e.g., in [46, 66, 69, 49, 67, 68, 61]. In [49, 67, 68]hidden Markov models for the underlying economic environment were proposed that allow fortaking (multiple) exogenous, even not directly observable, economic factors into account. Whilethese authors study the dividend maximization problem from a theoretical perspective, we areinterested in computing the risk measure. However, classical numerical methods fail when theproblem becomes high-dimensional, that is for example when exogenous economic factors aretaken into account. In this paper we propose a novel deep neural network algorithm to solvehigh-dimensional problems. As an application we use it to solve the dividend maximizationproblem in the model from [68] in high dimensions numerically.Classical algorithms for solving semilinear (degenerate) elliptic PDEs like finite difference orfinite element methods suffer from the so-called curse of dimensionality – the computationalcomplexity for solving the discretized equation grows exponentially in the dimension. In high-dimensions (say >
10) one has to resort to costly quadrature methods such as multilevel-MonteCarlo or the quasi-Monte Carlo-based method presented in [47]. In recent years, deep neural net-work (DNN) algorithms for high-dimensional PDEs have been studied extensively. Prominentexamples are [33, 22], where semilinear parabolic PDEs are associated with backward stochasticdifferential equations (BSDEs) through the (non-linear) Feynman-Kac formula and a DNN al-gorithm is proposed that solves these PDEs by solving the associated BSDEs. In the literaturethere exists a variety of DNN approaches for solving PDEs, in particular (degenerate) parabolicones. Great literature overviews are given, e.g., in [32, 8], out of which we list some contribu-tions here: [5, 6, 7, 10, 11, 12, 16, 17, 21, 23, 25, 26, 28, 34, 35, 36, 39, 43, 50, 51, 52, 53, 57, 59, 65].While in mathematical finance control problems (e.g., investment problems) are studied over rel-atively short time horizons, leading to (degenerate) parabolic PDEs, in insurance mathematicsthey are often considered over the whole lifetime of the insurance company, leading to (degener-ate) elliptic PDEs. 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 derivativesof 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 solutionto certain elliptic PDEs is shown.In the present article we propose a novel DNN algorithm for a large class of semilinear (de-generate) elliptic PDEs. For this, we adopt the approach from [33] for (degenerate) parabolicPDEs. The difference here is that we use the correspondence between the PDEs we seek to solveand BSDEs with random terminal time. This correspondence was first presented in [55], andelaborated afterwards, e.g., in [19, 15, 56, 62, 18].As these results are not as standard as the BSDE correspondence to parabolic PDEs, we sum-marize the theory in Section 2 for the convenience of the reader. In Section 3 we present theDNN algorithm, and test it in Sections 4.1 and 4.2. In Section 4.3 we present the model from[68] in which we seek to solve the dividend maximization problem and hence to compute therisk measure. That this method works also in high dimensions is demonstrated at the end ofSection 4.3, where numerical results are presented.2he method presented here can be applied to many other high-dimensional semilinear (degener-ate) elliptic PDE problems in insurance mathematics, such as the calculation of ruin probabili-ties, but we emphasize that its application possibilities are not limited to insurance problems.
This section contains a short survey on scalar backward stochastic differential equations withrandom terminal times and on how they are related to a certain type of semilinear elliptic partialdifferential equations.
Let (Ω , F , P , ( F t ) t ∈ [0 , ∞ ) ) be a filtered probability space satisfying the usual conditions andlet 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 rowor column vectors x , let | x | denote their Euclidean norm. We need the following notations anddefinitions for BSDEs. Definition 2.1. A BSDE with random terminal time is a triple ( τ, ξ, f), where • the terminal time τ : Ω → [0 , ∞ ] is an ( F t ) t ∈ [0 , ∞ ) -stopping time, • the generator f : Ω × [0 , ∞ ) × R × R × d → R is a process which satisfies that for all( y, z ) ∈ R × R × d , the process t (cid:55)→ f( t, y, z ) is progressively measurable, • the terminal condition ξ : Ω → R is an F τ -measurable random variable with ξ = 0 on { τ = ∞} . Definition 2.2. A solution to the BSDE ( τ, ξ, f) is a pair of progressively measurable processes( Y, Z ) = (cid:0) ( Y t ) t ∈ [0 , ∞ ) , ( Z t ) t ∈ [0 , ∞ ) (cid:1) with values in R × R × d , where • Y is continuous P -a.s. and for all T ∈ (0 , ∞ ), the trajectories t (cid:55)→ Z t belong to L ([0 , T ] , R × d ),and t (cid:55)→ f( t, Y t , Z t ) is in L ([0 , T ]), • for all T ∈ (0 , ∞ ) and all t ∈ [0 , T ] it holds a.s. that Y t = Y T + (cid:90) T ∧ τt ∧ τ f( s, Y s , Z s ) ds − (cid:90) T ∧ τt ∧ τ Z s dW s , (1) • Y t = ξ and Z t = 0 on { t ≥ τ } .Results on existence of solutions of BSDEs with random terminal time can be found in Pardoux’seminal article [55] (see [55, Theorem 3.2]), in [18] for generators with quadratic growth (see[18, Theorem 3.3]), and, e.g., in [19, 15, 56, 14, 62]; many of them cover multidimensional statespaces for the Y -process.Optimal control problems which can be treated using a BSDE setting have for example beenstudied in [18, Section 6]. For this they consider generators of the forward-backward formf( t, y, z ) = F ( X t , y, z ) = inf { g ( X t , u ) + zr ( X t , u ) : u ∈ U } − λy, (2)where X is a forward diffusion (see also the notation in the following subsection), U is a Banachspace, r is a Hilbert space-valued function (in their setting z takes values in the according dualspace) 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.3 .2 Semilinear elliptic PDEs and BSDEs with random terminal time In this subsection we recall the connection between semilinear elliptic PDEs and BSDEs withrandom (and possibly infinite) terminal time. The relationship between the theories is based ona nonlinear extension of the Feynman-Kac formula, see [55, Section 4].We define the forward process X as the stochastic process satisfying a.s., X t = x + (cid:90) t µ ( X s ) ds + (cid:90) t σ ( X s ) dW s , t ∈ [0 , ∞ ) , (3)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. Definition 2.3. • A semilinear (degenerate) elliptic PDE on the whole R d is of the form L u + F ( · , u, ( ∇ u ) σ ) = 0 , (4)where the differential operator L acting on C ( R d ) is given by L := 12 d (cid:88) i,j =1 ( σσ (cid:62) ) i,j ( x ) ∂ ∂x i ∂x j + d (cid:88) i =1 µ i ( x ) ∂∂x i , (5)and F : R d × R × R × d → R is such that the process ( t, y, z ) (cid:55)→ F ( X t , y, z ) is a generatorof a BSDE in the sense of Definition 2.1. • We say that a function u satisfies equation (4) with Dirichlet boundary conditions on theopen, bounded domain G ⊆ R d , if L u + F ( · , u, ( ∇ u ) σ ) = 0 , x ∈ G,u ( x ) = g ( x ) , x ∈ ∂G, (6)where g : R d → R is a bounded, continuous function. Definition 2.4.
1. A
BSDE associated to the PDE (4) on the whole R d is given by thetriplet ( τ, ξ, f), where τ ≡ ∞ , ξ = 0, f( t, y, z ) = F ( X t , y, z ), X is as in (3), and the solutionsatisfies a.s. for all T ∈ (0 , ∞ ) that Y t = Y T + (cid:90) Tt F ( X s , Y s , Z s ) ds − (cid:90) Tt Z s dW s , t ∈ [0 , T ] . (7)2. A BSDE associated to the PDE (6) with Dirichlet boundary conditions is given by thetriplet ( τ, g ( X τ ) , f), where τ = inf { t ∈ [0 , ∞ ) : X t / ∈ G } , f( t, y, z ) = F ( X t , y, z ), X is as in(3), and the solution satisfies a.s. for all T ∈ (0 , ∞ ) that Y t = Y T + (cid:90) T ∧ τt ∧ τ F ( X s , Y s , Z s ) ds − (cid:90) T ∧ τt ∧ τ Z s dW s , t ∈ [0 , T ] ,Y t = g ( X τ ) , Z t = 0 , t ≥ τ. (8)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.4 efinition 2.5. • A function u ∈ C ( R d ) is called viscosity subsolution of (4), if for all ϕ ∈ C ( R d ) and all points x ∈ R d where u − ϕ has a local maximum, L ϕ ( x ) + F ( x, u ( x ) , ( ∇ ϕ ( x )) σ ( x )) ≥ . • A function u ∈ C ( R d ) is called viscosity supersolution of (4), if for all ϕ ∈ C ( R d ) and allpoints x ∈ R d where u − ϕ has a local minimum, L ϕ ( x ) + F ( x, u ( x ) , ( ∇ ϕ ( x )) σ ( x )) ≤ . • A function u ∈ C ( R d ) is called viscosity solution of (4), if it is a viscosity sub- andsupersolution.A similar definition of viscosity solutions can be given for the case of Dirichlet boundary condi-tions (6), see [55].For later use, note that (8) can be rewritten in forward form as Y t = Y − (cid:90) t F ( X s , Y s , Z s ) ds + (cid:90) t Z s dW s , t ∈ [0 , τ ) ,Y t = g ( X τ ) , Z t = 0 , t ≥ τ. (9)The following theorems link the semilinear elliptic PDEs (4) and (6) to the associated BSDEs. Theorem 2.6 ([55, Theorem 4.1]) . Let ( t, y, z ) (cid:55)→ F ( X t , y, z ) meet the assumptions of [55,Theorem 3.2] and let u ∈ C ( R d ) satisfy E (cid:20)(cid:90) ∞ e λt | (( ∇ u ) σ )( X t ) | dt (cid:21) < ∞ with λ as in [55, Theorem 3.2]. If u is a classical solution of the PDE (4) , then Y t = u ( X t ) , Z t = (( ∇ u ) σ )( X t ) solve the BSDE (7) . An equivalent statement can be established for the system with boundaryconditions (6) and equation (8) , see [55]. Note that for all x ∈ R d , Y and Z are stochastic processes adapted to ( F t ) t ∈ [0 , ∞ ) . Therefore Y , Z are F -measurable and hence a.s. deterministic. For us, the connection between PDEs andBSDEs is of relevance because of the converse result, where x (cid:55)→ Y delivers a solution to therespective PDE. Theorem 2.7 ([55, Theorem 4.3]) . Assume that for some
K, K (cid:48) , p ∈ (0 , ∞ ) , γ ∈ ( −∞ , thefunction F satisfies for all x, y, y (cid:48) , z, z (cid:48) , (i) | F ( x, y, z ) | ≤ K (cid:48) (1 + | x | p + | y | + | z | ) ,(ii) (cid:104) y − y (cid:48) , F ( x, y, z ) − F ( x, y (cid:48) , z ) (cid:105) ≤ γ | y − y (cid:48) | ,(iii) | F ( x, y, z ) − F ( x, y, z (cid:48) ) | ≤ K | z − z (cid:48) | .Then [55, Theorem 3.2] can be applied to the generator ( t, y, z ) (cid:55)→ F ( X t , y, z ) , showing that thefunction u given by u ( x ) = Y is a viscosity solution to (4) , where Y is the first component ofthe unique solution to (7) in the class of solutions from [55, Theorem 3.2]. G and theexit time τ from (8). We refer to [55, Theorem 4.3]. A corresponding result for BSDEs withquadratic generator is [18, Theorem 5.2].To conclude, the correspondence between PDE (6) and BSDE (8) is given by Y t = u ( X t ), Z t = (( ∇ u ) σ )( X t ), ξ = g ( X τ ). For tackling elliptic PDEs which are degenerate (as it is thecase for our insurance mathematics example) we need to take the relationship a little furtherin 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) anddefine f ( x, y, ζ ) := F ( x, y, ζσ ( x )) , 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 solvingthe corresponding BSDEs forward in time (cf. (9)) Y t = Y − (cid:90) t f ( X s , Y s , Z s ) ds + (cid:90) t Z s σ ( X s ) dW s , t ∈ [0 , τ )for Y by approximating the paths of Z = ∇ u ( X ) by a DNN, see Section 3. Doing so, we obtainan estimate of a solution value u ( x ) for a given x ∈ R d . The idea of the proposed algorithm is inspired by [33], where the authors use the correspondencebetween BSDEs and semilinear parabolic PDEs to construct a DNN algorithm for solving thelatter. In the same spirit, we construct a DNN algorithm based on the correspondence to BSDEswith random terminal time for solving semilinear elliptic PDEs.The details of the algorithm are described in three steps of increasing specificity. First we ex-plain the DNN algorithm mathematically. This is done below. Second, Algorithm 1 at the endof this section provides a pseudocode. Third, our program code is provided on Github under acreative commons license. The algorithm is implemented in a generic manner so that it can bereused for other elliptic PDE problems.The goal of the algorithm is to calculate solution values u ( x ) of the semilinear (degenerate)elliptic PDE of interest. For the construction of the algorithm we use the correspondence to aBSDE with random terminal time. Recall from Section 2 that such a BSDE is given by a triplet( τ, g ( X τ ) , f ) that can be determined from the given PDE and by X t = x + (cid:90) t µ ( X s ) ds + (cid:90) t σ ( X s ) dW s (10)and Y t = Y − (cid:90) t f ( X s , Y s , Z s ) ds + (cid:90) t Z s σ ( X s ) dW s . (11)Furthermore, recall that we have identified Y t = u ( X t ), where u is the solution of the PDE weare interested in.The first step for calculating u is to approximate (10) up to the stopping time τ . To make thiscomputationally feasible, we choose T large and stop at τ ∧ T , hence at T at the latest. Now let Since Z σ ( X ) = Z is solvable for Z , f is well-defined. https://github.com/stefankremsner/elliptic-pdes t < t < · · · < t N = T , ∆ t n = t n +1 − t n . We simulate M paths ω , . . . , ω M of the Brownianmotion W . With this we approximate the forward process using the Euler-Maruyama scheme,that is X = x and X t n +1 ≈ X t n + µ ( X t n )∆ t n + σ ( X t n )∆ W n . (12)In the next step we compute Z . For all t n , Z t n = ∇ u ( X t n ) are approximated by DNNs, eachmapping G to R d . As noted above, the implementation of this (and all other steps) is provided.Now, we initialize u ( x ) and use the above approximations to compute the solution to the BSDEforward in time by approximating (11): u ( X t n +1 ) ≈ u ( X t n ) − (0 ,τ ) ( t n ) f ( X t n , u ( X t n ) , ∇ u ( X t n )) ∆ t n + (0 ,τ ) ( t n ) ∇ u ( X t n ) σ ( X t n )∆ W n . (13)Note that due to this construction, indirectly u ( X t n +1 ) is also approximated by a DNN as acombination of DNNs.For the training of the involved DNNs, we compare u ( X τ ∧ T ) with the terminal value ξ . Thisdefines the loss function for the training:1 M M (cid:88) k =1 | u ( X τ ∧ T ( ω k )) − ξ ( ω k ) | . After a certain number of training epochs the loss function is minimized and we obtain an ap-proximate solution value u ( x ) of the PDE. Remark 3.1.
Several approximation errors arise in the proposed algorithm:1. the approximation error of the Euler-Maruyama method, which is used for sampling theforward equation,2. the error of approximating the expected loss,3. the error of cutting off the potentially unbounded random terminal time at time T ,4. the approximation error of the deep neural network model for approximating Z t n for each t n .It is well known that for any continuous function we can find DNNs that approximate the functionarbitrarily well, see [38, 37]. This is, however, not sufficient to make any statement about theapproximation quality. Results on convergence rates are required. Though this question isalready studied in the literature (see, e.g., [8, 9, 13, 24, 27, 30, 31, 32, 42, 40, 41, 45, 48, 60]),results on convergence rates for given constructions are yet scarce and hence many questionsremain open while the number of proposed DNN algorithms grows.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 ) iscalculated as the mean value of u ( x ) in the last 3 network training epochs of each run,sampled according to the validation size (see below).7 r b N T E M validation size time per eight points .
75 500 0.5 200 64 256 119.17 s100 0.5 0 .
75 500 0.01 200 64 256 613.86 sTable 1: Parameters for the Poisson equation. • We choose a non-equidistant time grid in order to get a higher resolution for earlier (andhence 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. In this section we apply the proposed algorithm to three examples. The first one serves as avalidity check, the second one as an academic example with a non-linearity. Finally, we applythe algorithm to solve the dividend maximization problem under incomplete information.
The first example we study is the Poisson equation – a linear PDE.Let r ∈ (0 , ∞ ), G = (cid:8) x ∈ R d : | x | < r (cid:9) , ∂G = (cid:8) x ∈ R d : | x | = r (cid:9) , b ∈ R , and∆ u ( x ) = − b, x ∈ G,u ( x ) = 0 , x ∈ ∂G. (14)Solving (14) is equivalent to solving the BSDE with dX t = √ dW t , X = x,f ( x, y, ζ ) = b, ξ = 0 , up to the stopping time τ = inf { t ∈ [0 , T ] : | x | > r } .To obtain a reference solution for this linear BSDE we use an analytic formula for the expectationof τ , see [54, Example 7.4.2, p. 121]. This yields u ( x ) = b d (cid:0) r − | x | (cid:1) . We compute u ( x ) on the R and the R for 15 different values of x . Figure 1 shows theapproximate solution of u obtained by the DNN algorithm on the diagonal points { ( x, . . . , x ) ∈ R d : x ∈ [ − r, r ] } (in blue) and the analytical reference solution (in green). Table 1 contains theparameters we use.Note that as the expected value of τ decreases linearly in d , we adapt the cut off time T for d = 100 accordingly. The numerical examples were run on a Lenovo Thinkpad notebook with an Intel Core i7 processor (2.6 GHz)and 16 GB memory, without CUDA. lgorithm 1 Elliptic PDE Solver for a BSDE ( f, ξ ) with stopping time τ Require: number of training epochs E , maximal time T , step-size ∆ t , number of timesteps N , numberof sample paths M , number of hidden layer neurons dim, initial (random) starting values ( θ ( u )0 , θ ( ζ )0 ) function TrainableVariables (dim , θ ) (cid:46) see Pytorch or Tensorflow return a trainable variable with dimension 1 × dim initialized by θ . end function function Subnetwork ( x ) (cid:46) allowing x to be a tensor containing M rows (samples) return a trainable DNN, evaluated at x . end function for i = 0 , . . . , N do t i = timesteps( i ) (cid:46) Initialize non-equidistant timesteps end for for j = 1 , . . . , M do Sample Brownian motion trajectory (cid:16) w ( j ) t i (cid:17) ≤ i ≤ N Sample path from forward process (cid:16) x ( j ) t i (cid:17) ≤ i ≤ N calculate stopping time τ ( j ) calculate terminal value ξ ( j ) set x ( j ) t = x ( j ) τ ( j ) for all t > τ ( j ) end for u = TrainableVariables (1 , θ ( u )0 ) (cid:46) Initialize u ∇ u = TrainableVariables ( d, θ ( ζ )0 ) (cid:46) Initialize Z for j = 1 , . . . , M do u ( j ) = u ∇ u ( j ) = ∇ u end for for e = 1 , . . . , E do for i = 1 , . . . , N − do for j = 1 , . . . , M do u ( j ) = u ( j ) − f ( x ( j ) t i , u ( j ) , ∇ u ( j ) )( t i +1 − t i ) + ∇ u ( j ) σ ( x ( j ) t i )( w ( j ) t i +1 − w ( j ) t i ) if t i +1 > τ ( j ) then break end if end for ∇ u = Subnetwork ( x t i +1 ) end for update all trainable variables and the subnetwork’s weights according to the loss function1 M M (cid:88) j =1 ( u ( j ) − ξ ( j ) ) end forreturn ( u , ∇ u ) R (left) and on the R (right). d r N T E M validation size time per eight points2 1 100 5 500 64 256 204.58 s100 1 100 0.1 500 64 256 321.13 sTable 2: Parameters for the equation with quadratic gradient.
The second example is a semilinear PDE with a quadratic gradient term.Let r ∈ (0 , ∞ ), G = (cid:8) x ∈ R d : | x | < r (cid:9) , and ∂G = (cid:8) x ∈ R d : | x | = r (cid:9) . We consider the PDE∆ u ( x ) + |∇ u ( x ) | = 2 e − u ( x ) , x ∈ G,u ( x ) = log (cid:18) r + 1 d (cid:19) , x ∈ ∂G. (15)corresponding to the BSDE dX t = √ dW t , X = x, (16) f ( x, y, ζ ) = | ζ | − e − y , ξ = log (cid:18) | r | + 1 d (cid:19) . (17)In addition, this example we have an analytic reference solution given by u ( x ) = log (cid:18) | x | + 1 d (cid:19) . As in the previous example we compute u ( x ) for 15 different values of x on the R and the R .Figure 2 shows the approximate solution of u obtained by the DNN algorithm on the diagonalpoints { ( x, . . . , x ) ∈ R d : x ∈ [ − r, r ] } (in blue) and the analytical reference solution (in green).Table 2 contains the parameters we use.While classical numerical methods for PDEs would be a much better choice in the case d = 2,their application would not be feasible in the case d = 100.10igure 2: Approximate solution (blue) and reference solution (green) for the equation withquadratic gradient on the R (left) and on the R (right). The goal of this paper was to show how to use the proposed DNN algorithm to solve high-dimensional control problems that arise in insurance mathematics. We finally arrived at thepoint where we are ready to do so.Our example comes from [68], where the author studies De Finetti’s dividend maximizationproblem in a setup with incomplete information about the current state of the market. The hid-den market-state process determines the trend of the surplus process of the insurance companyand is modeled as a d -state Markov chain. Using stochastic filtering, in [68] they achieve totransform the one-dimensional problem under incomplete information to a d -dimensional prob-lem 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 [68] in orderto avoid ambiguities.The probability that the Markov chain modeling the market-state is in state i ∈ { , . . . , d − } is given by π i ( t ) = x i + (cid:90) t q d,i + d − (cid:88) j =1 ( q j,i − q d,i ) π j ( s ) ds + (cid:90) t π i ( s ) a i − ν s ρ dB s , (18)where ν t = a d + d − (cid:88) j =1 ( a j − a d ) π j ( t ) , (19) x i ∈ (0 , B is a one-dimensional Brownian motion, a , . . . , a d ∈ R are the values of the surplustrend in the respective market-states of the hidden Markov chain, and ( q i,j ) i,j ∈{ ,...,d } ∈ R d × d denotes the intensity matrix of the chain.Let ( (cid:96) t ) t ∈ [0 , ∞ ) be the dividend rate process. The surplus of the insurance company is given by (cid:101) X dt = x d + (cid:90) t ( ν s − (cid:96) s ) ds + ρB t , t ∈ [0 , ∞ ) , (20)where x d , ρ ∈ (0 , ∞ ). For later use we define also the dividend free surplus X dt = x d + (cid:90) t ν s ds + ρB t , t ∈ [0 , ∞ ) . (21)11he processes (18) and (20) form the d -dimensional state process underlying the optimal controlproblem we aim to solve.The goal of the insurance company is to determine its value by maximizing the discounteddividends payments until the time of ruin η = inf { t ∈ (0 , ∞ ] : (cid:101) X dt < } , that is it seeks to find u ( x , . . . , x d ) = sup ( (cid:96) t ) t ∈ [0 , ∞ ) ∈ A E x ,...,x d (cid:20)(cid:90) η e − δt (cid:96) t dt (cid:21) , (22)where δ ∈ (0 , ∞ ) is a discount rate, A is the set of admissible controls, and E x ,...,x d [ · ] denotesthe expectation under the initial conditions π i (0) = x i for i ∈ { , . . . , d − } and (cid:101) X d = x d .Admissible controls are ( F X d t ) t ≥ -progressively measurable, [0 , K ]-valued for K ∈ (0 , ∞ ), andfulfill (cid:96) t ≡ t > η , cf. [68].In order to tackle the problem, we solve the corresponding Hamilton-Jacobi-Bellmann (HJB)equation from [68], ( L − δ ) u + sup (cid:96) ∈ [0 ,K ] ( (cid:96) (1 − u x d )) = 0 , (23)where L is the second order degenerate elliptic operator L u = a d u x d + d − (cid:88) i =1 ( a i − a d ) x i u x d + q di + d − (cid:88) j =1 ( q ji − q di ) x i u x i + x i ( a i − ν ) u x d x i + 12 d − (cid:88) j =1 (cid:18)(cid:18) x i a i − νρ (cid:19) (cid:18) x j a j − νρ (cid:19) u x i x j (cid:19) + 12 ρ u x d x d . The supremum in (23) is attained at (cid:96) = (cid:40) K, u x d ≤ , u x d > . Plugging this into (23) we end up with a d -dimensional semilinear degenerate elliptic PDE:( L − δ ) u + K (1 − u x d ) { u xd ≤ } = 0 . (24)The boundary conditions in x d direction are given by u ( x , . . . , x d ) = (cid:40) K/δ, x d → ∞ , x d = 0 . No boundary conditions are required for the other variables, cf. [68].In [68, Corollary 3.6] it is proven that the unique viscosity solution to (24) solves the optimalcontrol problem (22). Hence, we can indeed solve the control problem by solving the HJB equa-tion.For the numerical approximation we cut off x d at r ∈ (0 , ∞ ). Hence, G = (cid:8) x ∈ R d : 0 < x d < r (cid:9) and ∂G = (cid:8) x ∈ R d : x d ∈ { , r } (cid:9) . For abbreviation we use u x d for ∂u∂x d etc. dX t = ( dπ ( t ) , . . . , dπ d − ( t ) , dX dt ) (cid:62) , X = x, that is dX t = µ ( X t ) dt + σ ( X t ) dW t , X = x, where W = ( B, W , . . . , W d ) (cid:62) , x = ( x , . . . , x d ) (cid:62) , and µ ( x ) = q d, + d − (cid:88) j =1 ( q j, − q d, ) x j , . . . , q d,d − + d − (cid:88) j =1 ( q j,d − − q d,d − ) x d − , a d + d − (cid:88) j =1 ( a j − a d ) x j (cid:62) ,σ ( x ) = x a − a d + (cid:80) d − j =1 ( a j − a d ) x j ρ . . . . . . . . . x d − a d − − a d + (cid:80) d − j =1 ( a j − a d ) x j ρ . . . ρ . . . . We claim that the BSDE associated to (24) is given in forward form by u ( X t ) = u ( x ) − (cid:90) t [ K (1 − u x d ( X s )) { u xd ( X s ) ≤ } − δu ( X s )] dt + (cid:90) t ∇ u ( X s ) σ ( X s ) dW s . (25)Applying Itˆo’s formula to u ( X ) yields u ( X t ) = u ( x ) + (cid:90) t L u ( X s ) dt + (cid:90) t ∇ u ( X s ) σ ( X s ) dW s . (26)Combining (25) and (26) gives u ( x ) − (cid:90) t [ K (1 − u x d ( X s )) { u xd ( X s ) ≤ } − δu ( X s )] dt + (cid:90) t ∇ u ( X s ) σ ( X s ) dW s = u ( x ) + (cid:90) t L u ( X s ) dt + (cid:90) t ∇ u ( X s ) σ ( X s ) dW s . Canceling terms verifies (in a heuristic manner) (24).Hence, the corresponding BSDE has the parameters f ( x, y, ζ ) = K (1 − ζ d ) { ζ d ≤ } − δy (27)and ξ = (cid:40) K/δ, X dτ = r, , X dτ = 0 , if τ < ∞ . 13 r K δ ρ a i N T
E M validation size time per eight points2 5 1.8 0.5 1 (cid:0) − id (cid:1)
100 5 500 64 256 317.42 s100 5 1.8 0.5 1 (cid:0) − id (cid:1)
100 5 500 64 256 613.15 sTable 3: Parameters for the dividend problem.case i = j even i = j odd i = j + 1 even i = j + 1 ≥ i = 1, j = d otherwise q i,j − . − .
25 0 . As for this example we have no analytic reference solution at hand, we use the solution from[68] 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 inreasonable computation time.As in the previous examples we compute u ( x ) on the R and on the R for 15 different valuesof x . Figure 3 shows the approximate solution of the HJB equation and hence the value of theinsurance company obtained by the DNN algorithm (in blue) and the reference solution from[68] (in green) for the case d = 2. For d = 100 we have no reference solution at hand. Figure 4shows the loss for a fixed value of x in the case d = 100. Tables 3 and 4 contain the parameterswe use. The goal of this paper was to compute the risk measure given by the expected discounted futuredividend payments in a complex high-dimensional economic environment. This demonstratesthe effectiveness of using DNN algorithms for solving some high-dimensional PDE problems ininsurance mathematics that cannot be solved by classical methods. In the literature the focus sofar was on parabolic PDE problems; however, in insurance mathematics we often face problemsup to an unbounded random terminal time, e.g., the time of ruin of the insurance company,leading to (degenerate) elliptic problems.Hence, we have proposed a novel deep neural network algorithm for a large class of semilinear(degenerate) elliptic PDEs associated to infinite time horizon control problems in high dimen-sions. The method extends the DNN algorithm proposed by Han, Jetzen, and E [33], which wasdeveloped for parabolic PDEs, to the case of (degenerate) elliptic semilinear PDEs. We haveattacked the problem inspired by a series of results by Pardoux [55].Of course, in low dimensions one would not use the proposed DNN algorithm – classical meth-ods are more efficient. However, recent models are frequently high dimensional, in which caseclassical methods fail due to the curse of dimensionality. Then the DNN algorithm presentedhere can be applied to compute the desired quantity.We emphasize that the method presented here can also be applied to many other high-dimensionalsemilinear (degenerate) elliptic PDE problems in insurance mathematics and beyond.An implementation of the algorithm is provided on Github under a creative commons license. https://github.com/stefankremsner/elliptic-pdes R for fixed π = π = 0 . R for fixed π = · · · = π = 0 .
01 (right,without reference solution).Figure 4: Interpolated loss for the case d = 100.15 cknowledgements The authors thank Gunther Leobacher for discussions and suggestions that improved the paper.S. Kremsner is supported by the Austrian Science Fund (FWF): Project F5508-N26, which ispart of the Special Research Program “Quasi-Monte Carlo Methods: Theory and Applications”.
References [1] H. Albrecher and S. Thonhauser. Optimality results for dividend problems in insurance.
RACSAM Revista de la Real Academia de Ciencias Exactas, Fisicas y Naturales. Serie A.Matematicas , 103(2):295–320, 2009.[2] S. Asmussen and M. Taksar. Controlled diffusion models for optimal dividend pay-out.
Insurance: Mathematics and Economics , 20(1):1–15, 1997.[3] B. Avanzi. Strategies for dividend distribution: A review.
North American ActuarialJournal , 13(2):217–251, 2009.[4] P. Azcue and N. Mular.
Stochastic Optimization in Insurance – A Dynamic Program-ming Approach . Springer Briefs in Quantitative Finance. Springer, New York, Heidelberg,Dordrecht, London, 2014.[5] C. Beck, S. Becker, P. Cheridito, A. Jentzen, and A. Neufeld. Deep splitting method forparabolic PDEs. arXiv:1907.03452 , 2019.[6] C. Beck, S. Becker, P. Grohs, N. Jaafari, and A. Jentzen. Solving stochastic differentialequations and Kolmogorov equations by means of deep learning. arXiv:1806.00421 , 2018.[7] C. Beck, W. E, and A. Jentzen. Machine learning approximation algorithms for high-dimensional fully nonlinear partial differential equations and second-order backwardstochastic differential equations.
Journal of Nonlinear Science , 29(4):1563–1619, 2019.[8] C. Beck, L. Gonon, and A. Jentzen. Overcoming the curse of dimensionality in the nu-merical approximation of high-dimensional semilinear elliptic partial differential equations. arXiv:2003.00596 , 2020.[9] C. Beck, F. Hornung, M. Hutzenthaler, A. Jentzen, and T. Kruse. Overcoming the curse ofdimensionality in the numerical approximation of Allen-Cahn partial differential equationsvia truncated full-history recursive multilevel Picard approximations. arXiv:1907.06729 ,2019.[10] S. Becker, P. Cheridito, and A. Jentzen. Deep optimal stopping.
Journal of MachineLearning Research , 20:74, 2019.[11] S. Becker, P. Cheridito, A. Jentzen, and T. Welti. Solving high-dimensional optimal stop-ping problems using deep learning. arXiv:1908.01602 , 2019.[12] J. Berg and K. Nystr¨om. A unified deep artificial neural network approach to partialdifferential equations in complex geometries.
Neurocomputing , 317:28–41, 2018.[13] J. Berner, P. Grohs, and A. Jentzen. Analysis of the generalization error: Empirical riskminimization over deep artificial neural networks overcomes the curse of dimensionality inthe numerical approximation of Black–Scholes partial differential equations.
SIAM Journalon Mathematics of Data Science , 2(3):631–657, 2020.1614] P. Briand, B. Delyon, Y. Hu, ´E. Pardoux, and L. Stoica. L p solutions of backward stochasticdifferential equations. Stochastic Processes and their Applications , 108(1):109–129, 2003.[15] P. Briand and Y. Hu. Stability of BSDEs with random terminal time and homogenizationof semi-linear elliptic PDEs.
Journal of Functional Analysis , 155(2):455–494, 1998.[16] Q. Chan-Wai-Nam, J. Mikael, and X. Warin. Machine learning for semi linear PDEs.
Journal of Scientific Computing , 79(3):1667–1712, 2019.[17] Y. Chen and J. W. Wan. Deep neural network framework based on backward stochasticdifferential equations for pricing and hedging american options in high dimensions.
Quan-titative Finance , pages 1–23, 2020.[18] F. Confortola and P. Briand. Quadratic BSDEs with random terminal time and ellipticPDEs in infinite dimension.
Electronic Journal of Probability , 13:1529–1561, 2008.[19] R. Darling and ´E. Pardoux. Backwards SDE with random terminal time and applicationsto semilinear elliptic PDE.
The Annals of Probability , 25(3):1135–1159, 1997.[20] B. de Finetti. Su un’impostazione alternativa della teoria collettiva del rischio.
Transactionsof the XVth International Congress of Actuaries , 2:433–443, 1957.[21] T. Dockhorn. A discussion on solving partial differential equations using neural networks. arXiv:1904.07200 , 2019.[22] W. E, J. Han, and A. Jentzen. Deep learning-based numerical methods for high-dimensionalparabolic partial differential equations and backward stochastic differential equations.
Com-munications in Mathematics and Statistics , 5(4):349–380, 2017.[23] W. E and B. Yu. The deep Ritz method: A deep learning-based numerical algorithm forsolving variational problems.
Commun. Math. Stat. , 6:1–12, 2018.[24] D. Elbr¨achter, P. Grohs, A. Jentzen, and C. Schwab. DNN expression rate analysis ofhigh-dimensional PDEs: application to option pricing. arXiv:1809.07669 , 2018.[25] A.-M. Farahmand, S. Nabi, and D. Nikovski. Deep reinforcement learning for partial dif-ferential equation control. , pages 3120–3127,2017.[26] M. Fujii, A. Takahashi, and M. Takahashi. Asymptotic expansion as prior knowledge in deeplearning method for high dimensional BSDEs.
Asia-Pacific Financial Markets , 26(3):391–408, 2019.[27] L. Gonon, P. Grohs, A. Jentzen, D. Kofler, and D. ˇSiˇska. Uniform error estimates forartificial neural network approximations. arXiv:1911.09647 , 2019.[28] L. Gouden`ege, A. Molent, and A. Zanette. Machine learning for pricing American optionsin high dimension. arXiv:1903.11275 , 2019.[29] P. Grohs and L. Herrmann. Deep neural network approximation for high-dimensional ellipticPDEs with boundary conditions. arXiv:2007.05384 , 2020.[30] P. Grohs, F. Hornung, A. Jentzen, and P. von Wurstemberger. A proof that artificialneural networks overcome the curse of dimensionality in the numerical approximation ofBlack-Scholes partial differential equations. arXiv:1809.02362 , 2018.1731] P. Grohs, F. Hornung, A. Jentzen, and P. Zimmermann. Space-time error estimates fordeep neural network approximations for differential equations. arXiv:1908.03833 , 2019.[32] P. Grohs, A. Jentzen, and D. Salimova. Deep neural network approximations for MonteCarlo algorithms. arXiv:1908.10828 , 2019.[33] J. Han, A. Jentzen, and W. E. Solving high-dimensional partial differential equations usingdeep learning.
Proceedings of the National Academy of Sciences , 115(34):8505–8510, 2018.[34] J. Han and J. Long. Convergence of the deep BSDE method for coupled FBSDEs.
Proba-bility, Uncertainty and Quantitative Risk , 5(1):1–33, 2020.[35] J. Han, M. Nica, and A. R. Stinchcombe. A derivative-free method for solving ellipticpartial differential equations with deep neural networks. ”Journal of Computational Physics ,419:109672, 2020.[36] P. Henry-Labord`ere. Deep primal-dual algorithm for BSDEs: Applications of machinelearning to CVA and IM.
Available at SSRN 3071506 , 2017.[37] K. Hornik. Approximation capabilities of multilayer feedforward networks.
Neural networks ,4(2):251–257, 1991.[38] K. Hornik, M. Stinchcombe, and H. White. Multilayer feedforward networks are universalapproximators.
Neural networks , 2(5):359–366, 1989.[39] C. Hur´e, H. Pham, and X. Warin. Some machine learning schemes for high-dimensionalnonlinear PDEs. arXiv:1902.01599 , 2019.[40] M. Hutzenthaler, A. Jentzen, T. Kruse, T. A. Nguyen, and P. von Wurstemberger. Over-coming the curse of dimensionality in the numerical approximation of semilinear parabolicpartial differential equations. arXiv:1807.01212 , 2018.[41] M. Hutzenthaler, A. Jentzen, and P. von Wurstemberger. Overcoming the curse of dimen-sionality in the approximative pricing of financial derivatives with default risks.
ElectronicJournal of Probability , 25:73 pp, 2020.[42] Martin Hutzenthaler, Arnulf Jentzen, Thomas Kruse, and Tuan Anh Nguyen. A proof thatrectified deep neural networks overcome the curse of dimensionality in the numerical approx-imation of semilinear heat equations.
SN Partial Differential Equations and Applications ,1:1–34, 2020.[43] A. Jacquier and M. Oumgari. Deep PPDEs for rough local stochastic volatility. arXiv:1906.02551 , 2019.[44] M. Jeanblanc-Piqu´e and A. N. Shiryaev. Optimization of the flow of dividends.
RussianMathematical Surveys , 50(2):257–277, 1995.[45] A. Jentzen, D. Salimova, and T. Welti. A proof that deep artificial neural networks overcomethe curse of dimensionality in the numerical approximation of Kolmogorov partial differ-ential equations with constant diffusion and nonlinear drift coefficients. arXiv:1809.07321 ,2018.[46] Z. Jiang and M. Pistorius. Optimal dividend distribution under markov regime switching.
Finance and Stochastics , 16(3):449–476, 2012.1847] P. Kritzer, G. Leobacher, M. Sz¨olgyenyi, and S. Thonhauser. Approximation methods forpiecewise deterministic markov processes and their costs.
Scandinavian Actuarial Journal ,2019(4):308–335, 2019.[48] G. Kutyniok, P. Petersen, M. Raslan, and R. Schneider. A theoretical analysis of deepneural networks and parametric PDEs. arXiv:1904.00377 , 2019.[49] G. Leobacher, M. Sz¨olgyenyi, and S. Thonhauser. Bayesian dividend optimization and finitetime ruin probabilities.
Stochastic Models , 30(2):216–249, 2014.[50] Z. Long, Y. Lu, X. Ma, and B. Dong. PDE-Net: Learning PDEs from Data.
In Proceedingsof the 35th International Conference on Machine Learning , pages 3208–3216, 2018.[51] L. Lu, X. Meng, Z. Mao, and G. E. Karniadakis. DeepXDE: A deep learning library forsolving differential equations. arXiv:1907.04502 , 2019.[52] K. O. Lye, S. Mishra, and D. Ray. Deep learning observables in computational fluid dy-namics.
Journal of Computational Physics , 410:109339, 2020.[53] M. Magill, F. Qureshi, and H. de Haan. Neural networks trained to solve differential equa-tions learn general representations.
Advances in Neural Information Processing Systems ,pages 4071–4081, 2018.[54] B. Øksendal.
Stochastic differential equations . Springer, 2003.[55] ´E. Pardoux. Backward stochastic differential equations and viscosity solutions of systems ofsemilinear parabolic and elliptic PDEs of second order. In
Stochastic Analysis and RelatedTopics VI , pages 79–127. Springer, 1998.[56] ´E. Pardoux. BSDEs, weak convergence and homogenization of semilinear PDEs. In
Non-linear analysis, differential equations and control , pages 503–549. Springer, 1999.[57] H. Pham and X. Warin. Neural networks-based backward scheme for fully nonlinear PDEs. arXiv:1908.00412 , 2019.[58] R. Radner and L. Shepp. Risk vs. profit potential: A model for corporate strategy.
Journalof Economic Dynamics and Control , 20(8):1373–1393, 1996.[59] M. Raissi. Deep hidden physics models: Deep learning of nonlinear partial differentialequations.
J. Mach. Learn. Res. , 19(1):932–955, 2018.[60] C. Reisinger and Y. Zhang. Rectified deep neural networks overcome the curse of di-mensionality for nonsmooth value functions in zero-sum games of nonlinear stiff systems. arXiv:1903.06652 , 2019.[61] A. M. Reppen, J.-C. Rochet, and H. M. Soner. Optimal dividend policies with randomprofitability.
Mathematical Finance , 30(1):228–259, 2020.[62] M. Royer. BSDEs with a random terminal time driven by a monotone generator and theirlinks with PDEs.
Stochastics and stochastic reports , 76(4):281–307, 2004.[63] H. Schmidli.
Stochastic Control in Insurance . Probability and its Applications. Springer,London, 2008.[64] S. E. Shreve, J. P. Lehoczky, and D. P. Gaver. Optimal consumption for general diffu-sions with absorbing and reflecting barriers.
SIAM Journal on Control and Optimization ,22(1):55–75, 1984. 1965] J. Sirignano and K. Spiliopoulos. DGM: A deep learning algorithm for solving partialdifferential equations.
Journal of computational physics , 375:1339–1364, 2018.[66] L. Sotomayor and A. Cadenillas. Classical and singular stochastic control for the optimaldividend policy when there is regime switching.
Insurance: Mathematics and Economics ,48:344–354, 2011.[67] M. Sz¨olgyenyi. Bayesian dividend maximization: A jump diffusion model. In M. Vanmaele,G. Deelstra, A. De Schepper, J. Dhaene, W. Schoutens, S. Vanduffel, and D. Vyncke, ed-itors,
Handelingen Contactforum Actuarial and Financial Mathematics Conference, Inter-play between Finance and Insurance, February 7-8, 2013 , pages 77–82. Koninklijke VlaamseAcademie van Belgi¨e voor Wetenschappen en Kunsten, Brussel, 2013.[68] M. Sz¨olgyenyi. Dividend maximization in a hidden Markov switching model.
Statistics &Risk Modeling , 32(3-4):143–158, 2016.[69] J. Zhu and F. Chen. Dividend optimization for regime-switching general diffusions.