A numerical algorithm for fully nonlinear HJB equations:an approach by control randomization

# A numerical algorithm for fully nonlinear HJB equations: an approach by control randomization

Idris Kharroubi111The research of the author benefited from the support of the French ANR research grant LIQUIRISK (ANR-11-JS01-0007).
CEREMADE, CNRS UMR 7534,
Université Paris Dauphine
and CREST,
Nicolas Langrené
Laboratoire de Probabilités et Modèles Aléatoires,
Université Paris Diderot
and EDF R&D
langrene at math.univ-paris-diderot.fr
Huyên Pham
Laboratoire de Probabilités et Modèles Aléatoires,
Université Paris Diderot
and CREST-ENSAE
pham at math.univ-paris-diderot.fr
###### Abstract

We propose a probabilistic numerical algorithm to solve Backward Stochastic Differential Equations (BSDEs) with nonnegative jumps, a class of BSDEs introduced in [9] for representing fully nonlinear HJB equations. In particular, this allows us to numerically solve stochastic control problems with controlled volatility, possibly degenerate. Our backward scheme, based on least-squares regressions, takes advantage of high-dimensional properties of Monte-Carlo methods, and also provides a parametric estimate in feedback form for the optimal control. A partial analysis of the error of the scheme is provided, as well as numerical tests on the problem of superreplication of option with uncertain volatilities and/or correlations, including a detailed comparison with the numerical results from the alternative scheme proposed in [7].

Key words: Backward stochastic differential equations, control randomization, HJB equation, uncertain volatility, empirical regressions, Monte-Carlo.

MSC Classification: 60H10, 65Cxx, 93E20.

## 1 Introduction

Consider the following general Hamilton-Jacobi-Bellman (HJB) equation:

 =0,(t,x)∈[0,T)×Rd (1.1) v(T,x) =g(x),x∈Rd

where is a bounded subset of . It is well known that the HJB equation (1.1) is the dynamic programming equation for the following stochastic control problem:

 v(t,x) =supα∈AEt,x[∫Ttf(Xαs,αs)ds+g(XαT)] (1.2) dXαs =b(Xαs,αs)ds+σ(Xαs,αs)dWs

Moreover, it is proved in [9] that this HJB equation admits a probabilistic representation by means of a BSDE with nonpositive jumps. We recall below this construction.

Introduce a Poisson random measure on with finite intensity measure associated to the marked point process , independent of , and consider the pure jump process , valued in , defined as follows:

 It=ζi,τi≤t<τi+1,

and interpreted as a randomization of the control process .

Next, consider the uncontrolled forward regime switching diffusion process

 dXs=b(Xs,Is)ds+σ(Xs,Is)dWs.

Observe that the pair process is Markov. Now, consider the following BSDE with jumps w.r.t. the Brownian-Poisson filtration .

 Yt=g(XT)+∫Ttf(Xs,Is,Ys,Zs)ds−∫TtZsdWs−∫Tt∫AUs(a)~μA(ds,da) (1.3)

where is the compensated measure of .

Finally, we constrain the jump component of the BSDE (1.3) to be nonpositive, i.e.

 Ut(a)≤0,dP⊗dt⊗λ(da)a.e.

We denote by an upper bound for the compact set of , i.e. for all , and we make the standing assumptions:

1. The functions and are Lipschitz: there exists s.t.

 |b(x1,a1)−b(x2,a2)|+|σ(x1,a1)−σ(x2,a2)| ≤ Lb,σ(|x1−x2|+|a1−a2|),

for all , .

2. The functions and are Lipschitz continuous: there exists s.t.

 |g(x1)−g(x2)| ≤ Lg|x1−x2| |f(x1,a1,y1,z1)−f(x2,a2,y2,z2)| ≤ Lf(|x1−x2|+|a1−a2|+|y1−y2|+|z1−z2|),

for all , .

Under these conditions, we consider the minimal solution of the following constrained BSDE:

 Yt = g(XT)+∫Ttf(Xs,Is,Ys,Zs)ds−∫TtZsdWs +KT−Kt−∫Tt∫AUs(a)~μA(ds,da),0≤t≤T,a.s.

subject to the constraint

 Ut(a)≤0,dP⊗dt⊗λ(da)a.e.onΩ×[0,T]×A (1.5)

By the Markov property of , there exists a deterministic function such that the minimal solution to (1)-(1.5) satisfies , .

###### Theorem 1.1.

[9] does not depend on : , and is a viscosity solution of the HJB equation (1.1):

 ∂y∂t+supa∈A{b(x,a).Dxy(t,x)+12tr(σσ⊤(x,a)D2xv(t,x))+f(x,a,y,σ⊤(x,a)Dxy)} =0 =g(x),x∈Rd

Now, the aim of this paper is to provide a numerical scheme for computing an approximation of the solution of the constrained BSDE (1)-(1.5). In light of Theorem 1.1, this will provide an approximation of the solution of the general HJB equation (1.1), which encompasses stochastic control problems such as the one described in equation (1.2), ie. problems where both the drift and the volatility of the underlying diffusion can be controlled, including degenerate diffusion coefficient.

The outline of the subsequent sections is the following.

First, Section 2 describes our scheme. We start from a time-discretization of the problem, proposed in [8], which gives rise to a backward scheme involving the simulation of the forward regime switching process , hence taking advantage of high-dimensional properties of Monte-Carlo methods. The final step towards an implementable scheme is to approximate the conditional expectations that arise from this scheme. Here we use empirical least-squares regression, as this method provides a parametric estimate in feedback form of the optimal control. A partial analysis of the impact of this approximation is provided, and the remaining obstacles towards a full analysis are highlighted.

Then, Section 3 is devoted to numerical tests of the scheme on various examples. The major application that takes advantage of the possibilities of our scheme is the problem of pricing and hedging contingent claims under uncertain volatility and (for multi-dimensional claims) correlation. Therefore most of this section is devoted to this specific application. To our knowledge, the only other Monte Carlo scheme for HJB equations that can handle continuous controls as well as controlled volatility is described in [7], where they make use of another generalization of BSDEs, namely second-order BSDEs. Therefore we compare the performance of our scheme to the results provided in their paper.

Finally, Section 4 concludes the paper.

## 2 Regression scheme

Define a deterministic time grid := for the interval , with mesh where . Denote by . The discretization of the constrained BSDE (1)-(1.5) can be written as follows:

 ⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩YN=g(XN)ΔiZi=Ei[Yi+1ΔW⊤i]Yi=Ei[Yi+1+f(Xi,Ii,Yi+1,Zi)Δi]Yi=esssupa∈AEi,a[Yi] (2.1)

where .

First, remark that, from the Markov property of , there exist deterministic functions and such that . Hence and can be seen as intermediate quantities towards the discrete-time approximation of the BSDE (1)-(1.5) , which do not depend on .

Formally, the jump constraint (1.5) states that a.s., meaning that the minimal solution satisfies .

Moreover, one can extract from the scheme if needed. Indeed, denoting , i.e. , then .

Finally remark that the numerical scheme (2.1) is explicit, as we choose to define as a function of and not of .

The convergence of the solution of the discretized scheme (2.1) towards the solution of the constrained BSDE (1)-(1.5) is thoroughly examined in [8]. In this paper, we start from the discrete version (2.1) and derive an implementable scheme from it.

Indeed, the discrete scheme (2.1) is in general not readily implementable because it involves conditional expectations that cannot be computed explicitly. It is thus necessary in practice to approximate these conditional expectations. Here we follow the empirical regression approach ([10, 3, 6, 18, 1]). In our context, apart from being easy to implement, the strong advantage of this choice is that, unlike other standard methods, it provides as a by-product a parametric feedback estimate for the optimal control.

The idea is to replace the conditional expectations from (2.1) by empirical regressions. This section is devoted to the analysis of the error generated by this replacement.

### 2.1 Localizations

The first step is to localize the discrete BSDE (2.1), i.e. to truncate it so that it admits a.s. deterministic bounds. Introduce and and define the following truncations of and :

 [Xi]X :=−RX∨Xi∧RX={−R1,X∨X1,i∧R1,X,…,−Rd,X∨Xd,i∧Rd,X}⊤ (2.2) [ΔWi]w :=−Rw√Δi∨ΔWi∧Rw√Δi={−Rw√Δi∨ΔW1,i∧Rw√Δi,…,−Rw√Δi∨ΔWq,i∧Rw√Δi}⊤ (2.3)

Define and define the localized version of the discrete BSDE (2.1), using the truncations (2.2) and (2.3).

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩YRN=g([XN]X)ΔiZRi=Ei[YRi+1[ΔW⊤i]w]YRi=Ei[YRi+1+f([Xi]X,Ii,YRi+1,ZRi)Δi]YRi=esssupa∈AEi,a[YRi] (2.4)

First, we check that this localized BSDE does admit a.s. bounds.

###### Lemma 2.1.

[almost sure bounds] For every and every , the following uniform bounds hold a.s.:

 ∣∣YRi∣∣,∣∣YRi∣∣ ≤ Cy=Cy(RX):=eC2T ⎷C2g(RX)+eC|π|L2fC2f(RX) ∣∣ZRi∣∣,∣∣ZRi∣∣ ≤ Cz=Cy(RX):=√q√ΔiCy

where , and

###### Proof.

First, as is continuous, there exists such that for all , . Hence

 (YRN)2=g2([XN]X)≤C2g(RX) (2.5)

Next,

 ΔiZRi=Ei[YRi+1[ΔWi]w]=Ei[(YRi+1−Ei[YRi+1])[ΔWi]w]

Thus, using the Cauchy-Schwarz inequality and dividing by :

 Δi(ZRi)2≤q(Ei[(YRi+1)2]−Ei[YRi+1]2)

Now, using Young’s inequality with :

Remark that

 ∣∣f([Xi]X,Ii,YRi+1,ZRi)∣∣ ≤ ∣∣f([Xi]X,Ii,YRi+1,ZRi)−f(0,0,0,0)∣∣+|f(0,0,0,0)| ≤ ≤ Cf(RX)+Lf(∣∣YRi+1∣∣+∣∣ZRi∣∣)

where . Hence

 (YRi)2 ≤ (1+γΔi)Ei[YRi+1]2+3(Δi+1γ)Δi(C2f(RX)+L2fEi[(YRi+1)2]+L2fEi[∣∣ZRi∣∣]) ≤

Thus, for every , one can group together the terms involving and using Jensen’s inequality:

 (YRi)2 ≤ (1+θ(3,γ)Δi)Ei[(YRi+1)2]+3(|π|+1γ)C2f(RX)Δi

where . Finally:

 (YRi)2≤esssupa∈AEi,a[(YRi+1)2](1+θ(3,γ)Δi)+3(|π|+1γ)C2f(RX)Δi (2.6)

Using equations (2.5) and (2.6), one obtains by induction that:

 (YRi)2≤ΓN−1i(3,γ)C2g(RX)+3(|π|+1γ)C2f(RX)N−1∑k=iΓki(3,γ)Δk (2.7)

where . Finally remark that

 ln(Γji(c,γ))=j∑k=iln(1+θ(c,γ)Δk)≤j∑k=iθ(c,γ)Δk=θ(c,γ)(tj+1−ti)

Thus

 Γji(c,γ)≤exp(θ(c,γ)(tj+1−ti)) (2.8)

And

 N−1∑k=iΓki(c,γ)Δk ≤ N−1∑k=ieθ(c,γ)(tj+1−ti)Δk≤eθ(c,γ)|π|N−1∑k=ieθ(c,γ)(tj−ti)Δk (2.9) ≤ eθ(c,γ)|π|∫tNtieθ(c,γ)(t−ti)dt=eθ(c,γ)|π|θ(c,γ)(eθ(c,γ)(tN−ti)−1)

Finally, combine equations (2.7), (2.8) and (2.9) with and to obtain the following a.s. bound for :

 (YRi)2 ≤ eθ(3,γ)(tN−ti)C2g(RX)+3(|π|+1γ)C2f(RX)eθ(3,γ)|π|θ(3,γ)(eθ(3,γ)(tN−ti)−1) ≤ eθ(3,γ)T{C2g(RX)+3eθ(3,γ)|π|θ(3,γ)(|π|+1γ)C2f(RX)}

In particular, for and :

 (YRi)2 ≤ eCT⎧⎨⎩C2g(RX)+eC|π|L2fC2f(RX)⎫⎬⎭=:C2y

where . The same inequality holds for . For , use the Cauchy-Schwarz inequality to obtain:

 (ZRi)2≤qΔiEi[(YRi+1)2]≤qΔiC2y=:C2z

and the same inequality holds for .∎

###### Lemma 2.2.

For , define where is a Gaussian random variable with mean and variance . Then:

 TR≤√2π1Re−R22
###### Proof.

Developing the square yields

 TR=2R2P(N>R)−4RE[N1{N>R}]+2E[N21{N>R}]

Then the two expectations can be explicited as follows

 E[N1{N>R}] = e−R22√2π E[N21{N>R}] = R√2πe−R22+P(N>R)

Finally, the use of Mill’s ratio inequality concludes the proof. ∎

Then, we can estimate bounds between the BSDEs (2.1) and (2.4).

###### Proposition 2.1.

The following bounds hold:

 (Yi−YRi)2 ≤ eCT{Lg(|ΔXN|2)∗+CN−1∑k=iΔk(|ΔXk|2)∗+2qCTC2yTRw}

where , and , , is the solution of the following linear constrained BSDE:

 ⎧⎨⎩Yk=(Xk−[Xk]X)2Yj=esssupa∈AEj,a[Yj+1],j=k−1,…,i
###### Proof.

Define , , , and . First

 |ΔYN|=∣∣g(XN)−g([XN]X)∣∣≤Lg∣∣ΔXpN∣∣

Then

 ΔiΔZi = Ei[Yi+1ΔW⊤i−YRi+1[ΔW⊤i]w] = Ei[ΔYi+1ΔW⊤i+YRi+1{ΔWi−[ΔWi]w}⊤] Ei[(ΔYi+1−Ei[ΔYi+1])ΔW⊤i]+Ei[YRi+1{ΔWi−[ΔWi]w}⊤]

Hence

 Δi(ΔZi)2≤2q(Ei[(ΔYi+1)2]−Ei[ΔYi+1]2)+2qC2yTRw

Then

 ΔYi=Ei[ΔYi+1+{f(Xi,Ii,Yi+1,Zi)−f([Xi]X,Ii,YRi+1,ZRi)}Δi]

Using Jensen’s inequality and Young’s inequality with parameter , :

 (ΔYi)2 ≤ (1+γΔi)Ei[ΔYi+1]2+(1+1γΔi)Δ2i3L2fEi[(ΔXi)2+(ΔYi+1)2+(ΔZi)2] ≤ Ei[ΔYi+1]2(Δi+1γ)(γ−6qL2f)+Ei[(ΔYi+1)2](Δi+1γ)3L2f(Δi+2q) +(Δi+1γ)Δi3L2f{(ΔXi)2+2qC2yTRw}

Now, for any , one can group together the terms in and using Jensen’s inequality:

 (ΔYi)2≤Ei[(ΔYi+1)2]{1+θ(3,γ)Δi}+3L2f(|π|+1γ)Δi{(ΔXi)2+2qC2yTRw}

where, as in Lemma 2.1, . Hence, using that for any random variables and ,

the following holds:

 (ΔYi)2≤{1+θ(3,γ)Δi}esssupa∈AEi,a[(ΔYi+1)2]+3L2f(|π|+1γ)Δi{(ΔXi)2+2qC2yTRw}

By induction

 (ΔYi)2 ≤ LgΓN−1i(3,γ)(|ΔXN|2)∗ +3L2f(|π|+1γ)N−1∑k=iΔkΓki(3,γ){(|ΔXk|2)∗+2qC2yTRw}

where, as in Lemma 2.1, . Finally, take to obtain the desired bound. ∎

### 2.2 Projections

In its current form, the scheme (2.4) is not readily implementable, because its conditional expectations cannot be computed in general. Therefore, there is a need to approximate these conditional expectations. For handiness and efficiency, we choose, in the spirit of [10] and [6], to approximate them by empirical least-squares regression.

First, we will study the impact of the replacement of the conditional expectations by theoretical least-squares regressions. We will see that the resulting scheme is not easy to analyze. Therefore, we will study a stronger version of it, and discuss their practical differences. As it is already a daunting task for standard BSDEs (cf. [10]), and in view of the difficulties already raised at theoretical regression level, we leave the study of the final replacement of these regressions by their empirical counterparts for further research.

Hence, for each , consider and that are non-empty closed convex subsets of , as well as the corresponding projection operators and . Using the above projection operators in lieu of the conditional expectations in (2.4), we obtain the following approximation scheme:

 (2.10)

where and are truncation operators that ensure that the a.s. upper bounds for from Lemma 2.1 will also hold for .

To be more specific, choose the subsets and as follows:

 SYi ={λ.pYi(Xi,Ii);λ∈RBYi} SZ,ki ={λ.pZ,ki(Xi,Ii);λ∈RBZ,ki},k=1,…,q

where , , and , , are predefined sets of deterministic functions from into . Hence, for any random variable in , is defined as follows:

 ^λYi(U) :=arginfλ∈RBYiE[(λ.pYi(Xi,Ii)−U)2] (2.11) PYi(U) :=^λYi(U).pY