Branching diffusions with jumps, and valuation with systemic counterparties

Christoph Belak; Daniel Hoffmann; Frank Seifried

1 Introduction

The objective of this paper is to derive probabilistic representations of solutions of a certain class of nonlinear parabolic partial differential equations (PDEs) with nonlocal terms in the nonlinearity. The representation is based on a branching diffusion mechanism with jumps at branching times and makes it possible to compute solutions by direct (nonnested) Monte Carlo simulation, leading to a numerical algorithm that does not suffer from the curse of dimensionality.

The class of PDEs under consideration in this paper takes the form

\partial_{t}u(t,x)+\mathcal{A}[u](t,x)+\int_{\varXi}f(t,x,\xi,\mathcal{J}[u](t% ,x,\xi))\gamma(\mathrm{d}\xi)=0,

(1.1)

where $\mathcal{A}$ denotes the infinitesimal generator of an Itô diffusion, ie, a (possibly degenerate) linear partial differential operator of second order; $\mathcal{J}$ is a nonlocal operator; and the nonlinearity $f$ is analytic in the jump terms.

In recent years, starting with Henry-Labordère (2012), there has been significant progress in the realm of probabilistic representations of PDEs with analytic nonlinearities acting on zeroth- and first-order derivatives. We refer in particular to Henry-Labordère et al (2014) for the zeroth-order case and to Henry-Labordère et al (2019) for the first-order case. We also refer to Bouchard et al (2017) for an extension of the branching diffusion approach to the case of locally analytic nonlinearities, to Bouchard et al (2019) for the case of Lipschitz nonlinearities, to Henry-Labordère and Touzi (2018) for higher-order PDEs and to Agarwal and Claisse (2020) for an extension to elliptic equations. The main contribution of the present paper is to extend the branching diffusion approach to the case of nonlocal terms inside the nonlinearity. This extension is achieved by introducing jump marks in the branching diffusion underlying the probabilistic representation result: we consider a branching diffusion similar to the one introduced in Henry-Labordère et al (2014) with the additional feature that, at each branching time, a subset of offspring particles jumps away from their parent’s position. We refer to the resulting object as a branching diffusion with jumps.¹¹ 1 This is not to be confused with the straightforward notion of branching jump-diffusions, where jump terms appear in the infinitesimal generator $\mathcal{A}$ ; see also the discussion below Theorem 2.6.

Our main motivation behind the derivation of a probabilistic representation is to open up the possibility of applying numerical algorithms for nonlocal nonlinear PDEs in high dimensions via Monte Carlo simulation. The effectiveness and efficiency of such algorithms is showcased in an example of a nonlocal nonlinear PDE, which we solve in dimensions up to $100$ . In addition, we show how our Monte Carlo methodology can be used in the pricing of (or equivalently, the computation of credit valuation adjustments for ) financial positions where the counterparty is a systemically important financial institution whose default causes a devaluation of the underlying. This jump-at-default model represents a particularly realistic setup for wrong-way risk (see, for example, Brigo et al 2019; Mercurio and Li 2015; Pykhtin and Sokol 2013). Our Monte Carlo approach complements existing methods by making it possible to price systemic defaultable positions for settings where the underlying dynamics are not tractable by PDE methods (unlike in, for example, Carr and Ghamami (2017) or Kim and Leung (2016)), including high-dimensional applications.

The remainder of the paper is organized as follows. Section 2 provides the stochastic construction of the branching diffusion with jumps and derives a probabilistic representation of classical solutions of PDE (2.1) in terms of it. In Section 3 we conversely establish sufficient conditions for the branching diffusion representation to yield a viscosity solution of (2.1). Section 4 illustrates how the branching diffusion representation allows for the efficient simulation of solutions via the Monte Carlo method in a stylized high-dimensional example. In Section 5 we showcase our methodology in the context of pricing with a systemically important, defaultable counterparty and a high-dimensional underlying (to wit, the Standard & Poor’s 500 (S&P 500) stock index). Section 6 concludes.

2 Branching diffusion representations for nonlocal partial differential equations

Throughout this paper, we fix a time horizon $T>0$ . The goal of this section is to provide a stochastic representation of a classical solution $u\colon[0,T]\times\mathbb{R}^{d}\to\mathbb{R}$ of a nonlocal PDE of the form

	$\displaystyle\partial_{t}u(t,x)+\mathcal{A}[u](t,x)+\int_{\varXi}f(t,x,\xi,% \mathcal{J}[u](t,x,\xi))\gamma(\mathrm{d}\xi)$	$\displaystyle=0,$		(2.1)
	$\displaystyle u(T,x)$	$\displaystyle=g(x),$		(2.2)

where

\mathcal{A}[u](t,x)\triangleq\mu(t,x)^{\mathrm{T}}\mathrm{D}_{x}u(t,x)+\tfrac{% 1}{2}\operatorname{tr}[\sigma(t,x)\sigma(t,x)^{\mathrm{T}}\mathrm{D}_{x}^{2}u(% t,x)]

is the infinitesimal generator of a diffusion process; $\gamma$ is a distribution on an abstract measurable space $(\varXi,\mathfrak{B})$ ; the nonlinearity

f\colon[0,T]\times\mathbb{R}^{d}\times\varXi\times\mathbb{R}^{m}\to\mathbb{R},% \quad f(t,x,\xi,y)\triangleq\sum_{i\in\mathcal{I}}c_{i}(t,x,\xi)y^{i},

is (multivariate) analytic in $y$ with measurable coefficients $c_{i}\colon[0,T]\times\mathbb{R}^{d}\times\varXi\to\mathbb{R}$ for $i\in\mathcal{I}$ , where $m\in\mathbb{N}$ and $\mathcal{I}\subseteq\mathbb{N}_{0}^{m}$ is a nonempty set of multi-indexes;²² 2 We use standard multi-index notation and write $\smash{y^{i}\triangleq\prod_{\ell=1}^{m}y_{\ell}^{i_{\ell}}}$ for $\smash{y\in\mathbb{R}^{m}}$ and $\smash{|i|\triangleq\sum_{\ell=1}^{m}i_{\ell}}$ , $\smash{i\in\mathbb{N}_{0}^{m}}$ . The case of a multivariate polynomial obtains if $\mathcal{I}$ is finite. In particular, note that, for $(t,x,\xi)\in[0,T]\times\mathbb{R}^{d}\times\varXi$ , $f(t,x,\xi,\mathcal{J}[u](t,x,\xi))=\sum_{i\in\mathcal{I}}c_{i}(t,x,\xi)u(t,% \varGamma_{1}(t,x,\xi))^{i_{1}}\cdots u(t,\varGamma_{m}(t,x,\xi))^{i_{m}}.$ and the nonlocal operator $\mathcal{J}$ inside the nonlinearity is given by $\mathcal{J}[u]\colon[0,T]\times\mathbb{R}^{d}\times\varXi\to\mathbb{R}^{m}$ with coordinates

\displaystyle\mathcal{J}_{\ell}[u](t,x,\xi)\triangleq u(t,\varGamma_{\ell}(t,x% ,\xi)),

where the evaluation loci $\varGamma_{\ell}\colon[0,T]\times\mathbb{R}^{d}\times\varXi\to\mathbb{R}^{d}$ , $\ell\in[1:m]$ , are measurable functions. For the branching diffusion with jumps constructed below, the operator $\mathcal{A}$ describes the dynamics of individual particles between branching times and $\mathcal{J}_{\ell}$ specifies the initial positions of particles with jump mark $\ell$ at branching times.

2.1 Preliminaries

For $\smash{n\in\mathbb{N}}$ , we denote by $\smash{\bm{N}_{n}\triangleq\bigcup_{\nu=1}^{n}\mathbb{N}^{\nu}}$ the set of all $\smash{\mathbb{N}}$ -words with length at most $\smash{n}$ and by $\smash{\bm{N}\triangleq\bigcup_{n\in\mathbb{N}}\bm{N}_{n}}$ the set of finite $\smash{\mathbb{N}}$ -words. We work on a probability space $\smash{(\varOmega,\mathfrak{A},\mathbb{P})}$ that supports the following random variables, all of which are taken to be mutually independent.

•

A family $\smash{\{W^{(k)}\}_{k\in\bm{N}}}$ of independent $\mathbb{R}^{d}$ -valued Brownian motions that serve as driving noise for the emerging diffusion processes with infinitesimal generator $\mathcal{A}$ .
•

A family $\smash{\{\Delta^{(k)}\}_{k\in\bm{N}}}$ of independent and identically distributed (iid) $\varXi$ -valued random variables with distribution $\gamma$ .
•

A family $\smash{\{\tau^{(k)}\}_{k\in\bm{N}}}$ of iid $\smash{\mathbb{R}_{+}}$ -valued random variables serving as lifetimes of the particles underlying the branching mechanism. We assume that the distribution of $\smash{\tau^{(k)}}$ , $k\in\bm{N}$ , admits a continuous density $\smash{\rho\colon\mathbb{R}_{+}\to\mathbb{R}_{+}}$ such that

$F(T)\triangleq\int_{T}^{\infty}\rho(s)\mathrm{d}s>0.$ (2.3)
•

A family $\smash{\{I^{(k)}\}_{k\in\bm{N}}}$ of iid $\mathcal{I}$ -valued random variables modeling the number of offspring of each particle as well as marks for their initial positions. We assume that

$p_{i}\triangleq\mathbb{P}[I^{(k)}=i]>0\quad\text{for all }i\in\mathcal{I},~{}k% \in\bm{N}\quad\text{and}\quad\sum_{i\in\mathcal{I}}|i|p_{i}<\infty.$

In the following, we describe the branching mechanism and the spatial dynamics separately to finally obtain the branching diffusion representation of PDE (2.1).

2.2 Branching mechanism

Figure 1: Illustration of the branching mechanism (without jump marks).

We fix an initial time $t\in[0,T]$ . The branching mechanism is defined by recursion on successive generations for each $\omega\in\varOmega$ .

Given a particle of generation $\smash{n\in\mathbb{N}}$ labeled $\smash{k=(k_{1},\dots,k_{n})\in\mathbb{N}^{n}}$ , we denote its parent particle by $\smash{k-\triangleq(k_{1},\ldots,k_{n-1})\in\mathbb{N}^{n-1}}$ . The branching time of particle $\smash{k}$ is given by $\smash{T^{(k)}_{t}\triangleq(T_{t}^{(k-)}+\tau^{(k)})\wedge T}$ ; on the event $\smash{\{T^{(k)}_{t}<T\}}$ the particle $\smash{k}$ is removed at time $\smash{T^{(k)}_{t}}$ and branches into $\smash{|I^{(k)}|}$ descendants, which are labeled by $\smash{(k_{1},\dots,k_{n},k_{n+1})\in\mathbb{N}^{n+1}}$ for $\smash{k_{n+1}\in[1{:}|I^{(k)}|]}$ . For each $k_{n+1}\in[1{:}|I^{(k)}|]$ we define the jump mark of particle $(k,k_{n+1})$ by

J^{(k,k_{n+1})}\triangleq\ell\quad\text{for~{}}\sum_{j=1}^{\ell-1}I^{(k)}_{j}<% k_{n+1}\leq\sum_{j=1}^{\ell}I^{(k)}_{j},

ie, we attach the jump mark $1$ to the first $I^{(k)}_{1}$ descendants of $k$ , the mark $2$ to the following $I^{(k)}_{2}$ descendants, etc. This iteration is well defined and uniquely determines the branching dynamics if we assume that the mechanism starts with a single particle with label $(1)$ of generation $1$ at time $t$ and

(1)-\triangleq(\phantom{0})\triangleq\emptyset\quad\text{and}\quad T^{% \emptyset}_{t}\triangleq t.

In the following, we only refer to $k\in\bm{N}$ as a particle if either $k=(1)$ or $k-$ is a particle and $k_{n}\in[1{:}|I^{(k-)}|]$ . Figure 1 visualizes this branching mechanism.

We denote the random set of all particles of generation $n\in\mathbb{N}$ alive at time $s\in[t,T]$ by

\mathcal{K}^{n}_{t}(s)\triangleq\begin{cases}\{k\in\mathbb{N}^{n}\colon k\text% { is a particle and }T_{t}^{(k-)}\leq s<T^{(k)}_{t}\}&\text{if }s\in[t,T),\\ \\ \{k\in\mathbb{N}^{n}\colon k\text{ is a particle and }T^{(k)}_{t}=T\}&\text{if% }s=T.\end{cases}

The set of all particles of generation $n$ alive before or at time $s$ is given by³³ 3 Since $\{\mathcal{K}^{n}_{t}(s)\}_{s\in[t,T]}$ is piecewise constant, the set $\bigcup_{r\in[t,s]}\mathcal{K}^{n}_{t}(r)=\bigcup_{\begin{subarray}{c}r\in[t,s% ]\cap\mathbb{Q},\\ r=s\text{~{}or~{}}r=t\end{subarray}}\mathcal{K}^{n}_{t}(r)$ is measurable.

\bar{\mathcal{K}}^{n}_{t}(s)\triangleq\bigcup_{r\in[t,s]}\mathcal{K}^{n}_{t}(r).

Finally, the set of all particles alive at time $s$ and the set of all particles alive before or at time $s$ are defined as

\mathcal{K}_{t}(s)\triangleq\bigcup_{n\in\mathbb{N}}\mathcal{K}^{n}_{t}(s)% \quad\text{and}\quad\bar{\mathcal{K}}_{t}(s)\triangleq\bigcup_{n\in\mathbb{N}}% \bar{\mathcal{K}}^{n}_{t}(s),

respectively. For ease of notation, we subsequently write

\mathcal{K}^{n}_{t}\triangleq\mathcal{K}^{n}_{t}(T),\qquad\mathcal{K}_{t}% \triangleq\mathcal{K}_{t}(T),\qquad\bar{\mathcal{K}}^{n}_{t}\triangleq\bar{% \mathcal{K}}^{n}_{t}(T),\qquad\bar{\mathcal{K}}_{t}\triangleq\bar{\mathcal{K}}% _{t}(T).

As in Henry-Labordère et al (2019, Proposition 2.4), the total number of particles is almost surely (a.s.) finite, ie,

\#\bar{\mathcal{K}}_{t}<\infty

(see also Athreya and Ney 1972, Theorem IV.1.1; Harris 1963, Chapter VI, Section 12ff).

2.3 Branching diffusion dynamics

The next step is to specify the dynamics of the individual particles. We first impose some standard regularity assumptions on the coefficient functions in the infinitesimal generator $\mathcal{A}$ .

Assumption 2.1.

The functions

\mu\colon[0,T]\times\mathbb{R}^{d}\to\mathbb{R}^{d}\quad\text{and}\quad\sigma% \colon[0,T]\times\mathbb{R}^{d}\to\mathbb{R}^{d\times d}

are measurable and satisfy the following Lipschitz and linear-growth conditions: there exists a constant $L>0$ such that

	$\displaystyle\|\mu(t,x)-\mu(t,y)\|+\|\sigma(t,x)-\sigma(t,y)\|$	$\displaystyle\leq L\|x-y\|,$	$\displaystyle t$	$\displaystyle\in[0,T],$	$\displaystyle x,y$	$\displaystyle\in\mathbb{R}^{d},$
	$\displaystyle\|\mu(t,x)\|^{2}+\|\sigma(t,x)\|^{2}$	$\displaystyle\leq L^{2}(1+\|x\|^{2}),$	$\displaystyle t$	$\displaystyle\in[0,T],$	$\displaystyle x$	$\displaystyle\in\mathbb{R}^{d}.$

Under this assumption, classical results such as Theorem 3.21 in Pardoux and Răşcanu (2014) imply that the stochastic differential equation

\left.\begin{aligned} \displaystyle\bar{X}^{t,x}_{s}&\displaystyle=x,&&% \displaystyle s\in[0,t],\\ \\ \displaystyle\mathrm{d}\bar{X}^{t,x}_{s}&\displaystyle=\mu(s,\bar{X}^{t,x}_{s}% )\mathrm{d}s+\sigma(s,\bar{X}^{t,x}_{s})\mathrm{d}\bar{W}_{s},&&\displaystyle s% \in[t,T],\end{aligned}\right\}

(2.4)

admits a pathwise unique strong solution for each starting configuration $(t,x)\in[0,T]\times\mathbb{R}^{d}$ . Here, $\smash{\bar{W}}$ is an $\smash{\mathbb{R}^{d}}$ -valued Brownian motion on $\smash{(\varOmega,\mathfrak{A},\mathbb{P})}$ that is independent of all other random variables defined so far in the paper. The natural filtration of $\smash{\bar{X}^{t,x}}$ augmented by all $\mathbb{P}$ -null sets is denoted by $\smash{\bar{\mathfrak{F}}^{t,x}=\{\bar{\mathfrak{F}}^{t,x}_{s}\}_{s\in[0,T]}}$ . An application of the results in Pardoux and Răşcanu (2014, Chapter 3.7) yields the following.

Lemma 2.2 (Properties of the diffusion).

The random field $\smash{\{\bar{X}^{t,x}_{s}\}_{s,t\in[0,T],x\in\mathbb{R}^{d}}}$ defined via (2.4) can be chosen such that it satisfies the following conditions.

(i) Continuity with respect to initial data.

The map

\bar{X}\colon[0,T]\times\mathbb{R}^{d}\times[0,T]\to\mathbb{R}^{d},\quad(t,x,s% )\mapsto\bar{X}^{t,x}_{s},

is a.s. continuous.

(ii) Flow property.

For all $(t,x)\in[0,T]\times\mathbb{R}^{d}$ and any $[t,T]$ -valued $\smash{\bar{\mathfrak{F}}^{t,x}}$ -stopping time $\tau$ ,

\bar{X}^{\tau,\bar{X}^{t,x}_{\tau}}_{\tau+s}\bm{1}_{\{\tau+s\leq T\}}=\bar{X}^% {t,x}_{\tau+s}\bm{1}_{\{\tau+s\leq T\}},\quad s\in[0,\infty).

(iii) Strong Markov property.

For all $\smash{(t,x)\in[0,T]\times\mathbb{R}^{d}}$ and $\smash{s\in[0,\infty)}$ , for any $\smash{[t,T]}$ -valued $\smash{\bar{\mathfrak{F}}^{t,x}}$ -stopping time $\smash{\tau}$ and for every bounded measurable function $\smash{h\colon[0,T]\times\mathbb{R}^{d}\to\mathbb{R}^{d}}$ , we have

\mathbb{E}[h(\tau+s,\bar{X}^{t,x}_{\tau+s})\bm{1}_{\{\tau+s\leq T\}}\mid\bar{% \mathfrak{F}}^{t,x}_{\tau}]=\mathbb{E}[h(\tau+s,\bar{X}^{t,x}_{\tau+s})\bm{1}_% {\{\tau+s\leq T\}}\mid(\tau,\bar{X}^{t,x}_{\tau})].

The branching diffusion is constructed by attaching to each particle in the branching mechanism a diffusion with the same dynamics as $\smash{\bar{X}}$ but with a different driving noise and a suitable initial condition. More precisely, we fix $x\in\mathbb{R}^{d}$ and define for each $\smash{k\in\bar{\mathcal{K}}_{t}}$ with $\smash{k=(1,k_{2},\dots,k_{n})\in\mathbb{N}^{n}}$ an associated diffusion

\smash{X^{(k)}=X^{k,t,x}=\{X^{k,t,x}_{s}\}_{s\in[T_{t}^{(k-)},T^{(k)}_{t}]}}

as the unique strong solution of

	$\displaystyle X^{(k)}_{T_{t}^{(k-)}}$	$\displaystyle=\varGamma_{J^{(k)}}(T_{t}^{(k-)},X^{(k-)}_{T_{t}^{(k-)}},\Delta^% {(k-)}),$
	$\displaystyle\mathrm{d}X^{(k)}_{s}$	$\displaystyle=\mu(s,X_{s}^{(k)})\mathrm{d}s+\sigma(s,X_{s}^{(k)})\mathrm{d}W^{% (k)}_{s},\quad s\in[T^{(k-)}_{t},T^{(k)}_{t}].$

Note that this iteration is well defined starting from $X^{(1)}_{t}=x$ . It follows that $X^{(k)}$ has the same dynamics as $\smash{\bar{X}}$ but with the different, independent driving noise $\smash{W^{(k)}}$ . The lifetime of $\smash{X^{(k)}}$ coincides with the lifetime of the particle $k$ . Moreover, the initial value of $\smash{X^{(k)}}$ is determined by the terminal value of the diffusion $\smash{X^{(k-)}}$ associated with its parent particle $k-$ , the map $\smash{\varGamma_{J^{(k)}}}$ corresponding to its mark $\smash{J^{(k)}}$ and the jump parameter $\smash{\Delta^{(k-)}}$ . For later reference, we encode all information available up to generation $n\in\mathbb{N}$ by setting

\mathfrak{F}^{n}\triangleq\sigma(W^{(k)},\tau^{(k)},\Delta^{(k)},I^{(k)}\colon k% \in\bm{N}_{n}).

For notational convenience, we furthermore write $\mathfrak{F}^{0}\triangleq\{\emptyset,\varOmega\}$ for the trivial $\sigma$ -algebra. Finally, for $n\in\mathbb{N}_{0}$ , we enlarge these $\sigma$ -algebras by the branching time information of the next generation, ie, we set

\mathfrak{G}^{n}\triangleq\mathfrak{F}^{n}\vee\sigma(\tau^{(k)}\colon k\in\bm{% N}_{n+1}).

For $n\in\mathbb{N}$ and $\smash{k\in\bar{\mathcal{K}}^{n}_{t}}$ , we observe that, conditional on $\mathfrak{G}^{n-1}$ , the laws of

X^{(k)}\quad\text{and}\quad\bar{X}^{T_{t}^{(k-)},X^{(k)}_{T_{t}^{(k-)}}}

on $[T_{t}^{(k-)},T^{(k)}_{t}]$ are identical. We stress that

X^{(k)}\quad\text{and}\quad\bar{X}^{T_{t}^{(k-)},X^{(k)}_{T_{t}^{(k-)}}}

do not coincide pathwise since the dynamics of $X^{(k)}$ are driven by $W^{(k)}$ , while those of

\bar{X}^{T_{t}^{(k-)},X^{(k)}_{T_{t}^{(k-)}}}

are driven by $\smash{\bar{W}}$ . Replacing the driving Brownian motion with a new, independent one for each offspring particle – without changing its distribution – will be the key step in the branching representation below.

2.4 Branching diffusion representation

We now address the branching diffusion representation of classical solutions of PDE (2.1). To begin with, we specify suitable boundedness assumptions on the coefficient functions $c_{i}$ and the terminal condition $g$ .

Assumption 2.3.

The functions $g\colon\mathbb{R}^{d}\to\mathbb{R}$ and $c_{i}\colon[0,T]\times\mathbb{R}^{d}\times\varXi\to\mathbb{R}$ , $i\in\mathcal{I}$ , in PDE (2.1) and boundary condition (2.2) are bounded and measurable.

In order for the possibly infinite series in the nonlinearity of PDE (2.1) to be defined unambiguously, we subsequently use the following standard definition of classical solutions.

Definition 2.4 (Classical solution).

Under Assumption 2.3, a continuous function

u\colon[0,T]\times\mathbb{R}^{d}\to\mathbb{R},\quad(t,x)\mapsto u(t,x)

is said to be a classical solution of PDE (2.1) with terminal condition (2.2) if

(i)

$u\in\mathcal{C}^{1,2}([0,T)\times\mathbb{R}^{d})$ ;
(ii)

for each $(t,x)\in[0,T)\times\mathbb{R}^{d}$ , it holds that

$\sum_{i\in\mathcal{I}}\int_{\varXi}|c_{i}(t,x,\xi)\mathcal{J}[u](t,x,\xi)^{i}|% \gamma(\mathrm{d}\xi)<+\infty;$
(iii)

$u$ satisfies PDE (2.1) for each $(t,x)\in[0,T)\times\mathbb{R}^{d}$ and BC (2.2) for each $x\in\mathbb{R}^{d}$ .

We are now in a position to establish the main result of this section, which allows us to represent a classical solution of PDE (2.1) by means of a functional of the branching diffusion. The key idea is to introduce randomization across subsequent generations and subsequently exploit the conditional independence structure.

Remark 2.5 (Randomization).

We fix a particle $k\in\bar{\mathcal{K}}^{n}_{t}$ of generation $n\in\mathbb{N}$ . By a slight abuse of notation, we drop any indexes pertaining to the initial position $(t,x)$ and write

\bar{X}\triangleq\bar{X}^{T^{(k-)},X^{(k)}_{T^{(k-)}}}\quad\text{and}\quad(% \Delta,\tau,I)\triangleq(\Delta^{(k)},\tau^{(k)},I^{(k)}).

For any classical solution $u$ of PDE (2.1), under suitable regularity and integrability assumptions, we have by a Feynman–Kac argument⁴⁴ 4 See (2.10) in the proof of Theorem 2.6 for details.

	$\displaystyle u(T^{(k-)},X^{(k)}_{T^{(k-)}})$
	$\displaystyle\quad=\mathbb{E}\bigg{[}g(\bar{X}_{T})+\int_{T^{(k-)}}^{T}f(s,% \bar{X}_{s},\Delta,\mathcal{J}[u](s,\bar{X}_{s},\Delta))\mathrm{d}s\biggm{\|}% \mathfrak{F}^{n-1}\bigg{]}$
	$\displaystyle\quad=\mathbb{E}\bigg{[}g(\bar{X}_{T})+\int_{T^{(k-)}}^{T}\sum_{i% \in\mathcal{I}}c_{i}(s,\bar{X}_{s},\Delta)\mathcal{J}[u](s,\bar{X}_{s},\Delta)% ^{i}\mathrm{d}s\biggm{\|}\mathfrak{F}^{n-1}\bigg{]}.$		(2.5)

The key idea underlying the branching diffusion representation is to represent the right-hand side recursively in terms of the branching diffusion $X^{(k)}$ , thus eliminating the integral, sum and nonlinearity within the conditional expectation. More precisely, we claim that

	$\displaystyle u(T^{(k-)},X^{(k)}_{T^{(k-)}})$
	$\displaystyle\qquad=\mathbb{E}\bigg{[}\bm{1}_{\{T^{(k)}=T\}}\frac{g(X^{(k)}_{T% })}{F(T-T^{(k-)})}+\bm{1}_{\{T^{(k)}<T\}}\frac{c_{I}(T^{(k)},X^{(k)}_{T^{(k)}}% ,\Delta)}{\rho(T^{(k)}-T^{(k-)})p_{I}}\mathcal{J}[u](T^{(k)},X^{(k)}_{T^{(k)}}% ,\Delta)^{I}\biggm{\|}\mathfrak{F}^{n-1}\bigg{]}.$		(2.6)

Let us start by considering the first summand in (2.6). Since

T^{(k)}=(T^{(k-)}+\tau)\wedge T

and $\smash{\bar{X}}$ and $\smash{X^{(k)}}$ have the same conditional distribution given $\mathfrak{G}^{n-1}$ , the tower property and the fact that $\mathfrak{F}^{n-1}\subseteq\mathfrak{G}^{n-1}$ yield

\mathbb{E}\bigg{[}\bm{1}_{\{T^{(k)}=T\}}\frac{g(X^{(k)}_{T})}{F(T-T^{(k-)})}% \biggm{|}\mathfrak{F}^{n-1}\bigg{]}=\mathbb{E}\bigg{[}\bm{1}_{\{\tau\geq T-T^{% (k-)}\}}\frac{g(\bar{X}_{T})}{F(T-T^{(k-)})}\biggm{|}\mathfrak{F}^{n-1}\bigg{]}.

Then, since $\smash{\tau}$ is independent of $\smash{\mathfrak{F}^{n-1}}$ and $\smash{\bar{X}}$ , we can simply integrate with respect to the density of $\tau$ and use the definition of $F$ in (2.3) to obtain

	$\displaystyle\mathbb{E}\bigg{[}\bm{1}_{\{T^{(k)}=T\}}\frac{g(X^{(k)}_{T})}{F(T% -T^{(k-)})}\biggm{\|}\mathfrak{F}^{n-1}\bigg{]}$	$\displaystyle=\mathbb{E}\bigg{[}\frac{g(\bar{X}_{T})}{F(T-T^{(k-)})}F(T-T^{(k-% )})\biggm{\|}\mathfrak{F}^{n-1}\bigg{]}$
		$\displaystyle=\mathbb{E}[g(\bar{X}_{T})\mid\mathfrak{F}^{n-1}]$

as in (2.5). The second term in (2.6) is slightly more involved, but can be handled similarly. First, we use the conditional identity, given $\smash{\mathfrak{G}^{n-1}}$ , of the laws of $\smash{\bar{X}}$ and $\smash{X^{(k)}}$ to obtain

	$\displaystyle\mathbb{E}\bigg{[}\bm{1}_{\{T^{(k)}<T\}}\frac{c_{I}(T^{(k)},X^{(k% )}_{T^{(k)}},\Delta)}{\rho(T^{(k)}-T^{(k-)})p_{I}}\mathcal{J}[u](T^{(k)},X^{(k% )}_{T^{(k)}},\Delta)^{I}\biggm{\|}\mathfrak{F}^{n-1}\bigg{]}$
	$\displaystyle\qquad\qquad=\mathbb{E}\bigg{[}\bm{1}_{\{T^{(k)}<T\}}\frac{c_{I}(% T^{(k)},\bar{X}_{T^{(k)}},\Delta)}{\rho(T^{(k)}-T^{(k-)})p_{I}}\mathcal{J}[u](% T^{(k)},\bar{X}_{T^{(k)}},\Delta)^{I}\biggm{\|}\mathfrak{F}^{n-1}\bigg{]}.$

Next, independence of $I$ and $\tau$ from all other objects involved allows us to integrate with respect to the associated probability mass and density functions, and we obtain

	$\displaystyle\mathbb{E}\bigg{[}\bm{1}_{\{T^{(k)}<T\}}\frac{c_{I}(T^{(k)},\bar{% X}_{T^{(k)}},\Delta)}{\rho(T^{(k)}-T^{(k-)})p_{I}}\mathcal{J}[u](T^{(k)},\bar{% X}_{T^{(k)}},\Delta)^{I}\biggm{\|}\mathfrak{F}^{n-1}\bigg{]}$
	$\displaystyle\qquad\qquad\qquad\qquad\qquad=\mathbb{E}\bigg{[}\int_{T^{(k-)}}^% {T}\sum_{i\in\mathcal{I}}c_{i}(s,\bar{X}_{s},\Delta)\mathcal{J}[u](s,\bar{X}_{% s},\Delta)^{i}\mathrm{d}s\biggm{\|}\mathfrak{F}^{n-1}\bigg{]}.$

Combining the above equations, we have therefore argued that

	$\displaystyle\mathbb{E}\bigg{[}\bm{1}_{\{T^{(k)}<T\}}\frac{c_{I}(T^{(k)},X^{(k% )}_{T^{(k)}},\Delta)}{\rho(T^{(k)}-T^{(k-)})p_{I}}\mathcal{J}[u](T^{(k)},X^{(k% )}_{T^{(k)}},\Delta)^{I}\biggm{\|}\mathfrak{F}^{n-1}\bigg{]}$
	$\displaystyle\qquad\qquad\qquad\qquad\qquad=\mathbb{E}\bigg{[}\int_{T^{(k-)}}^% {T}\sum_{i\in\mathcal{I}}c_{i}(s,\bar{X}_{s},\Delta)\mathcal{J}[u](s,\bar{X}_{% s},\Delta)^{i}\mathrm{d}s\biggm{\|}\mathfrak{F}^{n-1}\bigg{]}$

and thus (2.6) holds. The key advantage of (2.6) is that the nonlinearity is represented via

\mathcal{J}[u](T^{(k)},X^{(k)}_{T^{(k)}},\Delta)^{I}=\prod_{\ell=1}^{m}u(T^{(k% )},\varGamma_{\ell}(T^{(k)},X^{(k)}_{T^{(k)}},\Delta))^{I_{\ell}},

where it is possible to iterate over all generations of particles, using conditional independence across generations: indeed, the terms $\smash{\varGamma_{\ell}(T^{(k)},X^{(k)}_{T^{(k)}},\Delta^{(k)})}$ appearing as arguments in the function $u$ correspond to the initial positions of the particles of generation $n+1$ . This leads to the branching diffusion representation made rigorous in Theorem 2.6.

We next state the first main result of this paper: a branching diffusion representation of classical solutions of PDE (2.1).

Theorem 2.6.

(Branching representation of classical solutions) Let Assumptions 2.1 and 2.3 hold, $u\colon[0,T]\times\mathbb{R}^{d}\to\mathbb{R}$ be a classical solution of PDE (2.1) satisfying BC (2.2) and $(t,x)\in[0,T]\times\mathbb{R}^{d}$ . Suppose that

(i)

the family $\smash{\{\varPsi^{t,x}_{n}\}_{n\in\mathbb{N}_{0}}}$ defined in (2.8) is uniformly integrable and
(ii)

the functions $u$ and

$[0,T]\times\mathbb{R}^{d}\to\mathbb{R},\quad(t,x)\mapsto\sum_{i\in\mathcal{I}}% \int_{\varXi}|c_{i}(t,x,\xi)\mathcal{J}[u](t,x,\xi)^{i}|\gamma(\mathrm{d}\xi)$

satisfy a polynomial growth condition.⁵⁵ 5 $u\colon[0,T]\times\mathbb{R}^{d}\to\mathbb{R}$ is said to satisfy a polynomial growth condition if there are constants $K,\kappa>0$ such that $\smash{\sup_{t\in[0,T]}|u(t,x)|\leq K(1+|x|^{\kappa})}$ .

Then, $u$ admits the branching diffusion representation

\displaystyle u(t,x)=\mathbb{E}[\varPsi^{t,x}],\qquad\text{where}\varPsi^{t,x}% \triangleq\prod_{k\in\mathcal{K}}\frac{g(X^{(k)}_{T})}{F(T-T^{(k-)})}\prod_{k% \in\bar{\mathcal{K}}\setminus\mathcal{K}}\frac{c_{I^{(k)}}(T^{(k)},X^{(k)}_{T^% {(k)}},\Delta^{(k)})}{\rho(T^{(k)}-T^{(k-)})p_{I^{(k)}}}.

(2.7)

The representation (2.7) provides the basis for a direct (nonnested) Monte Carlo method by simulating a diffusion on $[0,T]$ that branches at an a.s. finite number of times with a.s. finitely many offspring; see Sections 4 and 5.

Proof of Theorem 2.6.

For $n\in\mathbb{N}_{0}$ , iteratively define the random variables,⁶⁶ 6 For notational convenience, we suppress the dependence on $\smash{(t,x)}$ on the right-hand sides of these definitions but highlight here that $\smash{X^{(k)}=X^{k,t,x}}$ and $\smash{T^{(k)}=T^{(k)}_{t}}$ and also that $\smash{\bar{\mathcal{K}}^{n}=\bar{\mathcal{K}}^{n}_{t}}$ and $\mathcal{K}^{n}=\mathcal{K}^{n}_{t}$ .

	$\displaystyle\mathcal{G}^{t,x}_{0}$	$\displaystyle\triangleq\mathcal{C}^{t,x}_{0}\triangleq 1,$
	$\displaystyle\mathcal{G}^{t,x}_{n}$	$\displaystyle\triangleq\mathcal{G}^{t,x}_{n-1}\prod_{k\in\mathcal{K}^{n}}\frac% {g(X^{(k)}_{T})}{F(T-T^{(k-)})},$
	$\displaystyle\mathcal{C}^{t,x}_{n}$	$\displaystyle\triangleq\mathcal{C}^{t,x}_{n-1}\prod_{k\in\bar{\mathcal{K}}^{n}% \setminus\mathcal{K}^{n}}\frac{c_{I^{(k)}}(T^{(k)},X^{(k)}_{T^{(k)}},\Delta^{(% k)})}{\rho(T^{(k)}-T^{(k-)})p_{I^{(k)}}}$

and

\displaystyle\mathcal{R}^{t,x}_{n}\triangleq\prod_{k\in\mathcal{K}^{n+1}}\frac% {g(X^{(k)}_{T})}{F(T-T^{(k-)})}\prod_{k\in\bar{\mathcal{K}}^{n+1}\setminus% \mathcal{K}^{n+1}}\frac{c_{I^{(k)}}(T^{(k)},X^{(k)}_{T^{(k)}},\Delta^{(k)})% \mathcal{J}[u](T^{(k)},X^{(k)}_{T^{(k)}},\Delta^{(k)})^{I^{(k)}}}{\rho(T^{(k)}% -T^{(k-)})p_{I^{(k)}}},

and set

\varPsi^{t,x}_{n}\triangleq\mathcal{G}^{t,x}_{n}\mathcal{C}^{t,x}_{n}\mathcal{% R}^{t,x}_{n}.

(2.8)

Since $\lim_{n\to\infty}\varPsi^{t,x}_{n}=\varPsi^{t,x}$ and $\{\varPsi^{t,x}_{n}\}_{n\in\mathbb{N}_{0}}$ is uniformly integrable, it follows from Vitali’s theorem that $\varPsi^{t,x}$ is integrable and $\lim_{n\to\infty}\mathbb{E}[\varPsi^{t,x}_{n}]=\mathbb{E}[\varPsi^{t,x}]$ . Hence, it suffices to show that

u(t,x)=\mathbb{E}[\varPsi^{t,x}_{n}]\quad\text{for each }n\in\mathbb{N}_{0}.

(2.9)

Step 1. Fix a particle $\smash{k\in\bar{\mathcal{K}}^{n}}$ of generation $\smash{n\in\mathbb{N}}$ and write

\bar{X}^{(k)}\triangleq\bar{X}^{T^{(k-)},X^{(k)}_{T^{(k-)}}}.

Setting $\smash{\tau_{\nu}\triangleq\inf\{r\geq T^{(k-)}\colon|\bar{X}^{(k)}_{r}|\geq% \nu\}\wedge T}$ for each $\nu\in\mathbb{N}$ , Itô’s lemma and the fact that $u$ solves PDE (2.1) yield

	$\displaystyle u(T^{(k-)},\bar{X}^{(k)}_{T^{(k-)}})$
	$\displaystyle\qquad=u(\tau_{\nu},\bar{X}^{(k)}_{\tau_{\nu}})-\int_{T^{(k-)}}^{% \tau_{\nu}}[\partial_{t}u(s,\bar{X}^{(k)}_{s})+\mathcal{A}[u](s,\bar{X}^{(k)}_% {s})]\mathrm{d}s-M_{\tau_{\nu}}$
	$\displaystyle\qquad=u(\tau_{\nu},\bar{X}^{(k)}_{\tau_{\nu}})+\int_{T^{(k-)}}^{% \tau_{\nu}}\int_{\varXi}f(s,\bar{X}^{(k)}_{s},\xi,\mathcal{J}[u](s,\bar{X}^{(k% )}_{s},\xi))\gamma(\mathrm{d}\xi)\mathrm{d}s-M_{\tau_{\nu}},$

where the local martingale part

\smash{M\triangleq\int_{T^{(k-)}}^{\cdot\wedge\tau_{\nu}}\mathrm{D}_{x}u(s,% \bar{X}^{(k)}_{s})\sigma(s,\bar{X}^{(k)}_{s})\mathrm{d}\bar{W}_{s}}

is in fact a martingale. Since $\smash{(T^{(k-)},X^{(k)}_{T^{(k-)}})}$ is $\smash{\mathfrak{F}^{n-1}}$ -measurable, taking conditional expectations yields

	$\displaystyle u(T^{(k-)},\bar{X}^{(k)}_{T^{(k-)}})$
	$\displaystyle\quad=\mathbb{E}\bigg{[}u(\tau_{\nu},\bar{X}^{(k)}_{\tau_{\nu}})+% \int_{T^{(k-)}}^{\tau_{\nu}}\int_{\varXi}f(s,\bar{X}^{(k)}_{s},\xi,\mathcal{J}% [u](s,\bar{X}^{(k)}_{s},\xi))\gamma(\mathrm{d}\xi)\mathrm{d}s\biggm{\|}% \mathfrak{F}^{n-1}\bigg{]}.$

It is clear that, as $\nu\to\infty$ ,

	$\displaystyle u(\tau_{\nu},\bar{X}^{(k)}_{\tau_{\nu}})$	$\displaystyle\to u(T,\bar{X}^{(k)}_{T})=g(\bar{X}^{(k)}_{T}),$
	$\displaystyle\int_{T^{(k)}}^{\tau_{\nu}}\int_{\varXi}f(s,\bar{X}^{(k)}_{s},\xi% ,\mathcal{J}[u]$	$\displaystyle(s,\bar{X}^{(k)}_{s},\xi))\gamma(\mathrm{d}\xi)\mathrm{d}s$
		$\displaystyle\to\int_{T^{(k)}}^{T}\int_{\varXi}f(s,\bar{X}^{(k)}_{s},\xi,% \mathcal{J}[u](s,\bar{X}^{(k)}_{s},\xi))\gamma(\mathrm{d}\xi)\mathrm{d}s$

and using (ii) it follows by a dominated convergence argument as in Karatzas and Shreve (1998, Theorem 5.7.6) that

	$\displaystyle u(T^{(k-)},X^{(k)}_{T^{(k-)}})$
	$\displaystyle\qquad=\mathbb{E}\bigg{[}g(\bar{X}^{(k)}_{T})+\int_{T^{(k-)}}^{T}% \int_{\varXi}f(s,\bar{X}^{(k)}_{s},\xi,\mathcal{J}[u](s,\bar{X}^{(k)}_{s},\xi)% )\gamma(\mathrm{d}\xi)\mathrm{d}s\biggm{\|}\mathfrak{F}^{n-1}\bigg{]},$

where we have used the fact that $\smash{\bar{X}^{(k)}_{T^{(k-)}}=X^{(k)}_{T^{(k-)}}}$ by definition of $\smash{\bar{X}^{(k)}}$ . Since $\smash{\Delta^{(k)}}$ has distribution $\smash{\gamma}$ and is independent of $\smash{\mathfrak{F}^{n-1}}$ , and $\smash{\bar{X}^{(k)}}$ and $\smash{\Delta^{(k)}}$ are independent, it follows that

	$\displaystyle u(T^{(k-)},X^{(k)}_{T^{(k-)}})$
	$\displaystyle\qquad=\mathbb{E}\bigg{[}g(\bar{X}^{(k)}_{T})+\int_{T^{(k-)}}^{T}% f(s,\bar{X}^{(k)}_{s},\Delta^{(k)},\mathcal{J}[u](s,\bar{X}^{(k)}_{s},\Delta^{% (k)}))\mathrm{d}s\biggm{\|}\mathfrak{F}^{n-1}\bigg{]},$		(2.10)

and by the argument in Remark 2.5 this can be further rewritten as

	$\displaystyle u(T^{(k-)},X^{(k)}_{T^{(k-)}})$
	$\displaystyle=\mathbb{E}\bigg{[}\bm{1}_{\{T^{(k)}=T\}}\frac{g(X^{(k)}_{T})}{F(% T-T^{(k-)})}+\bm{1}_{\{T^{(k)}<T\}}\frac{c_{I^{(k)}}(T^{(k)},X^{(k)}_{T^{(k)}}% ,\Delta^{(k)})}{\rho(T^{(k)}-T^{(k-)})p_{I^{(k)}}}\mathcal{J}[u](T^{(k)},X^{(k% )}_{T^{(k)}},\Delta^{(k)})^{I^{(k)}}\biggm{\|}\mathfrak{F}^{n-1}\bigg{]}.$		(2.11)

Step 2. We establish (2.9) by induction on $n$ . For $n=0$ , let $k=(1)$ be the only particle of generation $1$ , recall that $\mathfrak{F}^{0}$ is trivial and note that (2.11) can be rewritten as

	$\displaystyle u(t,x)$	$\displaystyle=\mathbb{E}\bigg{[}\bm{1}_{\{T^{(k)}=T\}}\frac{g(X^{(k)}_{T})}{F(% T-T^{(k-)})}+\bm{1}_{\{T^{(k)}<T\}}\frac{c_{I^{(k)}}(T^{(k)},X^{(k)}_{T^{(k)}}% ,\Delta^{(k)})\mathcal{J}[u](T^{(k)},X^{(k)}_{T^{(k)}},\Delta^{(k)})^{I^{(k)}}% }{\rho(T^{(k)}-T^{(k-)})p_{I^{(k)}}}\bigg{]}$
		$\displaystyle=\mathbb{E}[\mathcal{R}^{t,x}_{0}]=\mathbb{E}[\varPsi_{0}^{t,x}].$

Now, let $n\in\mathbb{N}$ and suppose that the claim is true for $n-1$ , ie,

u(t,x)=\mathbb{E}[\varPsi_{n-1}^{t,x}]=\mathbb{E}[\mathcal{G}^{t,x}_{n-1}% \mathcal{C}^{t,x}_{n-1}\mathcal{R}^{t,x}_{n-1}].

(2.12)

Let $\smash{k\in\bar{\mathcal{K}}^{n}}$ be an arbitrary particle of generation $n$ . On the event

\{k\in\bar{\mathcal{K}}^{n}\setminus\mathcal{K}^{n}\}=\{T^{(k)}<T\},

the particle $\smash{k}$ branches into $\smash{|I^{(k)}|}$ offspring particles $\smash{(k,k_{n+1})}$ , $\smash{k_{n+1}\in[1{:}|I^{(k)}|]}$ , of which the first $\smash{I^{(k)}_{1}}$ have mark $\smash{1}$ , ie, jump to $\smash{\varGamma_{1}(T^{(k)},X^{(k)}_{T^{(k)}},\Delta^{(k)})}$ , the next $\smash{I^{(k)}_{2}}$ have mark $\smash{2}$ , ie, jump to $\smash{\varGamma_{2}(T^{(k)},X^{(k)}_{T^{(k)}},\Delta^{(k)})}$ , and so forth. Thus, on the event $\smash{\{k\in\bar{\mathcal{K}}^{n}\setminus\mathcal{K}^{n}\}}$ , we have

	$\displaystyle\mathcal{J}[u](T^{(k)},X^{(k)}_{T^{(k)}},\Delta^{(k)})^{I^{(k)}}$
	$\displaystyle\quad=\prod_{\ell=1}^{m}u(T^{(k)},\varGamma_{\ell}(T^{(k)},X^{(k)% }_{T^{(k)}},\Delta^{(k)}))^{I^{(k)}_{\ell}}$
	$\displaystyle\quad=\prod_{k_{n+1}=1}^{\|I^{(k)}\|}u(T^{(k)},X_{T^{(k)}}^{(k,k_{n% +1})})=\prod_{\bar{k}\in\bar{\mathcal{K}}^{n+1},\bar{k}-=k}u(T^{(\bar{k}-)},X_% {T^{(\bar{k}-)}}^{(\bar{k})})$
	$\displaystyle\quad=\prod_{\bar{k}\in\bar{\mathcal{K}}^{n+1},\bar{k}-=k}\mathbb% {E}\bigg{[}\bm{1}_{\{T^{(\bar{k})}=T\}}\frac{g(X^{(\bar{k})}_{T})}{F(T-T^{(% \bar{k}-)})}+\bm{1}_{\{T^{(\bar{k})}<T\}}\frac{c_{I^{(\bar{k})}}(T^{(\bar{k})}% ,X^{(\bar{k})}_{T^{(\bar{k})}},\Delta^{(\bar{k})})\mathcal{J}[u](T^{(\bar{k})}% ,X^{(\bar{k})}_{T^{(\bar{k})}},\Delta^{(\bar{k})})^{I^{(\bar{k})}}}{\rho(T^{(% \bar{k})}-T^{(\bar{k}-)})p_{I^{(\bar{k})}}}\biggm{\|}\mathfrak{F}^{n}\bigg{]},$

where the final identity is due to (2.11). Thus, using the $\mathfrak{F}^{n}$ -conditional independence of the individual offspring particles, we have

	$\displaystyle\prod_{k\in\bar{\mathcal{K}}^{n}\setminus\mathcal{K}^{n}}\mathcal% {J}[u](T^{(k)},X^{(k)}_{T^{(k)}},\Delta^{(k)})^{I^{(k)}}$
	$\displaystyle=\mathbb{E}\bigg{[}\prod_{\bar{k}\in\bar{\mathcal{K}}^{n+1}}\bigg% {(}\bm{1}_{\{T^{(\bar{k})}=T\}}\frac{g(X^{(\bar{k})}_{T})}{F(T-T^{(\bar{k}-)})}$
	$\displaystyle\qquad\quad+\bm{1}_{\{T^{(\bar{k})}<T\}}\frac{{c_{I^{(\bar{k})}}(% T^{(\bar{k})},X^{(\bar{k})}_{T^{(\bar{k})}},\Delta^{(\bar{k})})}{\times% \mathcal{J}[u](T^{(\bar{k})},X^{(\bar{k})}_{T^{(\bar{k})}},\Delta^{(\bar{k})})% ^{I^{(\bar{k})}}}}{\rho(T^{(\bar{k})}-T^{(\bar{k}-)})p_{I^{(\bar{k})}}}\bigg{)% }\biggm{\|}\mathfrak{F}^{n}\bigg{]}$
	$\displaystyle=\mathbb{E}[\mathcal{R}_{n}^{t,x}\mid\mathfrak{F}^{n}].$		(2.13)

Then, using (2.13) and the definition of $\smash{\mathcal{R}_{n-1}^{t,x}}$ , we obtain

	$\displaystyle\mathcal{R}^{t,x}_{n-1}$	$\displaystyle=\frac{\mathcal{G}_{n}^{t,x}}{\mathcal{G}_{n-1}^{t,x}}\frac{% \mathcal{C}_{n}^{t,x}}{\mathcal{C}_{n-1}^{t,x}}\prod_{k\in\bar{\mathcal{K}}^{n% }\setminus\mathcal{K}^{n}}\mathcal{J}[u](T^{(k)},X^{(k)}_{T^{(k)}},\Delta^{(k)% })^{I^{(k)}}$
		$\displaystyle=\frac{\mathcal{G}_{n}^{t,x}}{\mathcal{G}_{n-1}^{t,x}}\frac{% \mathcal{C}_{n}^{t,x}}{\mathcal{C}_{n-1}^{t,x}}\mathbb{E}[\mathcal{R}^{t,x}_{n% }\mid\mathfrak{F}^{n}].$

Plugging this into (2.12) yields the claim since $\smash{(\mathcal{G}^{t,x}_{n},\mathcal{C}^{t,x}_{n})}$ is $\smash{\mathfrak{F}^{n}}$ -measurable and thus

u(t,x)=\mathbb{E}[\mathcal{G}^{t,x}_{n-1}\mathcal{C}^{t,x}_{n-1}\mathcal{R}^{t% ,x}_{n-1}]=\mathbb{E}[\mathcal{G}^{t,x}_{n}\mathcal{C}^{t,x}_{n}\mathbb{E}[% \mathcal{R}^{t,x}_{n}\mid\mathfrak{F}^{n}]]=\mathbb{E}[\varPsi_{n}^{t,x}].

∎

Remark 2.7.

A sufficient condition for (ii) in Theorem 2.6 to hold is that $u$ is bounded and $\smash{\sum_{i\in\mathcal{I}}\|c_{i}\|_{\infty}\|u\|_{\infty}^{|i|}<+\infty}$ ; the latter condition holds in particular if $\mathcal{I}$ is finite. If $u$ is bounded, similar arguments to those in Section 3.2 can be used to obtain sufficient conditions for (i). Note that in several applications, including valuation with counterparty risk as in Section 5, the condition that $u$ is bounded can be verified a priori.

The branching diffusion representation in Theorem 2.6 can be extended in several ways: combining our approach with that in Henry-Labordère et al (2019) allows the treatment of mixed local and nonlocal nonlinearities that include first-order derivatives. Moreover, the notion of branching diffusions with jumps can also be extended to that of branching jump-diffusions with jumps, allowing for an additional (linear) nonlocal term in the infinitesimal generator.

3 Viscosity solutions and branching diffusion representations

In the previous section, we started with a classical solution of PDE (2.1) and aimed to derive its branching diffusion representation. This result was achieved under the assumption that a classical solution exists. In this section, we study the converse question: can the branching diffusion representation be used to define a solution of the PDE? It is clear that in general the branching diffusion representation does not yield a sufficiently regular solution to qualify as a classical solution; hence, we subsequently work with the weaker concept of viscosity solutions.

The main result of this section gives sufficient conditions under which the branching diffusion representation defines a viscosity solution of PDE (2.1). We then derive conditions under which the family $\smash{\{\varPsi^{t,x}\}_{(t,x)\in[0,T]\times\mathbb{R}^{d}}}$ is uniformly integrable, which implies one of the key assumptions needed to obtain the viscosity property of the branching diffusion representation.

3.1 Viscosity solutions of nonlocal nonlinear PDEs

We first provide a definition of viscosity solutions of PDE (2.1) that is appropriate for dealing with the nonlocal terms in the nonlinearity.

Definition 3.1.

(Viscosity solution) Suppose Assumption 2.3 holds and let $u\colon[0,T]\times\mathbb{R}^{d}\to\mathbb{R}$ be a continuous function such that

\sum_{i\in\mathcal{I}}\int_{\varXi}|c_{i}(t,x,\xi)\mathcal{J}[u](t,x,\xi)^{i}|% \gamma(\mathrm{d}\xi)<+\infty\quad\text{for all }(t,x)\in[0,T)\times\mathbb{R}% ^{d}.

We say that

(1)

$u$ is a viscosity subsolution of PDE (2.1) if for all $(t,x)\in[0,T)\times\mathbb{R}^{d}$ and all test functions $\varphi\in\mathcal{C}^{1,2}([0,T]\times\mathbb{R}^{d})$ with $\varphi(t,x)=u(t,x)$ and $\varphi\geq u$ we have

$-\partial_{t}\varphi(t,x)-\mathcal{A}[\varphi](t,x)-\int_{\varXi}f(t,x,\xi,% \mathcal{J}[u](t,x,\xi))\gamma(\mathrm{d}\xi)\leq 0;$
(2)

$u$ is a viscosity supersolution of PDE (2.1) if for all $(t,x)\in[0,T)\times\mathbb{R}^{d}$ and all test functions $\varphi\in\mathcal{C}^{1,2}([0,T]\times\mathbb{R}^{d})$ with $\varphi(t,x)=u(t,x)$ and $\varphi\leq u$ we have

$-\partial_{t}\varphi(t,x)-\mathcal{A}[\varphi](t,x)-\int_{\varXi}f(t,x,\xi,% \mathcal{J}[u](t,x,\xi))\gamma(\mathrm{d}\xi)\geq 0;$
(3)

$u$ is a viscosity solution of PDE (2.1) if it is both a viscosity subsolution and a viscosity supersolution.

With this definition in place, we can state our second main result.

Theorem 3.2 (Viscosity property of the branching representation).

For all $(t,x)\in[0,T]\times\smash{\mathbb{R}^{d}}$ , let $\smash{\varPsi^{t,x}}$ be as given in (2.7) and define

u\colon[0,T]\times\mathbb{R}^{d}\to\mathbb{R},\quad(t,x)\mapsto u(t,x)% \triangleq\mathbb{E}[\varPsi^{t,x}].

(3.1)

In addition to Assumptions 2.1 and 2.3, assume that the stochastic differential equation coefficients $\mu$ , $\sigma$ and the PDE coefficients $c_{i}$ , $i\in\mathcal{I}$ , and $g$ are continuous. Moreover, suppose that the following conditions are satisfied for all $(t,x)\in[0,T]\times\mathbb{R}^{d}$ .

(i)

There exists $\varepsilon>0$ such that the family $\smash{\{\varPsi^{s,y}\}_{(s,y)\in\mathcal{B}_{\varepsilon}(t,x)}}$ is uniformly integrable.
(ii)

There exist $\delta>0$ and a measurable function $\zeta\colon\mathcal{I}\times\mathbb{R}^{d}\to\mathbb{R}$ such that

$|c_{i}(s,y,\xi)\mathcal{J}[u](s,y,\xi)^{i}|\leq\zeta(i,\xi)\quad\text{for all % }(s,y)\in\bar{\mathcal{B}}_{\delta}(t,x)$

with

$\sum_{i\in\mathcal{I}}\int_{\varXi}|\zeta(i,\xi)|\gamma(\mathrm{d}\xi)<+\infty.$
(iii)

It holds that

$\sum_{i\in\mathcal{I}}\mathbb{E}\bigg{[}\int_{t}^{T}\int_{\varXi}|c_{i}(s,\bar% {X}^{t,x}_{s},\xi)\mathcal{J}[u](s,\bar{X}^{t,x}_{s},\xi)^{i}|\gamma(\mathrm{d% }\xi)\mathrm{d}s\bigg{]}<+\infty.$

Then, $u$ is a viscosity solution of PDE (2.1).

Note that if $u$ is bounded and $\sum_{i\in\mathcal{I}}\|c_{i}\|_{\infty}\|u\|_{\infty}^{|i|}<+\infty$ (as in Remark 2.7), then conditions (ii) and (iii) are satisfied; sufficient conditions for (i), which simultaneously imply boundedness of $u$ , are given in Section 3.2 below.

Proof.

Note that the uniform integrability assumption (i) implies that $\varPsi^{t,x}\in\mathcal{L}^{1}(\mathbb{P})$ , and hence $u$ is well defined. Moreover, Lemma 2.2 and the continuity of $c_{i}$ , $i\in\mathcal{I}$ , and $g$ guarantee that $(t,x)\mapsto\varPsi^{t,x}$ is continuous, so $u$ is continuous by Vitali’s theorem. Finally, condition (ii) guarantees that $u$ satisfies the integrability conditions required on viscosity solutions in Definition 3.1.

Step 1. We fix $(t,x)\in[0,T)\times\mathbb{R}^{d}$ . As before, we drop the index $(t,x)$ in some of the random variables and processes to simplify notation. Moreover, we set

(X,I,\Delta)\triangleq(X^{(1)},I^{(1)},\Delta^{(1)}).

We first observe that

\bm{1}_{\{T^{(1)}=T\}}\varPsi^{t,x}=\bm{1}_{\{T^{(1)}=T\}}\frac{g(X_{T})}{F(T-% t)}

and also that

	$\displaystyle\bm{1}_{\{T^{(1)}<T\}}\varPsi^{t,x}$
	$\displaystyle\qquad=\bm{1}_{\{T^{(1)}<T\}}\frac{c_{I}(T^{(1)},X_{T^{(1)}},% \Delta)}{\rho(T^{(1)}-t)p_{I}}\prod_{k\in\mathcal{K}\setminus\{(1)\}}\frac{g(X% _{T}^{(k)})}{F(T-T^{(k-)})}\prod_{k\in\bar{\mathcal{K}}\setminus(\mathcal{K}% \cup\{(1)\})}\frac{c_{I^{(k)}}(T^{(k)},X^{(k)}_{T^{(k)}},\Delta^{(k)})}{\rho(T% ^{(k)}-T^{(k-)})p_{I^{(k)}}}$
	$\displaystyle\qquad=\bm{1}_{\{T^{(1)}<T\}}\frac{c_{I}(T^{(1)},X_{T^{(1)}},% \Delta)}{\rho(T^{(1)}-t)p_{I}}\prod_{k_{2}=1}^{\|I\|}\varPsi^{T^{(1)},X^{(1,k_{2% })}_{T^{(1)}}}.$

Using the definition of the branching diffusion and conditional independence structure, we find that

	$\displaystyle\mathbb{E}\bigg{[}\prod_{k_{2}=1}^{\|I\|}\varPsi^{T^{(1)},X^{(1,k_{% 2})}_{T^{(1)}}}\biggm{\|}\mathfrak{F}^{1}\bigg{]}$
	$\displaystyle\qquad\qquad=\prod_{k_{2}=1}^{\|I\|}\sum_{\ell=1}^{m}\bm{1}_{\{J^{(% 1,k_{2})}=\ell\}}\mathbb{E}[\varPsi^{s,y}]\mid_{(s,y)=(T^{(1)},\varGamma_{\ell% }(T^{(1)},X_{T^{(1)}},\Delta))}$
	$\displaystyle\qquad\qquad=\mathcal{J}[u](T^{(1)},X_{T^{(1)}},\Delta)^{I}.$

Putting these equations together and using the tower property of conditional expectation, it follows as in Remark 2.5 that

$\displaystyle u(t,x)$	$\displaystyle=\mathbb{E}[\mathbb{E}[\varPsi^{t,x}\mid\mathfrak{F}^{1}]]$
	$\displaystyle=\mathbb{E}\bigg{[}\bm{1}_{\{T^{(1)}=T\}}\frac{g(X_{T})}{F(T-t)}+% \bm{1}_{\{T^{(1)}<T\}}\frac{c_{I}(T^{(1)},X_{T^{(1)}},\Delta)\mathcal{J}[u](T^% {(1)},X_{T^{(1)}},\Delta)^{I}}{\rho(T^{(1)}-t)p_{I}}\bigg{]}$
	$\displaystyle=\mathbb{E}\bigg{[}g(\bar{X}_{T})+\int_{t}^{T}\int_{\varXi}f(s,% \bar{X}_{s},\xi,\mathcal{J}[u](s,\bar{X}_{s},\xi))\gamma(\mathrm{d}\xi)\mathrm% {d}s\bigg{]}.$	(3.2)

For $\varepsilon>0$ , we now introduce the stopping time

\tau_{\varepsilon}\triangleq\inf\{s\geq t\colon|\bar{X}_{s}-x|\geq\varepsilon% \}\wedge(t+\varepsilon)\wedge T.

From (3.2), the flow property and the strong Markov property of $\smash{\bar{X}}$ as noted in Lemma 2.2, in combination with the conditional version of Fubini’s theorem, which is applicable by (iii), it follows that $u$ satisfies the dynamic programming representation

u(t,x)=\mathbb{E}\bigg{[}u(\tau_{\varepsilon},\bar{X}_{\tau_{\varepsilon}})+% \int_{t}^{\tau_{\varepsilon}}\int_{\varXi}f(s,\bar{X}_{s},\xi,\mathcal{J}[u](s% ,\bar{X}_{s},\xi))\gamma(\mathrm{d}\xi)\mathrm{d}s\bigg{]}.

(3.3)

Step 2. From the dynamic programming representation (3.3), the viscosity property of $u$ follows by standard arguments. To keep the presentation self-contained, we provide a proof of the subsolution property (the supersolution property is established analogously). We fix a test function $\smash{\varphi\in\mathcal{C}^{1,2}([0,T]\times\mathbb{R}^{d})}$ with $\varphi(t,x)=u(t,x)$ and $\varphi\geq u$ . By the dynamic programming representation and Itô’s lemma, we find that

	$\displaystyle\varphi(t,x)$	$\displaystyle=u(t,x)$
		$\displaystyle=\mathbb{E}\bigg{[}u(\tau_{\varepsilon},\bar{X}_{\tau_{% \varepsilon}})+\int_{t}^{\tau_{\varepsilon}}\int_{\varXi}f(s,\bar{X}_{s},\xi,% \mathcal{J}[u](s,\bar{X}_{s},\xi))\gamma(\mathrm{d}\xi)\mathrm{d}s\bigg{]}$
		$\displaystyle\leq\mathbb{E}\bigg{[}\varphi(\tau_{\varepsilon},\bar{X}_{\tau_{% \varepsilon}})+\int_{t}^{\tau_{\varepsilon}}\int_{\varXi}f(s,\bar{X}_{s},\xi,% \mathcal{J}[u](s,\bar{X}_{s},\xi))\gamma(\mathrm{d}\xi)\mathrm{d}s\bigg{]}$
		$\displaystyle=\mathbb{E}\bigg{[}\varphi(t,x)+\int_{t}^{\tau_{\epsilon}}% \partial_{t}\varphi(s,\bar{X}_{s})+\mathcal{A}[\varphi](s,\bar{X}_{s})$
		$\displaystyle\qquad\qquad\qquad\qquad~{}~{}+\int_{\varXi}f(s,\bar{X}_{s},\xi,% \mathcal{J}[u](s,\bar{X}_{s},\xi))\gamma(\mathrm{d}\xi)\mathrm{d}s\bigg{]}.$

For $(s,y)\in[0,T]\times\mathbb{R}^{d}$ , we now define

I_{\varphi}(s,y)\triangleq\partial_{t}\varphi(s,y)+\mathcal{A}[\varphi](s,y)+% \int_{\varXi}f(s,y,\xi,\mathcal{J}[u](s,y,\xi))\gamma(\mathrm{d}\xi)

to arrive at

\mathbb{E}\bigg{[}\int_{t}^{\tau_{\epsilon}}I_{\varphi}(s,\bar{X}_{s})\mathrm{% d}s\bigg{]}\geq 0.

From condition (ii) and the dominated convergence theorem, it follows that $I$ is continuous. Therefore,

0\leq\mathbb{E}\bigg{[}\int_{t}^{\tau_{\epsilon}}I_{\varphi}(s,\bar{X}_{s})% \mathrm{d}s\bigg{]}\leq\mathbb{E}[\tau_{\varepsilon}-t]\max_{(s,y)\in\bar{% \mathcal{B}}_{\varepsilon}(t,x)}I_{\varphi}(s,y)

and thus, since $\mathbb{E}[\tau_{\varepsilon}-t]>0$ ,

\max_{(s,y)\in\bar{\mathcal{B}}_{\varepsilon}(t,x)}I_{\varphi}(s,y)\geq 0.

Letting $\varepsilon\downarrow 0$ allows us to conclude that $I_{\varphi}(t,x)\geq 0$ , so $u$ is a viscosity subsolution in the sense of Definition 3.1. ∎

3.2 Sufficient conditions for the uniform integrability of $\bm{\{\varPsi^{t,x}\}}$

In this section, we extend the results in Henry-Labordère et al (2019) to branching diffusions with jumps to give sufficient conditions for the uniform integrability of $\{\varPsi^{t,x}\}$ as required in Theorem 3.2.

Theorem 3.3 (Integrability conditions).

Let $\kappa\in(1,\infty)$ and define

C_{1}\triangleq\frac{\|g\|_{\infty}^{\kappa}}{F(T)^{\kappa-1}}\quad\text{and}% \quad C_{2}\triangleq\sup_{i\in\mathcal{I},t\in[0,T]}\biggl{(}\frac{\|c_{i}\|_% {\infty}}{\rho(t)p_{i}}\biggr{)}^{\kappa-1}.

Then, the family $\smash{\{\varPsi^{t,x}\}_{(t,x)\in[0,T]\times\mathbb{R}^{d}}}$ is bounded in $\mathcal{L}^{\kappa}(\mathbb{P})$ , and in particular is uniformly integrable, in either of the following two cases.

(i)

It holds that

$\frac{C_{1}}{F(T)}\vee C_{2}\leq 1.$
(ii)

The power series $\sum_{i\in\mathcal{I}}\|c_{i}\|_{\infty}x^{|i|}$ is nonzero⁷⁷ 7 Otherwise PDE (2.1) becomes linear; for that case the branching approach becomes unnecessary. and has an infinite radius of convergence, and the terminal time $T>0$ is sufficiently small in the sense that

$T<\int_{C_{1}}^{\infty}\bigg{(}C_{2}\sum_{i\in\mathcal{I}}\|c_{i}\|_{\infty}x^% {|i|}\bigg{)}^{\!-1}\mathrm{d}x.$ (3.4)

Proof.

Fix some $(t,x)\in[0,T]\times\mathbb{R}^{d}$ . By definition of $\varPsi^{t,x}$ , we have

|\varPsi^{t,x}|^{\kappa}=\prod_{k\in\mathcal{K}}\frac{|g(X_{T}^{(k)})|^{\kappa% }}{F(T-T^{(k-)})^{\kappa}}\prod_{k\in\bar{\mathcal{K}}\setminus\mathcal{K}}% \frac{|c_{I^{(k)}}(T^{(k)},X_{T^{(k)}}^{(k)},\Delta^{(k)})|^{\kappa}}{|\rho(T^% {(k)}-T^{(k-)})p_{I^{(k)}}|^{\kappa}}.

(3.5)

With this, under condition (i), and since $F$ is decreasing, we immediately find that

\mathbb{E}[|\varPsi^{t,x}|^{\kappa}]\leq\mathbb{E}\bigg{[}\prod_{k\in\mathcal{% K}}\frac{C_{1}}{F(T)}\prod_{k\in\bar{\mathcal{K}}\setminus\mathcal{K}}C_{2}^{% \kappa/(\kappa-1)}\bigg{]}\leq\mathbb{E}\bigg{[}\prod_{k\in\bar{\mathcal{K}}}% \frac{C_{1}}{F(T)}\vee C_{2}^{\kappa/(\kappa-1)}\bigg{]}\leq 1

and the proof is complete. Let us therefore subsequently assume that we are in case (ii); we follow the approach of Henry-Labordère et al (2019, Theorem 3.12). The basic idea is to identify an upper bound $\chi_{\infty}$ for $|\varPsi^{t,x}|^{\kappa}$ that can itself be regarded as a branching estimator for an ordinary differential equation (ODE) admitting a global solution.

First note that

\mathbb{E}[|\varPsi^{t,x}|^{\kappa}]\leq\mathbb{E}\bigg{[}\prod_{k\in\mathcal{% K}}\frac{C_{1}}{F(T-T_{t}^{(k-)})}\prod_{k\in\bar{\mathcal{K}}\setminus% \mathcal{K}}\frac{C_{2}\|c_{I^{(k)}}\|_{\infty}}{\rho(T^{(k)}-T^{(k-)})p_{I^{(% k)}}}\bigg{]}.

(3.6)

Now consider the following ODE to be solved backward in time:

\dot{\eta}(t)+C_{2}\sum_{i\in\mathcal{I}}\|c_{i}\|_{\infty}\eta(t)^{|i|}=0,% \quad t\in[0,T],\qquad\eta(T)=C_{1}.

(3.7)

Define the map

\varphi\colon[C_{1},\infty)\to[0,\infty],\quad y\mapsto\varphi(y)\triangleq% \int_{C_{1}}^{y}\bigg{(}C_{2}\sum_{i\in\mathcal{I}}\|c_{i}\|_{\infty}x^{|i|}% \bigg{)}^{\!-1}\mathrm{d}x.

Since the power series is nondegenerate, $\varphi$ is a continuous mapping that is strictly increasing where finite. Upon rearranging and integrating the ODE, note that $\eta\in\mathcal{C}^{1}([0,T])$ is a solution of (3.7) if and only if

-\int_{t}^{T}\bigg{(}C_{2}\sum_{i\in\mathcal{I}}\|c_{i}\|_{\infty}\eta(t)^{|i|% }\bigg{)}^{\!-1}\dot{\eta}(s)\mathrm{d}{s}=\int_{t}^{T}1\mathrm{d}s,\quad t\in% [0,T],\qquad\eta(T)=C_{1},

and the substitution $x\triangleq\eta(s)$ shows that this is the case if and only if

\varphi(\eta(t))=\int_{C_{1}}^{\eta(t)}\bigg{(}C_{2}\sum_{i\in\mathcal{I}}\|c_% {i}\|_{\infty}x^{|i|}\bigg{)}^{\!-1}\mathrm{d}x=T-t,\quad t\in[0,T].

Since $\varphi$ is strictly increasing where it is finite, the latter statement is equivalent to

T\in\operatorname{range}(\varphi)\cap\mathbb{R}=\bigg{\{}\xi\in\mathbb{R}% \colon 0\leq\xi\leq\int_{C_{1}}^{\infty}\bigg{(}C_{2}\sum_{i\in\mathcal{I}}\|c% _{i}\|_{\infty}z^{|i|}\bigg{)}^{\!-1}\mathrm{d}{z}\bigg{\}}.

Thus, condition (ii) on $T$ is both necessary and sufficient to guarantee the existence of a solution $\eta$ of ODE (3.7) on $[0,T]$ . With this, we are now able to define

\displaystyle\chi_{n}

\displaystyle\triangleq\prod_{k\in\bigcup_{\nu=1}^{n}\mathcal{K}^{\nu}}\frac{C% _{1}}{F(T-T^{(k-)})}\prod_{k\in\bigcup_{\nu=1}^{n}\bar{\mathcal{K}}^{\nu}% \setminus\mathcal{K}^{\nu}}\frac{C_{2}\|c_{I^{(k)}}\|_{\infty}}{\rho(T^{(k)}-T% ^{(k-)})p_{I^{(k)}}}\prod_{k\in\bar{\mathcal{K}}^{n+1}}\eta(T^{(k-)})

as well as

\chi_{\infty}\triangleq\lim_{n\to\infty}\chi_{n}=\prod_{k\in\mathcal{K}}\frac{% C_{1}}{F(T-T^{(k-)})}\prod_{k\in\bar{\mathcal{K}}\setminus\mathcal{K}}\frac{C_% {2}\|c_{I^{(k)}}\|_{\infty}}{\rho(T^{(k)}-T^{(k-)})p_{I^{(k)}}}.

By analogous arguments to those in Remark 2.5, we obtain

\displaystyle\eta(t)=\eta(T)+\int_{t}^{T}C_{2}\sum_{i\in\mathcal{I}}\|c_{i}\|_% {\infty}\eta(s)^{|i|}\mathrm{d}s=\mathbb{E}[\chi_{1}]=\cdots=\mathbb{E}[\chi_{% n}],\quad t\in[0,T],~{}n\in\mathbb{N}.

Then, thanks to (3.6) and Fatou’s lemma, it follows that

\mathbb{E}[|\varPsi^{t,x}|^{\kappa}]\leq\mathbb{E}[\chi_{\infty}]\leq\liminf_{% n\to\infty}\mathbb{E}[\chi_{n}]=\eta(t)\leq\sup_{t\in[0,T]}\eta(t)<\infty,

and the proof is complete. ∎

Under the conditions of Theorem 3.3, it follows in particular that

\{\varPsi^{t,x}\}_{(t,x)\in[0,T]\times\mathbb{R}^{d}}

is bounded in $\smash{\mathcal{L}^{1}(\mathbb{P})}$ , so the function $\smash{u\colon[0,T]\times\mathbb{R}^{d}\to\mathbb{R}}$ , $\smash{u(t,x)\triangleq\mathbb{E}[\varPsi^{t,x}]}$ is in fact bounded. Taking $\kappa=2$ , it follows that the branching estimator $\varPsi^{t,x}$ has finite variance; moreover, the second part of the proof of Theorem 3.3 suggests the following choice of branching parameters.

Remark 3.4 (Optimized branching parameters).

Suppose that the particles’ lifetimes are $\operatorname{Exp}(\beta)$ -distributed for some $\beta>0$ , and note that for $\kappa=2$ it follows from (3.5) that

|\varPsi^{t,x}|^{2}\leq\prod_{k\in\mathcal{K}}\frac{\mathrm{e}^{\beta T}\|g\|_% {\infty}^{2}}{F(T-T^{(k-)})}\prod_{k\in\bar{\mathcal{K}}\setminus\mathcal{K}}% \frac{\mathrm{e}^{\beta T}\|c_{I^{(k)}}\|_{\infty}^{2}/\beta p_{I^{(k)}}}{\rho% (T^{(k)}-T^{(k-)})p_{I^{(k)}}}.

(3.8)

As in (3.6) and (3.7), under suitable integrability conditions (see in particular (3.4)), the right-hand side of (3.8) constitutes a branching estimator for the ODE

\dot{\eta}(t)+\sum_{i\in\mathcal{I}}\frac{\|c_{i}\|_{\infty}^{2}}{\beta\mathrm% {e}^{-\beta T}p_{i}}\eta(t)^{|i|}=0,\quad t\in[0,T],\qquad\eta(T)=\mathrm{e}^{% \beta T}\|g\|_{\infty}^{2}.

(3.9)

Thus, we have $\smash{\mathbb{E}[|\varPsi^{t,x}|^{2}]\leq\eta(t)}$ ; in order to minimize this upper bound, we follow Henry-Labordère (2012, Section 4.3) and minimize the absolute value of the terminal slope in (3.9), ie, we seek

\min\sum_{i\in\mathcal{I}}\frac{\mathrm{e}^{\beta T(|i|+1)}\|c_{i}\|_{\infty}^% {2}\|g\|_{\infty}^{2|i|}}{\beta p_{i}}

such that

\beta>0\quad\text{and}\quad p_{i}\geq 0,~{}i\in\mathcal{I},\quad\sum_{i\in% \mathcal{I}}p_{i}=1.

The resulting optimized branching parameters $\{p_{i}^{\star}\}_{i\in\mathcal{I}}$ and $\beta^{\star}$ are thus given by

	$\displaystyle p_{i}^{\star}=\frac{\mathrm{e}^{(1/2)\beta^{\star}T(\|i\|+1)}\\|c_{% i}\\|_{\infty}\\|g\\|_{\infty}^{\|i\|}}{\sum_{j\in\mathcal{I}}\mathrm{e}^{(1/2)% \beta^{\star}T(\|j\|+1)}\\|c_{j}\\|_{\infty}\\|g\\|_{\infty}^{\|j\|}},\quad i\in% \mathcal{I},$
and
	$\displaystyle\beta^{\star}T=\frac{1}{1+\sum_{i\in\mathcal{I}}\|i\|p_{i}^{\star}}.$

4 Monte Carlo simulation: a high-dimensional example

The representation (3.1) derived in Theorem 3.2 makes it possible to compute solutions of PDE (2.1) by direct (nonnested, plain vanilla) Monte Carlo simulation. To illustrate the effectiveness and efficiency of the branching Monte Carlo algorithm, we consider the following stylized problem in $d\geq 1$ dimensions:

	$\displaystyle\partial_{t}u(t,x)+\frac{1}{2d^{2}}\Delta u(t,x)+\int_{\mathbb{R}% ^{d}}\sum_{i\in\mathcal{I}}c_{i}(t,x)u(t,x)^{i_{1}}u(t,x+\xi)^{i_{2}}\gamma(% \mathrm{d}\xi)=0,$
	$\displaystyle u(T,x)=\cos(1_{d}^{\mathrm{T}}x).$

Here, $\Delta$ denotes the Laplace operator in $\mathbb{R}^{d}$ ; the time horizon is $T=1$ ; $\gamma$ is the discrete uniform distribution supported on $\{-(\pi/2)e_{i}\in\mathbb{R}^{d}\colon i\in[1:d]\}$ ;⁸⁸ 8 Here, $e_{i}$ denotes the $i$ th unit vector in $\mathbb{R}^{d}$ . and the set of possible descendants is given by $\mathcal{I}\triangleq\{i\in\mathbb{N}_{0}^{2}\colon|i|\leq 2\}$ with coefficients

	$\displaystyle c_{(0,0)}(t,x)$	$\displaystyle=(\alpha+1/(2d))\cos(1_{d}^{\mathrm{T}}x)\exp\{\alpha(T-t)\}+\cos% (1_{d}^{\mathrm{T}}x)^{2}/d-1/(2d),$
	$\displaystyle c_{(1,0)}(t,x)$	$\displaystyle=(-1)\cos(1_{d}^{\mathrm{T}}x)\exp\{-\alpha(T-t)\}/d,$
	$\displaystyle c_{(0,1)}(t,x)$	$\displaystyle=(-1)c_{(1,0)}(t,x),$
	$\displaystyle c_{(2,0)}(t,x)$	$\displaystyle=c_{(0,2)}(t,x)=\exp\{-2\alpha(T-t)\}/(2d),$
	$\displaystyle c_{(1,1)}(t,x)$	$\displaystyle=(-2)c_{(2,0)}(t,x),$

where $\alpha=0.2$ . It is not hard to verify that the solution is given in closed form by

u(t,x)=\cos(1_{d}^{\mathrm{T}}x)\mathrm{e}^{\alpha(T-t)}=\cos\bigg{(}\sum_{i=1% }^{d}x_{i}\bigg{)}\mathrm{e}^{\alpha(T-t)},\quad(t,x)\in[0,T]\times\mathbb{R}^% {d}.

We refer to $u$ as the exact solution and use it as a benchmark to quantify the error of the estimates produced by our algorithm. All subsequent simulation results correspond to the initial configuration $(t,x)=(0,1_{d})$ , ie, we determine $u^{*}\triangleq u(0,1_{d})$ .

The branching parameters can be specified arbitrarily; as a general rule, see Remark 3.4 and Henry-Labordère (2012, Section 4.3): the aim is to choose the probability weights $p_{i}$ proportional to $\|c_{i}\|_{\infty}\|g\|^{|i|}_{\infty}$ . We adopt this route in Section 5; here, we simply use the ad hoc parameters reported in Table 1, which correspond to the number of descendants being distributed uniformly on $\{0,1,2\}$ , and, given that number, the jump marks are iid Bernoulli $(\frac{1}{2})$ across descendants.

Table 1: Parameters of the branching diffusion.

Parameter	Value
$\mathsf{Law}(\tau^{(\text{1})})$	$\varGamma(\kappa,\theta)$ with $\kappa=\text{0.5}$
	and $\theta=\text{2.5}$
$p_{(\text{0},\text{0})}$	$\text{1}/\text{3}$
$p_{(\text{1},\text{0})}$ , $p_{(\text{0},\text{1})}$ , $p_{(\text{1},\text{1})}$	$\text{1}/\text{6}$
$p_{(\text{2},\text{0})}$ , $p_{(\text{0},\text{2})}$	$\text{1}/\text{12}$

Our simulation study is conducted as follows. Given a spatial dimension $\smash{d\in\mathbb{N}}$ and a number of Monte Carlo simulations $\smash{N}$ , a simulation run consists of computing the estimator $\smash{\hat{u}_{d,N}}$ of $u^{*}$ as the average of $N$ iid copies $\smash{\{\varPsi^{0,1_{d}}_{n}\}_{n\in[1:N]}}$ of $\varPsi^{0,1_{d}}$ . All numerical computations are implemented in MATLAB.⁹⁹ 9 Hardware: Intel Core i7-6700 (3.40 GHz) with 31.2 GiB random access memory; Linux Mint operating system (18.1 Serena); no parallelization. Table 2 presents the simulation results for the estimator $\smash{\hat{u}_{d,N}}$ for several values of $d$ and $N$ .

Table 2: Simulation results for

\hat{u}_{d,N}

. [Standard deviations are given in parentheses.]

$?$	3	5	10	20	50	100
	$?$

$\text{10}^{\text{2}}$	$-$ 1.6463	0.3150	$-$ 0.9049	0.4225	1.2471	1.3105
	$-$ (0.255)	(0.130)	$-$ (0.131)	(0.068)	(0.122)	(0.109)
$\text{10}^{\text{3}}$	$-$ 1.2254	0.3240	$-$ 1.0134	0.4909	1.1319	1.0323
	$-$ (0.124)	(0.043)	$-$ (0.037)	(0.019)	(0.038)	(0.034)
$\text{10}^{\text{4}}$	$-$ 1.2643	0.3573	$-$ 1.0311	0.5013	1.1668	1.0488
	$-$ (0.054)	(0.014)	$-$ (0.012)	(0.006)	(0.012)	(0.011)
$\text{10}^{\text{5}}$	$-$ 1.1758	0.3492	$-$ 1.0260	0.4956	1.1797	1.0549
	$-$ (0.017)	(0.005)	$-$ (0.004)	(0.002)	(0.004)	(0.003)
$\text{10}^{\text{6}}$	$-$ 1.2064	0.3485	$-$ 1.0262	0.4969	1.1791	1.0536
	$-$ (0.007)	(0.001)	$-$ (0.001)	(0.001)	(0.001)	(0.001)
$u^{*}$	$-$ 1.2092	0.3465	$-$ 1.0248	0.4984	1.1786	1.0532

Table 3: Relative error

\text{err}_{\text{rel}}(\hat{u}_{d,N})

. [Standard deviations are given in parentheses. Estimates and standard deviations are based on 100 mutually independent simulation runs of

\hat{u}_{d,N}

. All values are in percent.]

$?$	3	5	10	20	50	100
	$?$

$\text{10}^{\text{2}}$	33.6 (36.1)	30.1 (26.7)	8.9 (7.2)	9.9 (7.8)	7.8 (5.7)	7.9 (6.3)
$\text{10}^{\text{3}}$	13.5 (14.5)	10.9 0(8.8)	3.1 (2.4)	3.3 (2.4)	2.5 (2.0)	2.6 (2.1)
$\text{10}^{\text{4}}$	04.4 0(4.3)	03.5 0(2.6)	0.8 (0.6)	1.0 (0.8)	0.8 (0.7)	0.7 (0.6)
$\text{10}^{\text{5}}$	01.6 0(1.4)	01.0 0(0.8)	0.3 (0.2)	0.3 (0.2)	0.3 (0.2)	0.2 (0.2)
$\text{10}^{\text{6}}$	00.6 0(0.4)	00.3 0(0.3)	0.1 (0.1)	0.1 (0.1)	0.1 (0.1)	0.1 (0.1)

To quantify the accuracy of the Monte Carlo algorithm, we further denote the relative error compared to the exact solution by

\mathrm{err}_{\mathrm{rel}}(\hat{u}_{d,N})\triangleq\frac{|\hat{u}_{d,N}-u^{*}% |}{|u^{*}|}

and report our simulation results in Table 3.

Table 4: Running time (in seconds) to compute

\hat{u}_{d,N}

. [Standard deviations are given in parentheses. Estimates and standard deviations are based on 100 mutually independent simulation runs of

\hat{u}_{d,N}

.]

	$?$

$?$	3	5	10	20	50	100
$\text{10}^{\text{2}}$	000.1 (0.0)	000.0 (0.0)	000.0 (0.0)	000.0 (0.0)	000.0 (0.0)	000.0 0(0.0)
$\text{10}^{\text{3}}$	000.1 (0.0)	000.2 (0.0)	000.2 (0.0)	000.1 (0.0)	000.2 (0.0)	000.2 0(0.0)
$\text{10}^{\text{4}}$	001.4 (0.1)	001.5 (0.1)	001.5 (0.1)	001.5 (0.1)	001.6 (0.1)	003.2 0(0.1)
$\text{10}^{\text{5}}$	014.3 (0.5)	014.4 (0.3)	014.4 (0.3)	014.6 (0.3)	015.2 (0.3)	016.7 0(0.3)
$\text{10}^{\text{6}}$	140.8 (1.1)	140.5 (1.4)	141.2 (1.3)	142.5 (1.0)	162.2 (4.6)	185.5 (14.9)

The algorithm is able to achieve high levels of accuracy even for high dimensions of $d=50$ and $d=100$ , provided $N$ is sufficiently large. Note that in this example the precision of the results actually increases with the dimension $d$ ; this is in line with Theorem 3.3, as in the above specification the norm $\|c_{i}\|_{\infty}$ appearing in the constant $C_{2}$ decreases with $d$ .

At the same time, the running times required to achieve these levels of accuracy are rather modest; the exact values are reported in Table 4.

Finally, Figure 2 displays the relative error, the standard deviation and the running time as functions of the number of Monte Carlo samples $N$ . As expected from the central limit theorem, the slope in the logarithmic plot of the standard deviation is approximately $-\frac{1}{2}$ . The running times also demonstrate that there is no curse-of-dimensionality effect.

5 Valuation with systemically important counterparties

In this section we illustrate the usefulness of the theory developed in this paper in the valuation of financial positions with systemic counterparty credit risk. Specifically, we assume that the counterparty in a given financial position is a systemically important bank (SIB);¹⁰¹⁰ 10 A list of banks classified as globally systemically important by the Financial Stability Board is available at https://bit.ly/3nRHNWm. its systemic importance is captured by jumps in the underlying risk factors, or equivalently, devaluations in risky asset prices, that occur upon the SIB’s default. This model setup was first proposed by Pykhtin and Sokol (2013) (see also Brigo et al 2019; Mercurio and Li 2015). We wish to stress that our focus here is on situations where finite-difference or fixed point methods for the pricing PDE (as developed by, for example, Carr and Ghamami (2017) or Kim and Leung (2016)) are not applicable.

5.1 Valuation with systemic risk

We consider a financial market that is free of arbitrage; the underlying probability measure (denoted by $\mathbb{P}$ in the preceding sections) is taken to be the relevant risk-neutral pricing measure $\mathbb{Q}$ . The financial market consists of a locally risk-free money market account $\smash{B=\{B_{t}\}_{t\in[0,T]}}$ and $d\in\mathbb{N}$ dynamically traded risky assets with prices $\smash{X=\{X_{t}\}_{t\in[0,T]}}$ given by an $\smash{\mathbb{R}^{d}_{+}}$ -valued semimartingale such that $X/B$ is a local $\mathbb{Q}$ -martingale. The financial position whose price is to be determined promises a time- $T$ payoff $g(X_{T})$ (where $g$ is a bounded measurable function of the underlyings $X$ ), provided the counterparty does not default before time $T$ .

Crucially, the counterparty in this financial position is a defaultable SIB. This means that, first, the financial position is subject to credit risk, and second, if and when the SIB counterparty defaults, there is a negative impact on risky asset prices $X$ , which simultaneously affects the value of the financial position $g(X_{T})$ . Assuming a fractional recovery of postdefault mark-to-market value as in Duffie and Singleton (1999), risk-neutral pricing yields

\frac{V_{t}}{B_{t}}=\mathbb{E}^{\mathbb{Q}}_{t}\bigg{[}\bm{1}_{\{\tau>T\}}% \frac{g(X_{T})}{B_{T}}+\bm{1}_{\{\tau\leq T\}}\frac{h(V_{\tau})}{B_{\tau}}% \bigg{]}\quad\text{on }\{\tau>t\},

(5.1)

where $V$ represents the value of the financial position provided the counterparty has not defaulted; $\tau$ is the (original) counterparty’s default time; and the recovery value the investor retrieves is given by $h(V_{\tau})$ , where $h(v)\triangleq Rv^{+}-v^{-}$ with a recovery rate $R\in[0,1]$ . Note that in (5.1) $V_{\tau}$ represents the time- $\tau$ mark-to-market price of an identical financial position with an SIB counterparty that has not defaulted, but is otherwise identical to the original one, immediately after the original counterparty’s default. In particular, $V_{\tau}$ is based on postdefault risky asset prices $X_{\tau}=X_{\tau-}+\Delta X_{\tau}$ , and is itself determined endogeneously by (5.1); in particular, a direct Monte Carlo approach to (5.1) would require nested simulations and thus become computationally infeasible.

Specifically, we assume that default events are modeled within a classical reduced-form framework and that the risk-neutral dynamics of $X$ are given by an $\mathbb{R}^{d}_{+}$ -valued jump-diffusion. Thus, the SIB counterparty’s default time is given by

\tau\triangleq\inf\{t\geq 0\colon Y_{t}\not=0\},

where $\smash{Y=\{Y_{t}\}_{t\in[0,T]}}$ is a Cox process with intensity $\smash{\{\lambda(t,X_{t})\}_{t\in[0,T]}}$ , and the risk-neutral dynamics of $B$ and $X$ are given by

	$\displaystyle\mathrm{d}B_{t}$	$\displaystyle=r(t,X_{t})B_{t}\mathrm{d}t,$
	$\displaystyle\mathrm{d}X_{t}$	$\displaystyle=\operatorname{diag}(X_{t-})[\mu(t,X_{t-})\mathrm{d}t+\sigma(t,X_% {t-})\mathrm{d}W_{t}+\Delta_{Y_{t}}\mathrm{d}Y_{t}],\quad X_{0}=x,$

where $\smash{\{\Delta_{n}\}_{n\in\mathbb{N}}}$ are iid $\smash{E\triangleq(-1,\infty)^{d}}$ -valued and independent of $Y$ and $W$ and

\mu(t,x)\triangleq r(t,x)1_{d}-\lambda(t,x)\mathbb{E}[\Delta_{1}].

In particular, the SIB counterparty’s default triggers a simultaneous devaluation in the financial position’s underlyings of size $\Delta_{1}$ . The corresponding pricing PDE is given by

	$\displaystyle\partial_{t}u(t,x)+\mathcal{A}^{\star}[u](t,x)-r(t,x)u(t,x)+% \lambda(t,x)$
	$\displaystyle\qquad\qquad\qquad\qquad\times\int_{E}(h\circ u(t,x+\operatorname% {diag}(\xi)x)-u(t,x))\nu(\mathrm{d}\xi)=0$

subject to the boundary condition $u(T,x)=g(x)$ , where $\nu$ denotes the distribution of $\Delta_{n}$ , $n\in\mathbb{N}$ , and the operator $\mathcal{A}^{\star}$ is given by

	$\displaystyle\mathcal{A}^{\star}[u](t,x)$	$\displaystyle\triangleq\mu(t,x)^{\mathrm{T}}\operatorname{diag}(x)\mathrm{D}_{% x}u(t,x)$
		$\displaystyle\qquad+\tfrac{1}{2}\operatorname{tr}[\operatorname{diag}(x)\sigma% (t,x)\sigma(t,x)^{\mathrm{T}}\operatorname{diag}(x)\mathrm{D}_{x}^{2}u(t,x)].$

Note that, similarly as in other credit risk valuation problems with recovery of mark-to-market value, the pricing formula (5.1) and the corresponding pricing equation are inherently implicit, with the price to be determined appearing inside the nonlinearity that represents the recovery value (see, for example, Carr and Ghamami (2017) or Kim and Leung (2016) and the references therein), with the additional complication that jumps at default imply that this term is also nonlocal.

5.2 Branching diffusion with jumps approach

Approximation of h by a polynomial. — Figure 3: Approximation of $h$ by a polynomial.

In order to solve the pricing equation, we follow Henry-Labordère (2012, Section 5) and approximate the recovery function $h$ by a polynomial (see Figure 3). Thus, we assume without loss of generality that $\|g\|_{\mathcal{L}^{\infty}}\leq 1$ (since $g$ is bounded, this can always be achieved by appropriate rescaling) and approximate the recovery value function using

v^{\pm}\approx\sum_{m=0}^{M}\alpha_{m}^{\pm}v^{m},\quad v\in[-1,1];

we then obtain

	$\displaystyle\partial_{t}u(t,x)+\mathcal{A}^{\star}[u](t,x)-[r(t,x)+\lambda(t,% x)]u(t,x)$
	$\displaystyle\qquad\qquad+\int_{E}\sum_{m=0}^{M}[R\alpha_{m}^{+}-\alpha_{m}^{-% }]\lambda(t,x)u(t,x+\operatorname{diag}(\xi)x)^{m}\nu(\mathrm{d}\xi)=0.$		(5.2)

Hence, setting $\mathcal{I}\triangleq\{(\ell,m)\in[0:M]^{2}\colon\ell=0\text{ or }m=0\}$ we can rewrite (5.2) as

	$\displaystyle\partial_{t}u(t,x)+\mathcal{A}^{\star}[u](t,x)+\int_{E}\sum_{i\in% \mathcal{I}}c_{i}(t,x)u(t,x)^{i_{1}}u(t,x+\operatorname{diag}(\xi)x)^{i_{2}}% \nu(\mathrm{d}\xi)=0,$		(5.3)
	$\displaystyle u(T,x)=g(x),$

where the coefficients $c_{(\ell,m)}$ , $(\ell,m)\in\mathcal{I}$ , are given by

	$\displaystyle c_{(0,m)}(t,x)$	$\displaystyle\triangleq[R\alpha_{m}^{+}-\alpha_{m}^{-}]\lambda(t,x)$		$\displaystyle\text{for }m\in[0:M],$
	$\displaystyle c_{(1,0)}(t,x)$	$\displaystyle\triangleq-[r(t,x)+\lambda(t,x)]\quad\text{and}\quad c_{(\ell,0)}% (t,x)\triangleq 0$		$\displaystyle\text{for }\ell\in[2:M].$

Since (5.3) is a special case of PDE (2.1), the branching diffusion with jumps approach developed in this paper can be applied to compute the solution via $u(0,x)=\mathbb{E}[\varPsi^{0,x}]$ with

	$\displaystyle\varPsi^{0,x}$	$\displaystyle\triangleq\prod_{k\in\mathcal{K}}\frac{g(X^{(k)}_{T})}{F(T-T^{(k-% )})}\times\prod_{k\in\bar{\mathcal{K}}\setminus\mathcal{K},I^{(k)}=(1,0)}-% \frac{r(T^{(k)},X^{(k)}_{T^{(k)}})+\lambda(T^{(k)},X^{(k)}_{T^{(k)}})}{\rho(T^% {(k)}-T^{(k-)})p_{(1,0)}}$
		$\displaystyle\qquad\times\prod_{m=0}^{M}\prod_{k\in\bar{\mathcal{K}}\setminus% \mathcal{K},I^{(k)}=(0,m)}\frac{[R\alpha_{m}^{+}-\alpha_{m}^{-}]\lambda(T^{(k)% },X^{(k)}_{T^{(k)}})}{\rho(T^{(k)}-T^{(k-)})p_{(0,m)}}.$

Here, $\smash{X^{(k)}\triangleq X^{k,0,x}}$ , $\smash{k\in\bar{\mathcal{K}}\triangleq\bar{\mathcal{K}}^{0}}$ , is a branching diffusion with jumps as specified in Section 2, where the dynamics of each individual particle $\smash{k\in\bar{\mathcal{K}}}$ are given by

\mathrm{d}X^{(k)}_{t}=\operatorname{diag}(X^{(k)}_{t-})[\mu(t,X^{(k)}_{t-})% \mathrm{d}t+\sigma(t,X^{(k)}_{t-})\mathrm{d}W_{t}],\quad t\in[T^{(k-)},T^{(k)}],

with initial conditions $X^{(1)}_{0}=x$ and

X^{(k)}_{T^{(k-)}}\triangleq\begin{cases}X^{(k-)}_{T^{(k-)}}&\text{if }I^{(k)}% =(1,0),\\ \\ X^{(k-)}_{T^{(k-)}}+\operatorname{diag}(\Delta^{(k-)})X^{(k-)}_{T^{(k-)}}&% \text{otherwise.}\end{cases}

5.3 Numerical results

While our setup and Monte Carlo methodology allow for general diffusion dynamics in the underlying risky assets, for illustration we use a multivariate Black–Scholes–Merton model, enhanced by the SIB’s credit risk with constant default intensity and devaluations captured by the model of Kou (2002). Thus, we take $Y$ as a Poisson process with intensity $\lambda\geq 0$ ; the risk-free rate $r\in\mathbb{R}$ is constant; and $g$ has $d=499$ underlying assets with volatility matrix $\sigma\in\mathbb{R}^{d\times d}$ calibrated to the S&P 500 stock index.¹¹¹¹ 11 The calibration period is January 1, 2019 to December 31, 2019; of the 505 constituents, six had to be omitted due to a lack of data. Jumps at default affect all assets equally, ie, $\Delta_{n}=\delta_{n}1_{d}$ , where $-\log(1+\delta_{n})$ is exponentially distributed with parameter $\eta\geq 0$ . We consider two financial positions representing a shifted equally weighted basket put and a shifted equally weighted basket discount call, respectively, ie,

g_{\mathrm{short}}(X_{T})\triangleq\bigg{(}K-\frac{1}{d}\sum_{i=1}^{d}X^{i}_{T% }\bigg{)}^{\!+}-L,\qquad g_{\mathrm{long}}(X_{T})=L-\bigg{(}K-\frac{1}{d}\sum_% {i=1}^{d}X^{i}_{T}\bigg{)}^{\!+}.

Figure 4 provides an illustration of the corresponding payoffs.

Payoff profiles ... and ... with K=2 and L=1. — Figure 4: Payoff profiles $g_{\text{short}}$ and $g_{\text{long}}$ with $K=\text{2}$ and $L=\text{1}$ .

To quantify the impact of systemic risk, we compare the valuations of the financial positions $g_{\mathrm{short}}$ and $g_{\mathrm{long}}$ in three benchmark scenarios that are identical in all respects except the choice of counterparty.

SIB counterparty.: The counterparty is systemically important, and its default triggers devaluations in risky asset prices (see above).
Non-SIB counterparty.: The counterparty is nonsystemic, but otherwise identical to the SIB; in particular, it is defaultable with the same credit risk characteristics as the SIB.
Default-free counterparty.: There is no counterparty credit risk.

In all three scenarios, the SIB is part of the model and will cause devaluations upon default; the scenarios thus differ only in the choice of counterparty and the resulting wrong- or right-way risk. The implementation of the branching diffusion with jumps is based on the parameters specified in Table 6, and uses the optimized branching parameters from Remark 3.4. In contrast to Section 4, where we deliberately employ a nonaccelerated standard Monte Carlo method for performance analysis, here we exploit standard techniques for variance reduction such as control variates and parallelization (see Glasserman 2003; Korn et al 2010).¹²¹² 12 The control variate is the default-free price of the corresponding geometric basket derivative. The relevant model and simulation parameters are reported in Tables 5 and 6; note in particular that for each underlying asset the expected devaluation upon the SIB’s default is $-50\%$ .

Figures 5–7 display our simulation results for different specifications of the default intensity.

Valuation of ... as a function of the default intensity (with 99%-confidence bands). — Figure 5: Valuation of $g_{\text{short}}$ as a function of the default intensity (with 99%-confidence bands).

Table 5: Market coefficients.

$?$	$?$	$?$	$\eta$	$x^{i}$	$?$	$?$
1	0.5%	40%	1	1	2	1

Table 6: Simulation parameters.

$?$	$?$	$\alpha^{\pm}_{\textbf{0}}$	$\alpha^{\pm}_{\textbf{1}}$	$\alpha^{\pm}_{\textbf{2}}$	$\alpha^{\pm}_{\textbf{3}}$	$\alpha^{\pm}_{\textbf{4}}$
$\text{5}\times\text{10}^{\text{6}}$	4	$-$ 0.06	$\pm$ 0.50	$-$ 0.82	$-$ 0.00	$-$ 0.41

Figure 6: Decomposition of the SIB spread into credit risk and systemic risk.

Figure 5 illustrates the impact of systemic interaction for the financial position $g_{\mathrm{short}}$ : devaluations imply a positive correlation between the SIB’s default events and underperforming risky asset prices, causing significant wrong-way risk for short positions. While this is qualitatively apparent, the quantitative size of this effect, particularly relative to the non-SIB counterparty, is remarkable.

In Figure 6 we thus decompose the implied SIB spread (ie, the difference between the value of an otherwise identical financial position with a default-free counterparty and the value of one with an SIB counterparty) into a pure credit risk component (green), which we identify with the spread between the nonsystemic defaultable counterparty and the default-free one, and a systemic risk component (red), which is present only because of the counterparty’s systemic importance. It is apparent that systemic risk is the main driver of the spread, accounting for 85–90% of the total spread across realistic default intensities.

Finally, Figure 7 demonstrates that the effect is reversed for long positions; we might refer to this as right-way risk, ie, negative correlation between the counterparty’s default and the value of the financial position. Technically, the systemic risk component turns negative, and the spread becomes significantly smaller than the spread for a nonsystemic counterparty.

6 Conclusion

To conclude, in this paper we developed a branching diffusion with jumps approach to solving parabolic PDEs with nonlocal nonlinearities. We showcased the performance of the resulting nonnested Monte Carlo methodology, and demonstrated how it applies to the valuation of financial positions with systemic counterparties and mark-to-market recovery. Several extensions to the above pricing model are conceivable within the setup of this paper, including stochastic recovery rates, shocks to state variables (eg, volatilities), idiosyncratic jumps in asset prices and margin periods of risk.

Declaration of interest

The authors report no conflicts of interest. The authors alone are responsible for the content and writing of the paper.

Acknowledgements

The authors thank one of the journal’s associate editors and the referees for valuable comments and suggestions. Daniel Hoffmann and Frank Seifried gratefully acknowledge financial support from Deutsche Forschungsgemeinschaft (DFG) as part of the Research Training Group 2126 (Algorithmic Optimization). We also thank Tuomo Lampinen (Silicon Cloud Technologies, LLC) for kindly providing historical S&P 500 data.

References

Agarwal, A., and Claisse, J. (2020). Branching diffusion representation of semi-linear elliptic PDEs and estimation using Monte Carlo method. Stochastic Processes and Their Applications 130(8), 5006–5036 (https://doi.org/10.1016/j.spa.2020.02.009).
Athreya, K. B., and Ney, P. E. (1972). Branching Processes. Springer (https://doi.org/10.1007/978-3-642-65371-1).
Bouchard, B., Tan, X., Warin, X., and Zou, Y. (2017). Numerical approximation of BSDEs using local polynomial drivers and branching processes. Monte Carlo Methods and Applications 23(4), 241–263 (https://doi.org/10.1515/mcma-2017-0116).
Bouchard, B., Tan, X., and Warin, X. (2019). Numerical approximation of general Lipschitz BSDEs with branching processes. ESAIM: Proceedings and Surveys 65, 309–329 (https://doi.org/10.1051/proc/201965309).
Brigo, D., Pede, N., and Petrelli, A. (2019). Multi currency credit default swaps: quanto effects and FX devaluation jumps. International Journal of Theoretical and Applied Finance 22(4), Paper 1950018 (https://doi.org/10.1142/S0219024919500183).
Carr, P., and Ghamami, S. (2017). Derivatives pricing under bilateral counterparty risk. The Journal of Risk 20(1), 77–107 (https://doi.org/10.21314/JOR.2017.395).
Duffie, D., and Singleton, K. (1999). Modeling term structure of defaultable bonds. Review of Financial Studies 12(4), 687–720 (https://doi.org/10.1093/rfs/12.4.687).
Glasserman, P. (2003). Monte Carlo Methods in Financial Engineering. Springer (https://doi.org/10.1007/978-0-387-21617-1).
Harris, T. E. (1963). The Theory of Branching Processes. Springer.
Henry-Labordère, P. (2012). Counterparty risk valuation: a marked branching diffusion approach. Risk 25(7), 133–137.
Henry-Labordère, P., and Touzi, N. (2018). Branching diffusion representation for nonlinear Cauchy problems and Monte Carlo approximation. Annals of Applied Probability 31(5), 2350–2375 (https://doi.org/10.1214/20-AAP1649).
Henry-Labordère, P., Tan, X., and Touzi, N. (2014). A numerical algorithm for a class of BSDEs via the branching process. Stochastic Processes and Their Applications 124(2), 1112–1140 (https://doi.org/10.1016/j.spa.2013.10.005).
Henry-Labordère, P., Oudjane, N., Tan, X., Touzi, N., and Warin, X. (2019). Branching diffusion representation of semilinear PDEs and Monte Carlo approximation. Annales de l’Institut Henri Poincaré: Probabilités et Statistiques 55(1), 184–210 (https://doi.org/10.1214/17-AIHP880).
Karatzas, I., and Shreve, S. (1998). Brownian Motion and Stochastic Calculus, 2nd edn. Springer (https://doi.org/10.1007/978-1-4612-0949-2).
Kim, J., and Leung, T. (2016). Pricing derivatives with counterparty risk and collateralization: a fixed point approach. European Journal of Operational Research 249(2), 525–539 (https://doi.org/10.1016/j.ejor.2015.06.055).
Korn, R., Korn, E., and Kroisandt, G. (2010). Monte Carlo Methods and Models in Finance and Insurance. CRC Press, Boca Raton, FL (https://doi.org/10.1201/9781420076196).
Kou, S. G. (2002). A jump-diffusion model for option pricing. Management Science 48(8), 1086–1101 (https://doi.org/10.1287/mnsc.48.8.1086.166).
Mercurio, F., and Li, M. (2015). Jumping with default: wrong-way risk modeling for CVA. Risk 28(12), 58–63.
Pardoux, E., and Răşcanu, A. (2014). Stochastic Differential Equations, Backward SDEs, Partial Differential Equations. Springer (https://doi.org/10.1007/978-3-319-05714-9).
Pykhtin, M., and Sokol, A. (2013). Exposure under systemic risk impact. Risk 26(9), 88–93.

[bib.bib1] Agarwal, A., and Claisse, J. (2020). Branching diffusion representation of semi-linear elliptic PDEs and estimation using Monte Carlo method. Stochastic Processes and Their Applications 130(8), 5006–5036 (https://doi.org/10.1016/j.spa.2020.02.009).

[bib.bib2] Athreya, K. B., and Ney, P. E. (1972). Branching Processes. Springer (https://doi.org/10.1007/978-3-642-65371-1).

[bib.bib3] Bouchard, B., Tan, X., Warin, X., and Zou, Y. (2017). Numerical approximation of BSDEs using local polynomial drivers and branching processes. Monte Carlo Methods and Applications 23(4), 241–263 (https://doi.org/10.1515/mcma-2017-0116).

[bib.bib4] Bouchard, B., Tan, X., and Warin, X. (2019). Numerical approximation of general Lipschitz BSDEs with branching processes. ESAIM: Proceedings and Surveys 65, 309–329 (https://doi.org/10.1051/proc/201965309).

[bib.bib5] Brigo, D., Pede, N., and Petrelli, A. (2019). Multi currency credit default swaps: quanto effects and FX devaluation jumps. International Journal of Theoretical and Applied Finance 22(4), Paper 1950018 (https://doi.org/10.1142/S0219024919500183).

[bib.bib6] Carr, P., and Ghamami, S. (2017). Derivatives pricing under bilateral counterparty risk. The Journal of Risk 20(1), 77–107 (https://doi.org/10.21314/JOR.2017.395).

[bib.bib7] Duffie, D., and Singleton, K. (1999). Modeling term structure of defaultable bonds. Review of Financial Studies 12(4), 687–720 (https://doi.org/10.1093/rfs/12.4.687).

[bib.bib8] Glasserman, P. (2003). Monte Carlo Methods in Financial Engineering. Springer (https://doi.org/10.1007/978-0-387-21617-1).

[bib.bib9] Harris, T. E. (1963). The Theory of Branching Processes. Springer.

[bib.bib10] Henry-Labordère, P. (2012). Counterparty risk valuation: a marked branching diffusion approach. Risk 25(7), 133–137.

[bib.bib11] Henry-Labordère, P., and Touzi, N. (2018). Branching diffusion representation for nonlinear Cauchy problems and Monte Carlo approximation. Annals of Applied Probability 31(5), 2350–2375 (https://doi.org/10.1214/20-AAP1649).

[bib.bib12] Henry-Labordère, P., Tan, X., and Touzi, N. (2014). A numerical algorithm for a class of BSDEs via the branching process. Stochastic Processes and Their Applications 124(2), 1112–1140 (https://doi.org/10.1016/j.spa.2013.10.005).

[bib.bib13] Henry-Labordère, P., Oudjane, N., Tan, X., Touzi, N., and Warin, X. (2019). Branching diffusion representation of semilinear PDEs and Monte Carlo approximation. Annales de l’Institut Henri Poincaré: Probabilités et Statistiques 55(1), 184–210 (https://doi.org/10.1214/17-AIHP880).

[bib.bib14] Karatzas, I., and Shreve, S. (1998). Brownian Motion and Stochastic Calculus, 2nd edn. Springer (https://doi.org/10.1007/978-1-4612-0949-2).

[bib.bib15] Kim, J., and Leung, T. (2016). Pricing derivatives with counterparty risk and collateralization: a fixed point approach. European Journal of Operational Research 249(2), 525–539 (https://doi.org/10.1016/j.ejor.2015.06.055).

[bib.bib16] Korn, R., Korn, E., and Kroisandt, G. (2010). Monte Carlo Methods and Models in Finance and Insurance. CRC Press, Boca Raton, FL (https://doi.org/10.1201/9781420076196).

[bib.bib17] Kou, S. G. (2002). A jump-diffusion model for option pricing. Management Science 48(8), 1086–1101 (https://doi.org/10.1287/mnsc.48.8.1086.166).

[bib.bib18] Mercurio, F., and Li, M. (2015). Jumping with default: wrong-way risk modeling for CVA. Risk 28(12), 58–63.

[bib.bib19] Pardoux, E., and Răşcanu, A. (2014). Stochastic Differential Equations, Backward SDEs, Partial Differential Equations. Springer (https://doi.org/10.1007/978-3-319-05714-9).

[bib.bib20] Pykhtin, M., and Sokol, A. (2013). Exposure under systemic risk impact. Risk 26(9), 88–93.

Journal of Computational Finance

1 Introduction

2 Branching diffusion representations for nonlocal partial differential equations

2.1 Preliminaries

2.2 Branching mechanism

2.3 Branching diffusion dynamics

Assumption 2.1.

Lemma 2.2 (Properties of the diffusion).

2.4 Branching diffusion representation

Assumption 2.3.

Definition 2.4 (Classical solution).

Remark 2.5 (Randomization).

Theorem 2.6.

Proof of Theorem 2.6.

Remark 2.7.

3 Viscosity solutions and branching diffusion representations

3.1 Viscosity solutions of nonlocal nonlinear PDEs

Definition 3.1.

Theorem 3.2 (Viscosity property of the branching representation).

Proof.

3.2 Sufficient conditions for the uniform integrability of {??,?}

Theorem 3.3 (Integrability conditions).

Proof.

Remark 3.4 (Optimized branching parameters).

4 Monte Carlo simulation: a high-dimensional example

5 Valuation with systemically important counterparties

5.1 Valuation with systemic risk

5.2 Branching diffusion with jumps approach

5.3 Numerical results

6 Conclusion

Declaration of interest

Acknowledgements

References

You are currently on corporate access.

3.2 Sufficient conditions for the uniform integrability of $\bm{\{\varPsi^{t,x}\}}$