Monotone methods in counterparty risk models with nonlinear Black–Scholes-type equations

Abstract

A nonlinear Black–Scholes-type equation is studied within counterparty risk models. The classical hypothesis on the uniform Lipschitz-continuity of the nonlinear reaction function allows for an equivalent transformation of the semilinear Black–Scholes equation into a standard parabolic problem with a monotone nonlinear reaction function and an inhomogeneous linear diffusion equation. This setting allows us to construct a scheme of monotone, increasing or decreasing, iterations that converge monotonically to the true solution. As typically any numerical solution of this problem uses most computational power for computing an approximate solution to the inhomogeneous linear diffusion equation, we discuss also this question and suggest several solution methods, including those based on Monte Carlo and finite differences/elements.

1 Introduction

Risk phenomena and their management have been an important topic of investigation since the financial crisis of 2007–2008. In this article we focus our attention on for options with risky values $\hat{V}(S,t)\in {\mathbb {R}}$ modelled by a nonlinear Black–Scholes-type equation:

$$\begin{aligned}&\frac{\partial \hat{V}}{\partial t} + {\mathcal {A}}_t\hat{V} - (r + \lambda _B + \lambda _C) \hat{V} = F(\hat{V}(S,t); S,t) \nonumber \\&\quad \hbox {for }\, (S,t)\in (0,\infty )\times (0,T); \end{aligned}$$

(1.1)

$$\begin{aligned}&\hat{V}(S,T) = h(S) \quad \hbox { for }\, S\in (0,\infty ) . \end{aligned}$$

(1.2)

The nonlinearity, $F(\,\cdot \, ; S,t):\,{\mathbb {R}}\rightarrow {\mathbb {R}},$ with ${\mathbb {R}}= (-\infty , +\infty )$ standing for the real line, is given by

$$\begin{aligned} F(M; S,t){\mathop {=}\limits ^{\mathrm{{def}}}}{}=\,&(R_B\lambda _B + \lambda _C)\, M^{-} - (\lambda _B + R_C\lambda _C)\, M^{+} + s_F\, M^{+}\nonumber \\&\quad \hbox {for } M\in {\mathbb {R}}\hbox { and }\, (S,t)\in (0,\infty )\times (0,T) , \end{aligned}$$

(1.3)

where we use the usual abbreviation $x^{+}{\mathop {=}\limits ^{\mathrm{{def}}}}\max \{ x,\,0\}$ and $x^{-}{\mathop {=}\limits ^{\mathrm{{def}}}}\max \{ -x,\,0\}$ for $x\in {\mathbb {R}}.$ Hence, $x = x^{+} - x^{-}$ and $|x| = x^{+} + x^{-}.$ These kinds of nonlinearities, often called “jumping nonlinearities, have a long tradition in Mathematical Modelling.

The parabolic partial differential equation (1.1) (PDE, for short) corresponds to the case when the nonlinearity $F(\,\cdot \,; S,t):\,M\mapsto F(M; S,t):\,{\mathbb {R}}\rightarrow {\mathbb {R}}$ on the right-hand side in Eq. (1.1) is taken with the mark-to-market value $M = \hat{V}(S,t).$ This case corresponds to a derivative contract $\hat{V}$ on an asset (stock) $S\in (0,\infty )$ between a  B and a  C that . The asset price S is not affected by a default of either B or C, and is assumed to follow the Markov process with the (time-dependent) generator (the Black–Scholes operator) ${\mathcal {A}}_t$ defined by

$$\begin{aligned}&({\mathcal {A}}_t V)(S,t){\mathop {=}\limits ^{\mathrm{{def}}}}\frac{1}{2}\, [\sigma (t)]^2 S^2\, \frac{\partial ^2 V}{\partial S^2} + [ q_S(t) - \gamma _S(t) ] S\, \frac{\partial V}{\partial S}\nonumber \\&\quad \hbox {for}\ V:\,(0,\infty )\times (0,T)\rightarrow {\mathbb {R}}:\,(S,t)\mapsto V(S,t) . \end{aligned}$$

(1.4)

As usual, we take the volatility, $\sigma ,$ to be a positive constant, $\sigma \in (0,\infty ).$ The value of $\gamma _S(t)$ reflects the rate of dividend income and the value of $q_S(t)$ is the net share position financing cost which depends on the risk-free rate r(t) and the repo-rate of S(t). “Typical” values for the terminal condition (1.2) are $h(S)\equiv V_T(S)$ where $V_T(S) = (S-K)^{+} = ( \mathrm {e}^X - K )^{+}$ for $X = \log ~S\in {\mathbb {R}}$ (in case of the European call option) and $V_T(S) = (S-K)^{-} = (K-S)^{+} = ( K - \mathrm {e}^X )^{+}$ (for the European put option).

The simple transformation $S\mapsto X = \log ~S:\,(0,\infty )\rightarrow {\mathbb {R}}$ of the asset (stock) price $S\in (0,\infty )$ into the logarithmic asset (stock) price $X = \log ~S\in {\mathbb {R}}$ has a “mathematical” justification in transforming the degenerate elliptic differential operator ${\mathcal {A}}_t$ in formula (1.4) into the regular differential operator ${\mathcal {A}}(\tau )$ in formula (2.3) below. Further connections between the asset price S and the logarithmic asset price $X = \log ~S$ can be found in the monograph by Fouque et al. [20, Sect. 1].

A frequently used alternative to our choice $M = \hat{V}(S,t)$ of the mark-to-market value M in the nonlinearity F(M; S, t) on the right-hand side in Eq. (1.1) is $M = V(S,t)$ where V denotes the same derivative between two parties that ; see e.g. Baustian et al. [6, Sect. 2] for numerical treatment. This risk-free value, V,  satisfies the classical (linear) Black–Scholes PDE (partial differential equation) with the prescribed terminal value $V(S,T) = h(S)$ for $S\in (0,\infty )$ at maturity time $t=T.$ Inserting this known value F(V(S, t); S, t) in Eq. (1.1) in place of $F(\hat{V}(S,t); S,t),$ we thus obtain an inhomogeneous linear equation for another (new) value of $\hat{V}(S,t).$ We refer to the works by Burgard and Kjaer [12, 13, Section 3], and [14] for details concerning modelling and to Arregui et al. [3] for numerical results. We warn the reader that Refs. [3, 12,13,14] use the convention $V = V^{+} + V^{-}$ with $V^{+}{\mathop {=}\limits ^{\mathrm{{def}}}}\max \{ V,\,0\}$ and $V^{-}{\mathop {=}\limits ^{\mathrm{{def}}}}\min \{ V,\,0\}$ (${}\le 0$) for $V\in {\mathbb {R}}$; nevertheless, we will stick with our notation $V = V^{+} - V^{-}$ with $V^{-}{\mathop {=}\limits ^{\mathrm{{def}}}}\max \{ -V,\,0\}$ (${}\ge 0$). We will not worry about this alternative any more and focus entirely on the nonlinear equation (1.1). Making use of Eq. (1.3), we arrive at the following equivalent form of Eq. (1.1), frequently used, cf. [13, Section 2, Eq. (1)]:

$$\begin{aligned}&\frac{\partial \hat{V}}{\partial t} + {\mathcal {A}}_t\hat{V} - r\hat{V} = {} - (1 - R_B)\lambda _B\, \hat{V}^{-} + (1 - R_C)\lambda _C\, \hat{V}^{+} + s_F\, \hat{V}^{+}\nonumber \\&\quad \hbox {for}\ (S,t)\in (0,\infty )\times (0,T) . \end{aligned}$$

(1.5)

This backward parabolic equation is supplemented by the terminal condition (1.2).

Models with are neither popular nor very frequent in Mathematical Finance. In the present article we treat a class of semilinear parabolic equations of type (1.5) with the standard linear diffusion operator $\frac{\partial }{\partial t} + {\mathcal {A}}_t$ and the nonlinear reaction function that is more general that the one on the right-hand side of (1.5) (only uniformly Lipschitz-continuous). As far as we know, this class was introduced in the work by Burgard and Kjaer [13, Section 3] and [14]. Another class of nonlinear models is based on a nonlinear Black–Scholes PDE with the quasilinear diffusion operator $\frac{\partial \hat{V}}{\partial t} + \frac{1}{2}\, \sigma ^2 S^2\, \frac{\partial ^2 V}{\partial S^2} + \cdots ,$ where the volatility $\sigma \equiv \sigma \genfrac(){}1{\partial ^2 V}{\partial S^2}$ depends on the second partial derivative, and with a “typical” linear reaction function (sometimes including also transaction costs). This class can be traced to Barles and Soner [5, Eq. (1.2), p. 372] with some additional analytic studies (on explicit solutions) performed in Bordag and Chmakova [9]. Some additional references to related numerical studies and simulations will be added in Sects. 4 and 5.

Last but not least, we would like to mention “modelling of incertitude in the environment” investigated in the works by Díaz et al. [16] and Díaz and Faghloumi [15, 17]. The stochastic formulation in these problems leads to degenerate obstacle problems closely related to parabolic problems with a free boundary that arise in Black–Scholes PDEs for American options. We believe that our current developments of Sattinger’s monotone methods [42] are applicable also to these kinds of Black–Scholes PDEs with a free boundary.

This article is organized as follows. We begin with a functional-analytic reformulation of the B–S equation (1.5) in the next section (Sect. 2). The terminal value problem (1.5), (1.2) will be transformed into an initial value Cauchy problem of parabolic type. This Cauchy problem is an initial value problem for the nonlinear (semilinear) B–S equation with a uniformly Lipschitz-continuous (nonlinear) reaction function, as well. In Sect. 3 we construct a of supersolutions of this B–S equation that converge as a monotone decreasing (i.e., nonincreasing) sequence to the solution from above; see our main result, Theorem 3.4. A closely related ramification of this monotone iteration scheme provides an increasing sequence of subsolutions of the B–S equation that converge to the solution from below; see Remark 3.5.

Numerical methods play an important role in Mathematical Finance. In Sect. 4 we discuss applications of two most common methods to Mathematical Finance, finite differences/elements and Monte Carlo. We discuss their advantages and problems when compared to each other. Finally, in Sect. 5 we derive an explicit formula for the solution of the inhomogeneous linear parabolic initial value problem for the B–S equation that serves for computing the monotone iteration scheme in Sect. 3. This formula is obtained by (with integrals over ${\mathbb {R}}^1$ and [0, T]) which makes it interesting for Monte Carlo computations. On the other hand, the solution of the inhomogeneous linear parabolic problem can be computed also by finite differences/elements.

2 Functional-analytic reformulation of the B–S equation

We wish to treat the terminal value problem (1.5) (or, equivalently, (1.1)) above, with the terminal condition (1.2), by standard analytic and numerical methods for semilinear parabolic initial value problems. To this end, we rewrite problem (1.5), (1.2) as the following general initial value problem for the unknown function $v:\,{\mathbb {R}}^1\times (0,T)\rightarrow {\mathbb {R}},$

$$\begin{aligned} \frac{\partial v}{\partial \tau } - {\mathcal {A}}(\tau ) v + r\, v&{} = \tilde{F}(v(x,\tau ); x, \tau )&\quad \hbox { for }\, (x,\tau )\in {\mathbb {R}}^1\times (0,T); \end{aligned}$$

(2.1)

$$\begin{aligned} v(x,0)&{}= v_0(x){\mathop {=}\limits ^{\mathrm{{def}}}}h(\mathrm {e}^x)&\quad \hbox { for }\, x\in {\mathbb {R}}^1 , \end{aligned}$$

(2.2)

where ${\mathcal {A}}(\tau )$ denotes the Black–Scholes operator defined by

$$\begin{aligned}&({\mathcal {A}}(\tau ) v)(x,\tau ){\mathop {=}\limits ^{\mathrm{{def}}}}({\mathcal {A}}_{T-\tau } v)(x,\tau ) \nonumber \\&\quad {} = \frac{1}{2}\, [\sigma (T-\tau )]^2 \,\frac{\partial ^2 v}{\partial x^2} + \left( q_S(T-\tau ) - \gamma _S(T-\tau ) - \frac{1}{2}\, [\sigma (T-\tau )]^2 \right) \frac{\partial v}{\partial x} \nonumber \\&\qquad \hbox {for }\; v:\,{\mathbb {R}}^1\times (0,T)\rightarrow {\mathbb {R}}:\,(x,\tau )\mapsto v(x,\tau ) , \end{aligned}$$

(2.3)

and the nonlinearity $\tilde{F}(\,\cdot \, ; x, \tau ):\,{\mathbb {R}}\rightarrow {\mathbb {R}}$ is given by

$$\begin{aligned} \tilde{F}(v; x, \tau )&{}{\mathop {=}\limits ^{\mathrm{{def}}}}{}- F(v; \mathrm {e}^x, T-\tau ) - (\lambda _B + \lambda _C)\, v \nonumber \\&{}= (1 - R_B)\lambda _B\, v^{-} - (1 - R_C)\lambda _C\, v^{+} - s_F\, v^{+} \nonumber \\&\quad \hbox { for }\, v\in {\mathbb {R}}\,\hbox { and }\, (x,\tau )\in {\mathbb {R}}^1\times (0,T) . \end{aligned}$$

(2.4)

Here, $\tau = T-t$ stands for the time to maturity and $x = \log ~S$ is the logarithmic asset (stock) price; we take $(x,\tau )\in {\mathbb {R}}^1\times (0,T).$ In the sequel we will never use the real time $t = T - \tau \in (0,T)$ any more, so we prefer to use the letter t in place of $\tau$ to denote the time to maturity, as it is usual in parabolic problems. According to this new notation, in Eq. (2.3) we replace the time-dependent coefficients $\sigma (T-\tau ),$ $q_S(T-\tau ),$ and $\gamma _S(T-\tau )$ by $\sigma (t),$ $q_S(t),$ and $\gamma _S(t),$ respectively, and thus forget about the original terminal value problem (1.5), (1.2):

$$\begin{aligned} ({\mathcal {A}}(t) v)(x,t)&{\mathop {=}\limits ^{\mathrm{{def}}}}\frac{\partial }{\partial x}\, \left[ \frac{1}{2}\, [\sigma (t)]^2 \,\frac{\partial v}{\partial x} + \left( q_S(t) - \gamma _S(t) - \frac{1}{2}\, [\sigma (t)]^2 \right) v(x,t) \right] \nonumber \\&\quad \hbox {for }\; v:\,{\mathbb {R}}^1\times (0,T)\rightarrow {\mathbb {R}}:\,(x,t)\mapsto v(x,t) . \end{aligned}$$

(2.5)

Next, in order to make the initial value problem (2.1), (2.2) compatible with the described in the article by Sattinger [42], we rewrite this problem as follows:

$$\begin{aligned} \frac{\partial v}{\partial t} - {\mathcal {A}}(t) v + (r + L_{\tilde{F}})\, v = \tilde{F}(v(x,t); x,t) + L_{\tilde{F}}\, v \quad \hbox {for }\, (x,t)\in {\mathbb {R}}^1\times (0,T); \end{aligned}$$

(2.6)

with the initial condition $v(x,0) = v_0(x){\mathop {=}\limits ^{\mathrm{{def}}}}h(\mathrm {e}^x)$ for $x\in {\mathbb {R}}^1$ in Eq. (2.2), where the constant $L_{\tilde{F}}\in {\mathbb {R}}_+ = [0,\infty )$ is defined by

$$\begin{aligned} L_{\tilde{F}} = \max \{ (1 - R_B)\lambda _B ,\, (1 - R_C)\lambda _C + s_F \} . \end{aligned}$$

(2.7)

According to [13], $s_F\equiv r_F - r,$ $\lambda _B\equiv r_B - r,$ and $\lambda _C\equiv r_C - r$ are some nonnegative constants and $R_B, R_C\in [0,1]$ are the recovery rates on the derivative positions of parties B and C,  respectively. As a consequence, the function $G(\,\cdot \, ; x,t):\,v\mapsto G(v; x,t):\,{\mathbb {R}}\rightarrow {\mathbb {R}},$ defined by

$$\begin{aligned} G(v; x,t)&= \tilde{F}(v; x,t) + L_{\tilde{F}}\, v \nonumber \\&= {}- \left[ L_{\tilde{F}} - (1 - R_B)\lambda _B \right] v^{-} + \left[ L_{\tilde{F}} - (1 - R_C)\lambda _C - s_F \right] v^{+} \nonumber \\&\quad \hbox {for }v\in {\mathbb {R}}\hbox { and }\, (x,t)\in {\mathbb {R}}^1\times (0,T) , \end{aligned}$$

(2.8)

is monotone increasing (i.e., nondecreasing) on ${\mathbb {R}}.$ Notice that both functions, $v\mapsto {}- v^{-}$ and $v\mapsto v^{+},$ are nondecreasing on ${\mathbb {R}}.$ Indeed, we have also

$$\begin{aligned}&\frac{\partial G}{\partial v} (v; x,t) = \frac{ \partial \tilde{F} }{\partial v} + L_{\tilde{F}} = \left\{ \quad \begin{array}{ll} L_{\tilde{F}} - (1 - R_B)\lambda _B &{}\quad \hbox {if } v < 0,\\ L_{\tilde{F}} - \left( (1-R_C)\lambda _C + s_F \right) &{}\quad \hbox {if } v > 0 ; \end{array} \right. \\&\hbox {with } 0\le \frac{\partial G}{\partial v} (v; x,t) \le L_{\tilde{F}} \quad \hbox {for all } v\in {\mathbb {R}}^1\setminus \{ 0\} \hbox { and } (x,t)\in {\mathbb {R}}\times (0,T) . \end{aligned}$$

In addition, the left-hand side of Eq. (2.6) clearly satisfies the weak maximum principle.

An alternative to Sattinger’s article [42] on monotone methods for “Nonlinear Elliptic and Parabolic Problems” is offered in a book form in the monograph Pao [39, Chapt. 1], §1.5 (pp. 20–26) and §1.7 (pp. 31–36), and further in Chapters 2 and 3. Numerous well-known details are included here, such as the relation between the one-sided Lipschitz condition expressed in the inequality on the left-hand side,

$$\begin{aligned} {}- L_{\tilde{F}} \le \frac{ \partial \tilde{F} }{\partial v} \le L_{\tilde{F}} \quad \hbox { for all }\, v\in {\mathbb {R}}^1\setminus \{ 0\} \,\hbox { and }\, (x,t)\in {\mathbb {R}}\times (0,T) , \end{aligned}$$

(2.9)

and the monotonicity hypothesis (G3) below.

We now specify our hypotheses on the general nonlinearity $G:\,{\mathbb {R}}\times {\mathbb {R}}^1\times (0,T)$ on the right-hand side of Eq. (2.6) treated in our present work. We assume that G satisfies the following hypotheses:

Hypotheses

(G1):

For each fixed $v\in {\mathbb {R}},$ the function $G(v; \,\cdot , \,\cdot \,):\,(x,t)\mapsto G(v;x,t):\,{\mathbb {R}}^1\times (0,T)\rightarrow {\mathbb {R}}$ is Lebesgue-measurable.

(G2):

For almost every pair $(x,t)\in {\mathbb {R}}^1\times (0,T),$ the function $G(\,\cdot \, ; x,t):\,v\mapsto G(v; x,t):\,{\mathbb {R}}\rightarrow {\mathbb {R}}$ is uniformly Lipschitz-continuous with a Lipschitz constant $L_G\in {\mathbb {R}}_+,$ that is, we have

$$\begin{aligned}&| G(v_1;x,t) - G(v_2;x,t) | \le L_G\, |v_1 - v_2| \nonumber \\&\quad \hbox { for all }\, v_1, v_2\in {\mathbb {R}} \,\hbox { and for almost all }\, (x,t)\in {\mathbb {R}}^1\times (0,T) . \end{aligned}$$

(2.10)

(G3):

For almost every pair $(x,t)\in {\mathbb {R}}^1\times (0,T),$ the function $G(\,\cdot \, ; x,t):\,v\mapsto G(v; x,t):\,{\mathbb {R}}\rightarrow {\mathbb {R}}$ is monotone increasing, that is, $v_1\le v_2$ in ${\mathbb {R}}$ implies $G(v_1;x,t)\le G(v_2;x,t).$

(G4):

There is a constant $C_0\in {\mathbb {R}}_+$ such that, at almost every time $t\in (0,T),$ the function $G( 0; \,\cdot \, ,t):\,x\mapsto G(0; x,t):\,{\mathbb {R}}^1\rightarrow {\mathbb {R}}$ satisfies the exponential growth restriction

$$\begin{aligned} | G(0;x,t) |\le C_0\cdot \exp (|x|)\quad \bigl ( {}\le C_0 (\mathrm {e}^x + \mathrm {e}^{-x}) \bigr ) \quad \hbox { for almost all }\, x\in {\mathbb {R}}^1 . \end{aligned}$$

(2.11)

(G5):

There are constants $C_1\in {\mathbb {R}}_+$ and $\vartheta _G\in (0,1)$ such that, for every $v\in {\mathbb {R}}$ and almost every $x\in {\mathbb {R}}^1,$ the function $G(v;x, \,\cdot \, ):\,t\mapsto G(v;x,t):\,(0,T)\rightarrow {\mathbb {R}}$ is Hölder-continuous with the Hölder exponent $\vartheta _G$ ($0< \vartheta _G < 1$) in the following sense,

$$\begin{aligned}&| G(v;x,t_1) - G(v;x,t_2) | \le C_1\, |v|\cdot |t_1 - t_2|^{\vartheta _G} \nonumber \\&\quad \hbox { for all }\, t_1, t_2\in (0,T) \,\hbox { and for almost all }\, (v,x)\in {\mathbb {R}}\times {\mathbb {R}}^1 . \end{aligned}$$

(2.12)

From now on, let us consider the following generalization of the initial value problem (2.6), (2.2):

$$\begin{aligned} \frac{\partial v}{\partial t} - {\mathcal {A}}(t) v + r_G\, v = G(v(x,t); x,t) \quad \hbox {for } (x,t)\in {\mathbb {R}}^1\times (0,T) ; \end{aligned}$$

(2.13)

with the initial condition $v(x,0) = v_0(x){\mathop {=}\limits ^{\mathrm{{def}}}}h(\mathrm {e}^x)$ for $x\in {\mathbb {R}}^1$ in Eq. (2.2), where the constant $r + L_{\tilde{F}}$ in Eq. (2.6) has been replaced by the new constant $r_G\in {\mathbb {R}}_+,$ owing to our monotonicity hypothesis (G3). Concerning hypotheses on the time-dependent coefficients that appear in the Black–Scholes operator ${\mathcal {A}}(t)$ defined in Eq. (2.5) (recall that $\tau = t$), we assume the following Hölder continuity:

Hypotheses

(BS1):

$\sigma :\,[0,T]\rightarrow (0,\infty )$ is a positive, Hölder-continuous function satisfying

$$\begin{aligned} | \sigma (t_1) - \sigma (t_2) | \le C_{\sigma }\, |t_1 - t_2|^{ \vartheta _{\sigma } } \quad \hbox { for all }\, t_1, t_2\in [0,T] , \end{aligned}$$

(2.14)

where $C_{\sigma }\in {\mathbb {R}}_+$ and $\vartheta _{\sigma }\in (0,1)$ are some constants independent from time $t\in [0,T].$

(BS2):

$q_S, \gamma _S:\,[0,T]\rightarrow {\mathbb {R}}$ is a pair of Hölder-continuous function satisfying

$$\begin{aligned} | q_S(t_1) - q_S(t_2) |&{} \le C_q\, |t_1 - t_2|^{ \vartheta _q } \quad \hbox { and }\quad \end{aligned}$$

(2.15)

$$\begin{aligned} | \gamma _S(t_1) - \gamma _S(t_2) |&{} \le C_{\gamma }\, |t_1 - t_2|^{ \vartheta _{\gamma } } \quad \hbox { for all }\, t_1, t_2\in [0,T] , \end{aligned}$$

(2.16)

where $C_q, C_{\gamma }\in {\mathbb {R}}_+$ and $\vartheta _q, \vartheta _{\gamma }\in (0,1)$ are some constants (independent from $t\in [0,T]$).

Remark 2.1

(Hölder exponents) In Hypotheses (G5), (BS1), and (BS2) we may and will replace the Hölder exponents $\vartheta _G,$ $\vartheta _{\sigma },$ $\vartheta _q,$ and $\vartheta _{\gamma }$ by their minimum $\vartheta _0,$ $\vartheta _0 = \min \{ \vartheta _G ,\, \vartheta _{\sigma } ,\, \vartheta _q ,\, \vartheta _{\gamma } \} ,\quad \vartheta _0\in (0,1) .$ $\square$

Clearly, from Hypothesis (BS1) we derive $\sigma (t)\ge \sigma _0 = \min _{t\in [0,T]} \sigma (t) > 0$ for all $t\in [0,T].$ This fact, combined with (BS2), guarantees the uniform ellipticity of the Black–Scholes operator ${\mathcal {A}}(t)$ independently from $t\in [0,T].$

Remark 2.2

(Risk-free interest rate) One may also suggest to replace the multiplicative constant $r_G\in {\mathbb {R}}$ on the left-hand side of Eq. (2.13) by the time-dependent risk-free interest rate $r:\,[0,T]\rightarrow {\mathbb {R}}$ satisfying a Hölder continuity condition analogous to those in Eqs. (2.15) and (2.16). However, this change would not make Eq. (2.13) more general in that it could be reduced to the present form (2.13) with the term $r_G\, v$ as follows:

First, define $r_G\in {\mathbb {R}}$ by $r_G = \max _{t\in [0,T]} r(t)$; then replace the function G(v; x, t) on the right-hand side of Eq. (2.13) by the sum $G_r(v;x,t) = G(v;x,t) + [r_G - r(t)]\, v$ for $(v;x,t)\in {\mathbb {R}}\times {\mathbb {R}}^1\times (0,T).$ Clearly, thanks to $r_G - r(t)\ge 0$ for every $t\in [0,T],$ the function $G_r(\,\cdot \,; \,\cdot , \,\cdot \,):\,(v;x,t)\mapsto G_r(v;x,t):\,{\mathbb {R}}^1\times (0,T)\rightarrow {\mathbb {R}}$ satisfies all Hypotheses (G1)–(G5) imposed on the function G. We conclude that the interest rate difference, $r_G - r(t),$ can be included in the reaction function G. We thus keep Eq. (2.13) in the present form with $r_G\in {\mathbb {R}}$ being a given constant. $\square$

Our last hypothesis in problem (2.13), (2.2) restricts the growth of the initial condition $v(x,0) = v_0(x){\mathop {=}\limits ^{\mathrm{{def}}}}h(\mathrm {e}^x)$ for $x\in {\mathbb {R}}^1$ in Eq. (2.2) as follows:

Hypothesis

(v$_{\mathbf {0}}$):

The function $v_0:\,{\mathbb {R}}\rightarrow {\mathbb {R}}$ is Lebesgue-measurable and there is a constant $C_h\in {\mathbb {R}}_+$ such that, for almost all $x\in {\mathbb {R}}^1,$ we have

$$\begin{aligned} |v_0(x)| = |h(\mathrm {e}^x)|\le C_h\cdot \exp (|x|) \quad \bigl ( {}\le C_h (\mathrm {e}^x + \mathrm {e}^{-x}) \bigr ) . \end{aligned}$$

(2.17)

As we have already indicated in our hypothesis (G4) on the exponential growth restriction of G,  we are going to look for (strong, weak or mild) solutions $v:\,{\mathbb {R}}^1\times (0,T)\rightarrow {\mathbb {R}}$ to the initial value problem (2.13), (2.2) satisfying an analogous exponential growth restriction of type $v( \,\cdot \, ,t)\in H_{\mathbb {C}}$ at every time $t\in (0,T),$ where $H_{\mathbb {C}} = L^2({\mathbb {R}};{\mathfrak {w}})$ denotes the complex Hilbert space of all complex-valued Lebesgue-measurable functions $f:\,{\mathbb {R}}\rightarrow \mathbb {C}$ with the finite norm

$$\begin{aligned} \Vert f\Vert _H{\mathop {=}\limits ^{\mathrm{{def}}}}\left( \int _{{\mathbb {R}}} |f(x)|^2\, {\mathfrak {w}}(x) \,{\mathrm {d}}x \right) ^{1/2} < \infty , \end{aligned}$$

where ${\mathfrak {w}}(x){\mathop {=}\limits ^{\mathrm{{def}}}}\mathrm {e}^{-\mu |x|}$ is a weight function with some constant $\mu \in (2,\infty ).$ This norm is induced by the inner product

$$\begin{aligned} (f,g)_H\equiv (f,g)_{ L^2({\mathbb {R}};{\mathfrak {w}}) } {\mathop {=}\limits ^{\mathrm{{def}}}}\int _{{\mathbb {R}}} f\, \bar{g}\cdot {\mathfrak {w}}(x) \,{\mathrm {d}}x \quad \hbox { for }\, f,g\in H_{\mathbb {C}} . \end{aligned}$$

As usual, the symbol $\bar{z}$ denotes the complex conjugate of a complex number $z\in \mathbb {C}$ where ${\mathbb {C}} = {\mathbb {R}}+ \mathrm {i}{\mathbb {R}}$ is the complex plane. We consider the complex Hilbert space $H_{\mathbb {C}}$ only for better understanding of our applications using holomorphic semigroups in $H_{\mathbb {C}}$ generated by the (unbounded) Black–Scholes operator ${\mathcal {A}}(t):\,H_{\mathbb {C}}\rightarrow H_{\mathbb {C}}$ in Eq. (2.13) above. Our solutions v(x, t) to the initial value problem (2.13), (2.2) will be always real-valued, i.e., $v( \,\cdot \, ,t)\in H$ at every time $t\in (0,T),$ where H denotes the closed real vector subspace of all real-valued functions $f:\,{\mathbb {R}}\rightarrow {\mathbb {R}}$ from $H_{\mathbb {C}}.$ The domain of the differential operator ${\mathcal {A}}(t),$ denoted by ${\mathcal {D}}({\mathcal {A}}(t)),$ is a complex vector subspace of $H_{\mathbb {C}}$ which is independent from time $t\in [0,T],$ i.e., ${\mathcal {D}}({\mathcal {A}}(t)) \equiv D_{\mathbb {C}}\subset H_{\mathbb {C}}$ for every $t\in [0,T].$ The vector space $D_{\mathbb {C}}$ consists of all functions $f\in H_{\mathbb {C}}$ whose weak (distributional) derivatives $f'= \frac{{\mathrm {d}}f}{{\mathrm {d}}x}$ and $f''= \frac{{\mathrm {d}}^2 f}{{\mathrm {d}}x^2}$ belong to $H_{\mathbb {C}},$ as well. We set $D = D_{\mathbb {C}}\cap H$ to denote the closed real vector subspace of all real-valued functions $f:\,{\mathbb {R}}\rightarrow {\mathbb {R}}$ from $D_{\mathbb {C}}.$ The vector space $D_{\mathbb {C}}$ becomes a Banach space under the norm

$$\begin{aligned} \Vert f\Vert _D{\mathop {=}\limits ^{\mathrm{{def}}}}\Vert f\Vert _H + \Vert f''\Vert _H \quad \hbox { for }\, f\in D_{\mathbb {C}} . \end{aligned}$$

This norm is equivalent with the stronger norm

$$\begin{aligned} {\left| \left| \left| f \right| \right| \right| }_D{\mathop {=}\limits ^{\mathrm{{def}}}}\Vert f\Vert _H + \Vert f'\Vert _H + \Vert f''\Vert _H \quad \hbox { for }\, f\in D_{\mathbb {C}} , \end{aligned}$$

by a simple interpolation inequality. We refer to the monograph by Engel and Nagel [18] for details concerning holomorphic semigroups and their (infinitesimal) generators, especially to [18, Chapt. II, Sect. 4a,  pp. 96–109].

We denote by $H_{\mathbb {C}}^1$ the complex interpolation space between $D_{\mathbb {C}}$ and $H_{\mathbb {C}}$ that consist of all functions $f:\,{\mathbb {R}}\rightarrow {\mathbb {R}}$ from $H_{\mathbb {C}}$ such that both $f, f'\in H_{\mathbb {C}}.$ $H_{\mathbb {C}}^1$ is a vector space which becomes a Banach space under the norm

$$\begin{aligned} \Vert f\Vert _{H^1}{\mathop {=}\limits ^{\mathrm{{def}}}}\Vert f\Vert _H + \Vert f'\Vert _H \quad \hbox { for }\, f\in H_{\mathbb {C}}^1 . \end{aligned}$$

Hence, we have the continuous imbeddings $D_{\mathbb {C}}\hookrightarrow H_{\mathbb {C}}^1\hookrightarrow H_{\mathbb {C}} .$ Moreover, given any fixed $t\in [0,T],$ $H_{\mathbb {C}}^1$ is the domain of the sesquilinear form

$$\begin{aligned}&{\mathcal {Q}}(t):\,H_{\mathbb {C}}^1\times H_{\mathbb {C}}^1\rightarrow \mathbb {C}:\,(f,g) \;\longmapsto \; {\mathcal {Q}}(t)(f,g){\mathop {=}\limits ^{\mathrm{{def}}}}{}- \int _{{\mathbb {R}}} [{\mathcal {A}}(t) f](x)\, \overline{g(x)} \cdot {\mathfrak {w}}(x) \,{\mathrm {d}}x \nonumber \\&\quad = {}- \int _{-\infty }^{+\infty } \left\{ \genfrac{}{}{}1{1}{2}\, [\sigma (t)]^2\, f''(x) + \left( q_S(t) - \gamma _S(t) - \genfrac{}{}{}1{1}{2}\, [\sigma (t)]^2 \right) f'(x) \right\} \overline{g(x)}\cdot {\mathfrak {w}}(x) \,{\mathrm {d}}x \nonumber \\&\quad = \int _{-\infty }^{+\infty } \left\{ \genfrac{}{}{}1{1}{2}\, [\sigma (t)]^2\, f'(x)\, \overline{g'(x)} - \genfrac{}{}{}1{\mu }{2}\, [\sigma (t)]^2\, \mathop {\mathrm {sign}}(x)\, f'(x)\, \overline{g(x)} \right. \nonumber \\&\qquad {} - \left. \left( q_S(t) - \gamma _S(t) - \genfrac{}{}{}1{1}{2}\, [\sigma (t)]^2 \right) f'(x)\, \overline{g(x)} \right\} \cdot {\mathfrak {w}}(x) \,{\mathrm {d}}x \nonumber \\&\quad {} = \genfrac{}{}{}1{1}{2}\, [\sigma (t)]^2 \cdot \int _{-\infty }^{+\infty } f'(x)\, \overline{g'(x)}\cdot {\mathfrak {w}}(x) \,{\mathrm {d}}x \nonumber \\&\qquad {} + \genfrac{}{}{}1{\mu }{2}\, [\sigma (t)]^2 \cdot \left( \int _{-\infty }^0 - \int _0^{+\infty } \right) f'(x)\, \overline{g(x)}\cdot {\mathfrak {w}}(x) \,{\mathrm {d}}x \nonumber \\&\qquad {} - \left( q_S(t) - \gamma _S(t) - \genfrac{}{}{}1{1}{2}\, [\sigma (t)]^2 \right) \int _{-\infty }^{+\infty } f'(x)\, \overline{g(x)}\cdot {\mathfrak {w}}(x) \,{\mathrm {d}}x \nonumber \\&\qquad \quad \hbox {defined first only for }\, f,g\in D_{\mathbb {C}}. \end{aligned}$$

(2.18)

The continuous extension of ${\mathcal {Q}}(t)(f,g)$ to all $f,g\in H_{\mathbb {C}}^1$ is immediate, thanks to $D_{\mathbb {C}}$ being a dense vector subspace of $H_{\mathbb {C}}^1.$ A few simple applications of the Cauchy–Schwarz inequality show that the (non-symmetric) sesquilinear form ${\mathcal {Q}}(t)$ on $H_{\mathbb {C}}^1$ is coercive. Indeed, with a help from Hypothesis (BS1) we have $\sigma (t)\ge \sigma _0 = \min _{t\in [0,T]} \sigma (t) > 0$ for all $t\in [0,T].$ Consequently, if the constant $\lambda _0\in (0,\infty )$ below is chosen sufficiently large, then we get the following more precise quantification of coercivity at every time $t\in [0,T]$:

$$\begin{aligned} {\mathcal {Q}}(t)(f,f) + \lambda \, (f,f)_H \ge \frac{\sigma _0}{4}\, \Vert f\Vert _{H^1}^2 + \Vert f\Vert _H^2 \quad \hbox { for every }\, \lambda \ge \lambda _0 , \end{aligned}$$

(2.19)

thanks to Hypotheses (BS1) and (BS2). We note that the constant $\lambda _0$ depends neither on time $t\in [0,T]$ nor on the number $\lambda \ge \lambda _0.$

Let $I\equiv I_H$ denote the identity mapping on $H_{\mathbb {C}}.$ Given any real number $\lambda \ge \lambda _0,$ from ineq. (2.19) we infer that the linear operator

$$\begin{aligned} {}- {\mathcal {A}}_{\lambda }(t) {\mathop {=}\limits ^{\mathrm{{def}}}}{}- {\mathcal {A}}(t) + \lambda I :\,D_{\mathbb {C}}\subset H_{\mathbb {C}}\rightarrow H_{\mathbb {C}} \end{aligned}$$

is an isomorphism of the Banach space $D_{\mathbb {C}}$ onto another Banach space $H_{\mathbb {C}},$ both, algebraically and topologically. Now we can apply the well-known results for abstract linear initial value problems of parabolic type, e.g., from Evans [19], Chapt. 7, §1.1, p. 352, or Lions [35], Chapt. IV, §1, p. 44, or [36], Chapt. III, eq. (1.11), p. 102, to conclude that the inhomogeneous linear parabolic initial value problem

$$\begin{aligned} \frac{\partial v}{\partial t} - {\mathcal {A}}(t) v + r_G\, v = f(x,t) \quad \hbox {for}\ (x,t)\in {\mathbb {R}}^1\times (0,T) ; \end{aligned}$$

(2.20)

with the initial condition $v(x,0) = v_0(x){\mathop {=}\limits ^{\mathrm{{def}}}}h(\mathrm {e}^x)$ for $x\in {\mathbb {R}}^1$ in Eq. (2.2), possesses a unique weak solution $v:\,[0,T]\rightarrow H,$ whenever the initial value $v_0\in H$ is given, such that $v_0:\,{\mathbb {R}}^1\rightarrow {\mathbb {R}}$ obeys ineq. (2.17) in Hypothesis (v$_{\mathbf {0}}$). The weak solution v is continuous as an H-valued function of time $t\in [0,T],$ that is, $v\in C([0,T]\rightarrow H).$ This is the “linear part” of the semilinear problem (2.13), (2.2) with a prescribed inhomogeneity $f:\,[0,T]\rightarrow H$ that is assumed to be strongly Lebesgue-measurable and (essentially) bounded on (0, T),  i.e., $f\in L^{\infty }((0,T)\rightarrow H).$ In our case, $f\in L^{\infty }((0,T)\rightarrow H)$ follows from our stronger hypothesis below:

Hypothesis

(f1):

$f:\,[0,T]\rightarrow H\subset H_{\mathbb {C}} = L^2({\mathbb {R}};{\mathfrak {w}})$ is a continuous function, i.e., $f\in C([0,T]\rightarrow H),$ where $f:\,{\mathbb {R}}^1\times [0,T]\rightarrow {\mathbb {R}}:\,(x,t)\mapsto f(x,t)$ satisfies $f(t)\equiv f( \,\cdot \, ,t)\in H$ for every $t\in [0,T].$

Although in the following sections we work only with weak solutions $v:\,[0,T]\rightarrow H,$ $v(0) = v_0\in H,$ to the “linear part” (2.20) of the semilinear equation (2.13), we would like to remark that the unique weak solution $v:\,[0,T]\rightarrow H$ to the linear initial value problem (2.20), (2.2) becomes a (unique) classical solution if $f:\,[0,T]\rightarrow H$ satisfies the following stronger, Hölder-continuity hypothesis with the Hölder exponent $\vartheta _f\in (0,1)$:

Hypothesis

(f1$^{\prime }$):

$f:\,[0,T]\rightarrow H\subset H_{\mathbb {C}} = L^2({\mathbb {R}};{\mathfrak {w}})$ is a $\vartheta _f$-Hölder-continuous function, i.e., there are constants $\vartheta _f\in (0,1)$ and $C_f\in {\mathbb {R}}_+$ such that

$$\begin{aligned} \Vert f( \,\cdot \, ,t_1) - f( \,\cdot \, ,t_2)\Vert _H \le C_f\cdot |t_1 - t_2|^{\vartheta _f} \quad \hbox {for all}\ t_1, t_2\in [0,T] . \end{aligned}$$

(2.21)

This is the case if the inhomogeneity $f:\,{\mathbb {R}}^1\times [0,T]\rightarrow {\mathbb {R}}$ satisfies $f( \,\cdot \, ,t)\in H$ for each $t\in [0,T]$ and there are some constants $\tilde{C}_f\in {\mathbb {R}}_+$ and $\kappa \in {\mathbb {R}},$ with $1\le \kappa < \mu /2,$ such that

$$\begin{aligned}&|f(x,t_1) - f(x,t_2)| \le \tilde{C}_f\, \mathrm {e}^{\kappa |x|}\cdot |t_1 - t_2|^{\vartheta _f}\nonumber \\&\quad \hbox { for a.e. (almost every) }x\in {\mathbb {R}}^1 \hbox { and for all }\, t_1, t_2\in [0,T] . \end{aligned}$$

(2.22)

Recall that $\mu$ ($\mu > 2$) is the constant in the weight function ${\mathfrak {w}}(x){\mathop {=}\limits ^{\mathrm{{def}}}}\mathrm {e}^{-\mu |x|}$ in the Hilbert space $H_{\mathbb {C}} = L^2({\mathbb {R}};{\mathfrak {w}}).$ We remark that for the inhomogeneous linear equation (2.20), the nonlinearity G on the right-hand side of Eq. (2.13) becomes $G(v;x,t)\equiv f(x,t)$; hence, we have $\vartheta _G = \vartheta _f$ in Hypothesis (G5). Let us recall that, by Remark 2.1, we have replaced the Hölder exponents $\vartheta _G,$ $\vartheta _{\sigma },$ $\vartheta _q,$ and $\vartheta _{\gamma }$ by their minimum $\vartheta _0$; hence, we may include also the value of $\vartheta _f$ in that minimum:

$$\begin{aligned} \vartheta _0 = \min \{ \vartheta _G ,\, \vartheta _{\sigma } ,\, \vartheta _q ,\, \vartheta _{\gamma } ,\, \vartheta _f \} ,\quad \vartheta _0\in (0,1) . \end{aligned}$$

(2.23)

Indeed, according to the existence, uniqueness, and regularity results for problem (2.20), (2.2) in Pazy [40, Chapt. 5, §5.7], Theorem 7.1 on p. 168, if (2.21) holds, then the unique weak solution $v:\,[0,T]\rightarrow H$ to the linear initial value problem (2.20), (2.2) described above happens to be a unique classical solution which, among other properties, is continuous as a function $v:\,[0,T]\rightarrow H,$ i.e., $v\in C([0,T]\rightarrow H),$ continuously differentiable on the time interval (0, T],  $v(t)\equiv v( \,\cdot \, ,t)\in D$ for every $t\in (0,T],$ i.e., ${\mathcal {A}}(t) v(t)\in H$ for $t\in (0,T],$ and v satisfies the abstract differential equation

$$\begin{aligned} \frac{\partial v}{\partial t} - {\mathcal {A}}(t) v + r_G\, v&{} = f(t)&\quad \hbox { in }H \hbox { for }\, t\in (0,T) ; \end{aligned}$$

(2.24)

$$\begin{aligned} \quad \hbox {with}\quad v(0)&{}= v_0\in H&\quad \hbox { (the initial condition).} \end{aligned}$$

(2.25)

Although in our parabolic evolutionary problems we make use of the monograph by Pazy [40], we would like to point to a closely related (and much newer) monograph by Vrabie [45] for interesting alternatives to results in [40].

In typical applications of our results to Mathematical Finance, in particular, to counterparty risk models treated in our present work, both, the coefficients in the Black–Scholes operator ${\mathcal {A}}(\tau )$ defined in Eq. (2.3) and the nonlinearity $\tilde{F}(\,\cdot \, ; x, \tau ):\,{\mathbb {R}}\rightarrow {\mathbb {R}}$ defined in Eq. (2.4), are assumed to be continuous, or even Hölder-continuous, in time $\tau \in [0,T].$ In contrast, less restrictive time regularity hypotheses are needed in Stochastic Control Theory, such as piecewise continuity in time $\tau \in [0,T].$ Because our present work does not treat problems in “Stochastic Control Theory”, we refer the interested reader to the monograph by Bensoussan and Lions [8] for this interesting topic and for methods how to relax our time regularity hypotheses.

3 Monotone methods for the nonlinear B–S equation

We make use of the inhomogeneous linear problem (2.20), (2.2) in order to describe an iterative scheme for approximating the unique weak solution $v:\,[0,T]\rightarrow H,$ $v(0) = v_0\in H,$ to the semilinear problem (2.13), (2.2).

3.1 Preliminary comparison results for parabolic problems

First, the so-called weak maximum principle for a classical solution v of the inhomogeneous linear problem (2.20), (2.2) is established e.g. in Friedman [22, Chapt. 2, Sect. 4, Theorem 9, p. 43]. A standard approximation procedure of a weak solution v by a sequence of classical solutions yields the corresponding weak maximum principle also for the weak solution v. More precisely, if the inequalities $v(x,0) = v_0(x)\ge 0$ and $f(x,t)\ge 0$ are valid for almost all $(x,t)\in {\mathbb {R}}^1\times (0,T),$ then also $v(x,t)\ge 0$ holds for almost all $(x,t)\in {\mathbb {R}}^1\times (0,T).$

Second, let $v:\,[0,T]\rightarrow H$ be a classical solution of the inhomogeneous linear problem

$$\begin{aligned} \frac{\partial v}{\partial t} - {\mathcal {A}}(t) v + r_G\, v&{} = g(t)&\quad \hbox { in }H \hbox { for }\, t\in (0,T) ; \end{aligned}$$

(3.1)

$$\begin{aligned} \quad \hbox { with }\quad v(0)&{}= v_g\in H&\quad \hbox { (the initial condition),} \end{aligned}$$

(3.2)

where $g:\,[0,T]\rightarrow H$ is a function continuous and bounded in (0, T),  i.e., $g\in C((0,T)\rightarrow H)$ with $\Vert g\Vert _{ L^{\infty }((0,T)\rightarrow H) }{\mathop {=}\limits ^{\mathrm{{def}}}}\sup _{t\in (0,T)} \Vert g(t)\Vert _H < \infty .$ We say that v is a supersolution to the inhomogeneous linear problem (2.20), (2.2), if the inequalities $v_g(x)\ge v_0(x)$ and $g(x,t)\ge f(x,t)$ hold for almost all $(x,t)\in {\mathbb {R}}^1\times (0,T).$ Analogously, a subsolution v to the inhomogeneous linear problem (2.20), (2.2) is a classical solution v of problem (3.1), (3.2) for which the inequalities $v_g(x)\le v_0(x)$ and $g(x,t)\le f(x,t)$ hold for almost all $(x,t)\in {\mathbb {R}}^1\times (0,T).$ More generally, if $v:\,(x,t)\in {\mathbb {R}}^1\times (0,T)\rightarrow {\mathbb {R}}$ is a weak solution to the inhomogeneous linear problem (2.20), (2.2), then the notions of super- and subsolution to the inhomogeneous linear problem (2.20), (2.2) are defined by means of an approximation procedure by a sequence of classical solutions again. A rigorous functional-analytic way to define super- and subsolution v to problem (2.20), (2.2) is to require that $v:\,(x,t)\in {\mathbb {R}}^1\times (0,T)\rightarrow {\mathbb {R}}$ have all regularity properties of a weak solution to problem (2.20), (2.2) stated in Evans [19], Chapt. 7, §1.1, p. 352, or Lions [35], Chapt. IV, §1, p. 44, or [36], Chapt. III, eq. (1.11), p. 102, and, in addition to these regularity properties, the following inequality is valid in the sense of distributions on ${\mathbb {R}}^1\times (0,T)$:

$$\begin{aligned} \frac{\partial v}{\partial t} - {\mathcal {A}}(t) v + r_G\, v&{} \ge f(t) \quad ({}\le f(t))&\quad \hbox {in }H \hbox { for } t\in (0,T) ; \end{aligned}$$

(3.3)

$$\begin{aligned} \quad \hbox {with}\quad v(0)&{}\ge v_0\in H \quad ({}\le v_0\in H)&\quad \hbox { (the initial condition).} \end{aligned}$$

(3.4)

Here, the inequalities with “$\,\ge \,$” (“$\,\le \,$”, respectively) specify a (weak) supersolution (a (weak) subsolution). The reader is referred to Friedman [21, Chapt. 3, Sect. 3, Theorem 8, p. 51] for details about positive (or nonnegative) distributions. Clearly, any function $v:\,(x,t)\in {\mathbb {R}}^1\times (0,T)\rightarrow {\mathbb {R}}$ which is simultaneously a (weak) supersolution and a (weak) subsolution) of the inhomogeneous linear problem (2.20), (2.2) is a (weak) solution to this problem. Combining these definitions of super- and subsolution, denoted by $\overline{v}, \underline{v}:\,{\mathbb {R}}^1\times (0,T)\rightarrow {\mathbb {R}},$ respectively, having the initial values satisfying $\underline{v}(0)\le v_0\le \overline{v}(0)$ a.e. in ${\mathbb {R}}^1,$ with the weak maximum principle for the difference $w = \overline{v} - \underline{v},$ we obtain the following auxiliary weak comparison result.

Lemma 3.1

(Weak comparison). Assume that $\overline{v}, \underline{v}:\,{\mathbb {R}}^1\times (0,T)\rightarrow {\mathbb {R}},$ respectively,  is a pair of (weak) super- and subsolutions of problem (2.20),  (2.2) satisfying $\underline{v}(0)\le \overline{v}(0)$ a.e. in ${\mathbb {R}}^1.$ Then,  at every time $t\in [0,T),$ we have $\underline{v}(x,t)\le \overline{v}(x,t)$ for a.e. $x\in {\mathbb {R}}^1.$

Observe that we have left the initial value $v_0\in H$ out of this lemma since we use it usually with either $\underline{v}(0) = v_0$ or $\overline{v}(0) = v_0$ as the initial condition attached to the differential equation

$$\begin{aligned} \frac{ \partial \underline{v} }{\partial t} - {\mathcal {A}}(t) \underline{v} + r_G\, \underline{v} = \underline{f}(x,t) \quad \hbox {for } (x,t)\in {\mathbb {R}}^1\times (0,T) , \end{aligned}$$

(3.5)

or

$$\begin{aligned} \frac{ \partial \overline{v} }{\partial t} - {\mathcal {A}}(t) \overline{v} + r_G\, \overline{v} = \overline{f}(x,t) \quad \hbox {for } (x,t)\in {\mathbb {R}}^1\times (0,T) , \end{aligned}$$

(3.6)

respectively, where $\underline{f}( \,\cdot \, ,t)\le \overline{f}( \,\cdot \, ,t)$ a.e. in ${\mathbb {R}}^1,$ at every time $t\in (0,T).$

Proof of Lemma 3.1

  We subtract Eq. (3.5) from (3.6), thus obtaining an analogous equation for the difference $w = \overline{v} - \underline{v}$ with the right-hand side equal to $g(x,t) = \overline{f}(x,t) - \underline{f}(x,t)\ge 0$ for a.e. $(x,t)\in {\mathbb {R}}^1\times (0,T).$ Then the desired result, $w( \,\cdot \, ,t)\ge 0$ a.e. in ${\mathbb {R}}^1,$ at every time $t\in [0,T),$ follows from Friedman [22, Chapt. 2, Sect. 4, Theorem 9, p. 43], cf. also [22, Chapt. 2, Sect. 6, Theorem 16, p. 52]. $\square$

We now give simple examples of super- and subsolutions of problem (2.13), (2.2).

Example 3.2

(Super- and subsolutions) Let us define the function

$$\begin{aligned} V(x,t) = K\, \mathrm {e}^{\lambda t} \left( \mathrm {e}^{\kappa x} + \mathrm {e}^{-\kappa x} \right) = 2K\, \mathrm {e}^{\lambda t}\cdot \cosh (\kappa x) \quad \hbox { for }\, (x,t)\in {\mathbb {R}}^1\times [0,T] , \end{aligned}$$

(3.7)

where $\kappa \in {\mathbb {R}}$ is a constant satisfying $1\le \kappa < \mu /2,$ and $K, \lambda \in {\mathbb {R}}$ with $K\ge 1$ and $\lambda \ge 0$ are some other constants (large enough) to be determined below:

The left-hand side of Eq. (2.13) with $v=V$ becomes

$$\begin{aligned}&\text {l.h.s.}(x,t) = \frac{\partial V}{\partial t} - {\mathcal {A}}(t) V + r_G\, V(x,t)\nonumber \\&\quad = \lambda \, V(x,t) - \frac{1}{2}\, \kappa ^2 [\sigma (t)]^2\, V(x,t)\nonumber \\&\qquad -\, \kappa \left[ q_S(t) - \gamma _S(t) - \frac{1}{2}\, [\sigma (t)]^2\right] \cdot \frac{ \mathrm {e}^{\kappa x} - \mathrm {e}^{-\kappa x} }{ \mathrm {e}^{\kappa x} + \mathrm {e}^{-\kappa x} }\, V(x,t) + r_G\, V(x,t)\nonumber \\&\quad \ge \left[ (\lambda + r_G) - \frac{1}{2}\, \kappa ^2 [\sigma (t)]^2 \right] V(x,t) - \kappa \cdot \left| q_S(t) - \gamma _S(t) - \frac{1}{2}\, [\sigma (t)]^2 \right| \, V(x,t)\nonumber \\&\quad \ge \left\{ (\lambda + r_G) - \genfrac{}{}{}1{1}{2}\, \kappa ^2 \Vert \sigma \Vert _{ L^{\infty }(0,T) }^2 - \kappa \left[ \Vert q_S - \gamma _S\Vert _{ L^{\infty }(0,T) } + \genfrac{}{}{}1{1}{2}\, \Vert \sigma \Vert _{ L^{\infty }(0,T) }^2 \right] \right\} V(x,t)\nonumber \\&\qquad \hbox {for } (x,t)\in {\mathbb {R}}^1\times [0,T] . \end{aligned}$$

(3.8)

As usual, $\Vert \sigma \Vert _{ L^{\infty }(0,T) }$ stands for the supremum norm of a continuous function $\sigma :\,[0,T]\rightarrow {\mathbb {R}}.$

On the other hand, the right-hand side of Eq. (2.13) with $v=V$ becomes

$$\begin{aligned} \text {r.h.s.}(x,t)&= G(V(x,t); x,t) = G(0; x,t) + \left[ G(V(x,t); x,t) - G(0; x,t)\right] \nonumber \\&\le C_0\, \exp (|x|) + L_G\, V(x,t)\le C_0'\, V(x,t) \quad \hbox {for } (x,t)\in {\mathbb {R}}^1\times [0,T] ,\nonumber \\&\quad \hbox {where }\ C_0'{\mathop {=}\limits ^{\mathrm{{def}}}}({C_0}/{K}) + L_G \quad ({} > 0) , \end{aligned}$$

(3.9)

and we have taken advantage of inequalities (2.10) and (2.11) in Hypotheses (G2) and (G4), respectively. Subtracting Eq. (3.9) from (3.8) we arrive at

$$\begin{aligned}&\text {l.h.s.}(x,t) - \text {r.h.s.}(x,t) = \frac{\partial V}{\partial t} - {\mathcal {A}}(t) V + r_G\, V(x,t) - G(V(x,t); x,t) \nonumber \\&\quad \ge \Bigl \{ \lambda {} + r_G - \genfrac{}{}{}1{1}{2}\, \kappa ^2 \Vert \sigma \Vert _{ L^{\infty }(0,T) }^2 \nonumber \\&\qquad - \kappa \left[ \Vert q_S - \gamma _S\Vert _{ L^{\infty }(0,T) } + \genfrac{}{}{}1{1}{2}\, \Vert \sigma \Vert _{ L^{\infty }(0,T) }^2 \right] - C_0'\Bigr \} V(x,t)\ge 0 \nonumber \\&\qquad \quad \hbox {for } (x,t)\in {\mathbb {R}}^1\times [0,T] , \end{aligned}$$

(3.10)

provided $\lambda \in {\mathbb {R}}$ satisfies

$$\begin{aligned} \lambda + r_G\ge \Lambda&{\mathop {=}\limits ^{\mathrm{{def}}}}{} \genfrac{}{}{}1{1}{2}\, \kappa ^2 \Vert \sigma \Vert _{ L^{\infty }(0,T) }^2\nonumber \\&\quad + \kappa \left[ \Vert q_S - \gamma _S\Vert _{ L^{\infty }(0,T) } + \genfrac{}{}{}1{1}{2}\, \Vert \sigma \Vert _{ L^{\infty }(0,T) }^2 \right] + C_0'\quad ({} > 0) . \end{aligned}$$

(3.11)

Recalling our Hypothesis (v$_{\mathbf {0}}$) with ineq. (2.17) on the growth of the initial condition, we take the constant $K\in {\mathbb {R}}$ such that $K\ge \max \{ 1,\, C_h\}$ which guarantees also

$$\begin{aligned} |v_0(x)|= & {} |h(\mathrm {e}^x)|\le C_h\cdot \exp (|x|) \le C_h\cdot \exp (\kappa |x|)\nonumber \\\le & {} V(x,0) = K\, \left( \mathrm {e}^{\kappa x} + \mathrm {e}^{-\kappa x} \right) = 2K\cdot \cosh (\kappa x) \quad \hbox {for}\ (x,t)\in {\mathbb {R}}^1\times [0,T] . \nonumber \\ \end{aligned}$$

(3.12)

It follows that the function $V:\,{\mathbb {R}}^1\times [0,T]\rightarrow {\mathbb {R}}$ defined in Eq. (3.7) is a supersolution of problem (2.13), (2.2).

Analogous arguments show that the function ${}- V:\,{\mathbb {R}}^1\times [0,T]\rightarrow {\mathbb {R}}$ is a subsolution of problem (2.13), (2.2). Notice that in this case, ineq. (3.9) has to be replaced by

$$\begin{aligned}&G(-V(x,t); x,t) = G(0; x,t) + \left[ G(-V(x,t); x,t) - G(0; x,t)\right] \nonumber \\&\quad \ge {}- C_0'\, V(x,t) \quad \hbox {for}\ (x,t)\in {\mathbb {R}}^1\times [0,T] . \end{aligned}$$

(3.13)

$\square$

3.2 Construction of monotone iterations

Recalling Lemma 3.1, our Hypothesis (v$_{\mathbf {0}}$) with ineq. (2.17), and applying Example 3.2 with $\kappa = 1$ $(< \mu /2),$ we are now ready to construct a for calculating a (weak) solution $v:\,{\mathbb {R}}^1\times (0,T)\rightarrow {\mathbb {R}}$ to the initial value problem (2.13), (2.2). We start by setting $\kappa = 1$ and fixing the constants $K, \lambda \in {\mathbb {R}}$ with $K\ge 1$ and $\lambda + r_G\ge \Lambda$ $(\, > 0)$ large enough, such that both inequalities, (3.11) and (3.12), are valid. It follows from Eq. (3.7) and ineq. (3.10) that the function $u_0:\,{\mathbb {R}}^1\times (0,T)\rightarrow {\mathbb {R}}$ defined by

$$\begin{aligned} u_0(x,t) = K\, \mathrm {e}^{\lambda t} \left( \mathrm {e}^x + \mathrm {e}^{- x} \right) = 2K\, \mathrm {e}^{\lambda t}\cdot \cosh (x) \quad \hbox {for}\ (x,t)\in {\mathbb {R}}^1\times [0,T] , \end{aligned}$$

(3.14)

is a supersolution of problem (2.13), (2.2). We remark that ${}- u_0$ happens to be a subsolution of this problem, by Example 3.2 with $\kappa = 1,$ as well. The first iterate, $u_1:\,{\mathbb {R}}^1\times [0,T]\rightarrow {\mathbb {R}},$ is constructed as the (weak) solution $u_1$ to the following analogue of the inhomogeneous linear initial value problem (2.20), (2.2):

$$\begin{aligned} \frac{\partial u_1}{\partial t} - {\mathcal {A}}(t) u_1 + r_G\, u_1&{} = G(u_0(x,t); x,t) \nonumber \\&\left( {}\le C_0'\, u_0(x,t) \right)&\quad \hbox { for }\, (x,t)\in {\mathbb {R}}^1\times [0,T] ; \end{aligned}$$

(3.15)

$$\begin{aligned} u_1(x,0)&{}= v_0(x){\mathop {=}\limits ^{\mathrm{{def}}}}h(\mathrm {e}^x)&\quad \hbox { for }\, x\in {\mathbb {R}}^1 . \end{aligned}$$

(3.16)

Since $u_0$ is a (weak) supersolution of problem (2.13), (2.2), by Example 3.2 with $\kappa = 1,$ we may apply Lemma 3.1 to conclude that $u_1(x,t)\le u_0(x,t)$ holds for a.e. $(x,t)\in {\mathbb {R}}^1\times (0,T).$ In addition, making use of Eq. (3.13), we get also

$$\begin{aligned} \begin{aligned} {}- u_0(x,t)\le u_1(x,t)\le u_0(x,t) \quad \left( {} = K\, \mathrm {e}^{\lambda t} \left( \mathrm {e}^{\kappa x} + \mathrm {e}^{-\kappa x} \right) \right) \\ \quad \hbox { for }\, (x,t)\in {\mathbb {R}}^1\times [0,T] . \end{aligned} \end{aligned}$$

(3.17)

Our next step is the following induction hypothesis. Let us assume that, for some integer $m\ge 1,$ in addition to $u_0$ and $u_1$ above, we have already constructed the first $(m+1)$ functions $u_0, u_1, u_2,\ldots , u_m:\,{\mathbb {R}}^1\times [0,T]\rightarrow {\mathbb {R}}$ with the following properties:

• (a) Every function $u_j:\,{\mathbb {R}}^1\times [0,T]\rightarrow {\mathbb {R}}$; $j = 0,1,2,\ldots , m,$ is Lebesgue-measurable and continuous in time as a function $u_j:\,[0,T]\rightarrow H,$ i.e., $u_j\in C([0,T]\rightarrow H).$

• (b) The inequalities

  $$\begin{aligned} \begin{aligned} {}- u_0(x,t)\le u_j(x,t)\le u_{j-1}(x,t)\le u_0(x,t) \quad \hbox { for }\, (x,t)\in {\mathbb {R}}^1\times [0,T] , \end{aligned} \end{aligned}$$

  (3.18)

  are valid for every $j = 1,2,3,\ldots , m.$

• (c) For each $j = 1,2,3,\ldots , m,$ the function $u_j:\,{\mathbb {R}}^1\times [0,T]\rightarrow {\mathbb {R}}$ is the (weak) solution to the following analogue of the inhomogeneous linear initial value problem (2.20), (2.2); cf. problem (3.15), (3.16) above:

  $$\begin{aligned} \frac{\partial u_j}{\partial t} - {\mathcal {A}}(t) u_j + r_G\, u_j&{} = G(u_{j-1}(x,t); x,t) \nonumber \\&\left( {}\le C_0'\, u_0(x,t) \right)&\quad \hbox { for }\, (x,t)\in {\mathbb {R}}^1\times [0,T] ; \end{aligned}$$

  (3.19)
  $$\begin{aligned} u_j(x,0)&{}= v_0(x){\mathop {=}\limits ^{\mathrm{{def}}}}h(\mathrm {e}^x)&\quad \hbox { for }\, x\in {\mathbb {R}}^1 . \end{aligned}$$

  (3.20)

In our last step (induction on the index $m\ge 1$) we construct the $(m+1)$-st iterate, $u_{m+1}:\,{\mathbb {R}}^1\times [0,T]\rightarrow {\mathbb {R}},$ to be the (weak) solution $u_{m+1}$ to the following problem; cf. problem (3.15), (3.16):

$$\begin{aligned}&\frac{\partial u_{m+1}}{\partial t} {} - {\mathcal {A}}(t) u_{m+1} + r_G\, u_{m+1} = G(u_m(x,t); x,t) \nonumber \\&\quad \left( {}\le C_0'\, u_0(x,t) \right) \quad \hbox {for}\ (x,t)\in {\mathbb {R}}^1\times [0,T] ; \end{aligned}$$

(3.21)

$$\begin{aligned}&u_{m+1}(x,0) = v_0(x){\mathop {=}\limits ^{\mathrm{{def}}}}h(\mathrm {e}^x) \quad \hbox {for}\ x\in {\mathbb {R}}^1 . \end{aligned}$$

(3.22)

By arguments analogous to those used in the construction of $u_1$ from $u_0$ above, we conclude that $u_{m+1}$ exists and satisfies

$$\begin{aligned} {}- u_0(x,t)\le u_{m+1}(x,t)\le u_m(x,t)\le u_0(x,t) \quad \hbox {for}\ (x,t)\in {\mathbb {R}}^1\times [0,T] . \end{aligned}$$

(3.23)

Here, we have used our monotonicity hypothesis (G3) to conclude that $u_m\le u_{m-1}$ a.e. in ${\mathbb {R}}^1\times [0,T]$ entails $G(u_m(x,t); x,t)\le G(u_{m-1}(x,t); x,t)$ for a.e. $(x,t)\in {\mathbb {R}}^1\times [0,T].$ Finally, we get $u_{m+1}\in C([0,T]\rightarrow H).$ This concludes the construction of the desired iterates.

Remark 3.3

(Inhomogeneous linear problem). An for calculating the weak solution $u_{m+1}:\,{\mathbb {R}}^1\times [0,T]\rightarrow {\mathbb {R}}$ to problem (3.21), (3.22) will be described later in Sect. 5, Corollary 5.2. Numerical methods for computing this solution will be discussed in Remark 5.4, as well. $\square$

As an obvious consequence of our construction we conclude that, in addition to $u_0,$ also each iterate $u_j$; $j = 1,2,3,\ldots ,$ is a (weak) supersolution of problem (2.13), (2.2).

Standard application of Lebesgue’s monotone (or dominated) convergence theorem yields the monotone pointwise convergence $u_m(x,t)\searrow v(x,t)$ as $m\nearrow \infty$ for a.e. $(x,t)\in {\mathbb {R}}^1\times [0,T],$ as well as the $L^2$-type norm convergence $\Vert u_m(t) - v(t)\Vert _H\searrow 0$ as $m\nearrow \infty ,$ for a.e. $t\in (0,T),$ where $v(t)\equiv v( \,\cdot \, ,t)\in H$ satisfies

$$\begin{aligned} {}- u_0(x,t)\le v(x,t)\le u_m(x,t)\le u_0(x,t) \quad \hbox {for}\ (x,t)\in {\mathbb {R}}^1\times [0,T] . \end{aligned}$$

(3.24)

Furthermore, we get another $L^2$-type norm convergence in the Lebesgue(–Hilbert) space $L^2([0,T]\rightarrow H),$ namely,

$$\begin{aligned} \Vert u_m-v\Vert _{ L^2([0,T]\rightarrow H) } {\mathop {=}\limits ^{\mathrm{{def}}}}\left( \int _0^T \Vert u_m(t) - v(t)\Vert _H^2 \,{\mathrm {d}}t \right) ^{1/2} \;\searrow \; 0 \quad \hbox {as}\ m\nearrow \infty . \end{aligned}$$

(3.25)

We combine this result with Theorem 1.2 (and its proof) in Pazy [40, Chapt. 6, §6.1, pp. 184–185] to conclude that $v\in C([0,T]\rightarrow H)$ and $v:\,{\mathbb {R}}^1\times [0,T]\rightarrow {\mathbb {R}}$ is a mild solution to the initial value problem (2.13), (2.2). Finally, applying the well-known results from Evans [19], Chapt. 7, §1.1, p. 352, or Lions [35], Chapt. IV, §1, p. 44, or [36], Chapt. III, eq. (1.11), p. 102, we find out that v is a (weak) solution to problem (2.13), (2.2).

We summarize the most important results from our monotone iteration scheme (3.14)–(3.23) for problem (2.13), (2.2) in the following theorem. Precise definitions of the Hölder spaces used in this theorem, $H^{\theta , \theta /2}\bigl ( D_{1+1}^{(T)} \bigr )$ (a local Hölder space) and $H^{2+\theta , 1 + (\theta /2)}(\overline{Q'})$ over the parabolic domain $D_{1+1}^{(T)} = {\mathbb {R}}^1\times (0,T)$ (an open strip in ${\mathbb {R}}^1\times {\mathbb {R}}$) and its compact subset $\overline{Q'} = [a,b]\times [\tau , T']\subset D_{1+1}^{(T)}$ (a compact rectangle), respectively, with $\theta \in (0,1),$ $-\infty< a< b < +\infty ,$ and $0< \tau< T'< T,$ can be found in Ladyzhenskaya et al. [32, Chapt. I, §1, pp. 7–8].

Theorem 3.4

(Monotone iterations). Let $v_0\in H$ obey Hypothesis (v$_{\mathbf {0}})$ with ineq. (2.17). Then the monotone iterations $u_0\ge u_1\ge \cdots \ge u_{j-1}\ge u_j\ge \cdots \ge {}- u_0 ,$ described in items (a), (b), and (c) above,  converge in the Lebesgue(–Hilbert) space $L^2([0,T]\rightarrow H)$ to a function $v:\,{\mathbb {R}}^1\times (0,T)$ according to formula (3.25). The limit function,  $v\in L^2([0,T]\rightarrow H),$ is a (weak) solution to problem (2.13),  (2.2). Furthermore,  there is a constant $\theta \in (0,1)$ such that $u_m\in H^{\theta , \theta /2}\bigl ( D_{1+1}^{(T)} \bigr )$ holds for every $m = 1,2,3,\ldots ,$ and $v\in H^{\theta , \theta /2}\bigl ( D_{1+1}^{(T)} \bigr ),$ as well.

Finally,  assume that the function

Then we get even $u_m\in H^{2+\theta , 1 + (\theta /2)}\bigl ( D_{1+1}^{(T)} \bigr )$ for every $m = 1,2,3,\ldots ,$ together with $v\in H^{2+\theta , 1 + (\theta /2)}\bigl ( D_{1+1}^{(T)} \bigr ) ,$ where the convergence $u_m\rightarrow v$ holds in the norm of the Hölder space $H^{2+\theta ', 1 + (\theta '/2)}(\overline{Q'})$ over any compact rectangle $\overline{Q'} = [a,b]\times [\tau , T']\subset D_{1+1}^{(T)}$ in the open strip $D_{1+1}^{(T)} = {\mathbb {R}}^1\times (0,T),$ with $-\infty< a< b < +\infty$ and $0< \tau< T'< T,$ and with any Hölder exponent $\theta '\in (0,\theta ).$ In particular,  each function $u_m$ $(m = 1,2,3,\ldots )$ is a strong (classical) solution to problem (3.19),  (3.20),  whereas v is a strong (classical) solution to problem (2.13),  (2.2).

Proof

The convergence of our monotone iterations in the Lebesgue(–Hilbert) space $L^2([0,T]\rightarrow H)$ has been established in Eq. (3.25) before this theorem. Standard regularity theory in Ladyzhenskaya et al. [32, Chapt. III, §10], Theorem 10.1 on p. 204, guarantees that there is a constant $\theta \in (0,1)$ such that $u_m\in H^{\theta , \theta /2}\bigl ( D_{1+1}^{(T)} \bigr )$ holds for every $m = 1,2,3,\ldots .$ Another regularity result in [32], Chapt. IV, §15, Theorem 15.1, pp. 405–406, then yields $u_m\in H^{2+\theta , 1 + (\theta /2)}\bigl ( D_{1+1}^{(T)} \bigr )$ for every $m = 1,2,3,\ldots ,$ with the Hölder norms $\Vert u_m\Vert _{ H^{2+\theta , 1 + (\theta /2)}(\overline{Q'}) }$ being uniformly bounded, whenever $\overline{Q'} = [a,b]\times [\tau , T']\subset D_{1+1}^{(T)}$ is a compact rectangle in the open strip $D_{1+1}^{(T)} = {\mathbb {R}}^1\times (0,T),$ with $-\infty< a< b < +\infty$ and $0< \tau< T'< T.$ A simple argument based upon the compactness theorem by Arzelà–Ascoli in $H^{2+\theta , 1 + (\theta /2)}(\overline{Q'})$ shows that $\Vert v\Vert _{ H^{2+\theta , 1 + (\theta /2)}(\overline{Q'}) } < +\infty ,$ as well, together with $v\in H^{2+\theta , 1 + (\theta /2)}\bigl ( D_{1+1}^{(T)} \bigr ) .$ In this context, Theorem 3.1 on pp. 983–984 in Sattinger [42] applies to our monotone iteration scheme (3.14)–(3.23) for problem (2.13), (2.2). $\square$

Two additional important results in Sattinger’s work [42], Theorems 3.3 and 3.4 on p. 986, can be applied to our monotone iteration scheme, (3.14)–(3.23), as well. We remark that the in [42] are formulated for nonlinear elliptic and parabolic boundary value problems in a bounded spatial domain $\Omega \subset {\mathbb {R}}^N.$ Nevertheless, they are applicable also to our terminal value problem (1.1), (1.2) (as evidenced by our proof of Theorem 3.4 above), even though we have no boundary conditions for the stock price $S\in (0,\infty )$; more precisely, not for the logarithmic stock price $x = \log ~S\in {\mathbb {R}}$ that we use throughout our present article; cf. [6, §3.3, p. 50]. In analogy with this scheme, (3.14)–(3.23), for problem (2.13), (2.2), which has resulted in the monotone (pointwise) decreasing sequence of iterations $u_0\ge u_1\ge \cdots \ge u_{j-1}\ge u_j\ge \cdots \ge {}- u_0$ a.e. in ${\mathbb {R}}^1\times (0,T)$ in Theorem 3.4 above, we are able to construct another monotone iteration scheme for problem (2.13), (2.2) that is monotone (pointwise) as follows:

Remark 3.5

(Increasing monotone iterations). We set

$$\begin{aligned} w_0(x,t) = {}- u_0(x,t) = {}- K\, \mathrm {e}^{\lambda t} \left( \mathrm {e}^x + \mathrm {e}^{- x} \right) \quad \hbox {for}\ (x,t)\in {\mathbb {R}}^1\times [0,T] , \end{aligned}$$

(3.26)

and recall that $w_0 = {}- u_0$ is a subsolution of problem (2.13), (2.2), by Example 3.2 with $\kappa = 1$ and some constants $K\ge 1$ and $\lambda \ge 0$ large enough, as specified in this example. We define the first iteration $w_1:\,{\mathbb {R}}^1\times [0,T]\rightarrow {\mathbb {R}}$ by replacing the functions $u_0$ and $u_1$ in Eqs. (3.15) and (3.16) by $w_0$ and $w_1,$ respectively. Thus, the inequalities in (3.17) have to be reversed:

$$\begin{aligned} {}- u_0(x,t) = w_0(x,t)\le w_1(x,t)\le u_0(x,t) \quad \hbox {for}\ (x,t)\in {\mathbb {R}}^1\times [0,T] . \end{aligned}$$

(3.27)

In order to construct the j-th iteration $w_j:\,{\mathbb {R}}^1\times [0,T]\rightarrow {\mathbb {R}}$ from $w_{j-1}:\,{\mathbb {R}}^1\times [0,T]\rightarrow {\mathbb {R}}$; $j = 1,2,3,\ldots ,$ we replace the pair $(u_{j-1}, u_j)$ by $(w_{j-1}, w_j)$ in (3.18)–(3.23), thus arriving at the desired monotone (pointwise) increasing sequence of iterations $w_0 = - u_0\le w_1\le \cdots \le w_{j-1}\le w_j\le \cdots \le u_0$ a.e. in ${\mathbb {R}}^1\times (0,T).$

The remaining part of Theorem 3.4 remains valid for $w_m$ in place of $u_m$; $m = 1,2,3,\ldots ,$ both, with and without hypothesis (G1’), including the monotone convergence $w_m\rightarrow v$ in the Lebesgue space $L^2([0,T]\rightarrow H)$ to the function $v:\,{\mathbb {R}}^1\times (0,T)$ described in Theorem 3.4 by formula (3.25) and $u_m\rightarrow v$ in the norm of the Hölder space $H^{2+\theta ', 1 + (\theta '/2)}(\overline{Q'}) ,$ respectively. $\square$

Numerical approximation methods for semilinear parabolic systems based on a standard finite difference discretization that produces a monotone iteration scheme can be found in the monographs by Leung [33, 34]. Accelerated Monotone Convergence of this iteration scheme is treated in [33], Chapt. VI, §6.3, pp. 292–300, in the context of a two-point boundary value problem followed immediately by $L^2$-Convergence for finite-difference solutions in two dimensional domains in §6.4, pp. 300–323, which can be viewed as a numerical implementation of our $L^2$-type norm convergence result in Eq. (3.25). A full, large parabolic system is treated numerically in [33, Chapt. VII], §7.1–§7.3, pp. 327–357. Further applications of monotone iteration schemes appear in the newer monograph by Leung [34], Chapt. 3, §3.2–§3.3, pp. 210–257.

An important question in this (cf. Sattinger [42]) is the in Eq. (3.25) above. A suitable norm for estimating this convergence in the Banach space $C([0,T]\rightarrow H)\hookrightarrow L^2((0,T)\rightarrow H)$ is the “weighted” supremum norm

$$\begin{aligned} {\left| \left| \left| u \right| \right| \right| }_{ L^{\infty }(0,T) }{\mathop {=}\limits ^{\mathrm{{def}}}}\sup _{0\le t\le T} \left( \mathrm {e}^{-\omega t}\, \Vert u(\,\cdot ,\, t)\Vert _H \right) \quad ({} < \infty ) , \end{aligned}$$

(3.28)

where $\omega > 0$ is a sufficiently large positive real number that guarantees the estimate

$$\begin{aligned}&{\left| \left| \left| u_{m+1} - u_j \right| \right| \right| }_{ L^{\infty }(0,T) } \le c_{\omega }\, {\left| \left| \left| u_m - u_{j-1} \right| \right| \right| }_{ L^{\infty }(0,T) } \le \cdots \\&\quad \le c_{\omega }^j\, {\left| \left| \left| u_{m+1-j} - u_0 \right| \right| \right| }_{ L^{\infty }(0,T) } \quad \hbox {for every}\ j = 1,2,3,\ldots , m; \end{aligned}$$

$m = 1,2,3,\ldots ,$ with some positive constant $c_{\omega } < 1.$ We refer to Pazy [40, Chapt. 2, §2.2, pp. 44–45] for details about how to estimate the constant $\omega > 0$ (from below) from the values of $r_G = r + L_{\tilde{F}}$ and $L_G$ in Eq. (2.6) and ineq. (2.10), respectively. Then also the estimate

$$\begin{aligned}&{\left| \left| \left| v - u_j \right| \right| \right| }_{ L^{\infty }(0,T) } \le c_{\omega }\, {\left| \left| \left| v - u_{j-1} \right| \right| \right| }_{ L^{\infty }(0,T) } \le \cdots \\&\quad \le c_{\omega }^j\, {\left| \left| \left| v - u_0 \right| \right| \right| }_{ L^{\infty }(0,T) } \quad \hbox {for every}\ j = 1,2,3,\ldots , \end{aligned}$$

follows immediately by letting $m\rightarrow \infty$ above.

4 Numerical methods: finite differences/elements compared to Monte Carlo

The Monte Carlo method enjoys broad popularity especially among the specialists in Probability interested in numerical simulations. Since its introduction into the toolbox of computational finance in 1977 by Boyle et al. [11], Monte Carlo simulation has continued to gain acceptance as a computational method of pricing more and more sophisticated options and as a risk management tool, as well. The Monte Carlo method is characterized by its very high flexibility and its ability to process a multidimensional problem. One of the main strengths of the Monte Carlo method is that it does not require discretization, and is therefore insensitive to space dimension. Consequently, the more dimensions the problem to be solved has, the more competitive the Monte Carlo method is when compared to methods based on space discretization. However, the rate of convergence of the method, apparently dependent on the space dimension of the problem (often referred to as “the curse of dimensionality” in case of a “high” space dimension, which is better known as degeneracy properties of likelihood ratios (LRs)), is known to be quite low; cf. Glynn and Rubinstein [23]. An important identified limitation of the Monte Carlo method is its inability to treat discretization of nonlinearintegral equations

directly. However, iteration methods that require the solution of a linear problem at each iteration step (such as ours in Sect. 3) have been studied intensively in the past, see e.g. Plotnikov [41]. The Monte Carlo method is applied there to solve linear problems in order to construct a sequence of iterates that are expected to converge to the solution of the given nonlinear integral equation. This is the strategy we have followed also in the work reported here (in Sect. 3).

Some recent works open perspectives to overcome the difficulties with discretization of nonlinear integral equations. Significant progress has been made on the construction of methods based on a . As the nonlinearities to be treated in numerical discretizations of PDEs in Finance are mostly Lipschitz-continuous (Lipschitzian, for short), often piecewise linear and continuous, the most efficient Monte Carlo algorithms focus on treating precisely such Lipschitzian nonlinearities. Among the most interesting earlier pioneering works, we would like to mention the article by Henry-Labordère in [25]. Here, the author approximates a simple piecewise linear nonlinearity by a polynomial in order to be able to take advantage of the probability-based computational technique, a branching diffusion process. This technique was first introduced by McKean [38] and Skorokhod [43] to provide a probabilistic representation for the solution of the Kolmogorov–Petrovskii–Piskunov PDE ( equation, for short) and, more generally, of semilinear PDEs whose nonlinearity is a power series with nonnegative coefficients from the interval [0, 1]. Since the KPP equation in [38, 43] has only a quadratic or cubic nonlinearity, numerical approximation using branching diffusion process is quite efficient. However, for applications to the in finance treated in [25], the nonlinearity that one is interested in is a polynomial. It is therefore interesting to investigate methods that can treat directly monotone nonlinearities that are not necessarily polynomial. As to the speed and precision of the PDE and Monte-Carlo methods, a brief comparison is given in the work by Loeper and Pironneau [37], together with a mixed PDE/Monte-Carlo method that provides better results for the Heston stochastic volatility model. An entirely different approach to Heston’s model [28], based on “orthogonal series expansion”, is employed in Alziary and Takáč [2, Sect. 11 (Appendix), §11.1–§11.2], pp. 48–50, in Takáč [44], and in the numerical simulations by Baustian et al. [7]. Replacing a piecewise linear nonlinearity by a polynomial introduces a significant error into the algorithm in [25].

These branching diffusion processes provide a generalization of the Feynman–Kac theory for nonlinear PDEs and lead to a representation of the solutions in the form of expectation and thus allow to use Monte Carlo methods for simulations. A significant progress in the application of branching diffusion process was made in Henry-Labordère et al. [27]. Therefore, more recent studies prefer to work directly with an arbitrary Lipschitzian nonlinearity when applying branching diffusion process; cf. Bouchard et al. [10] and Henry-Labordère et al. [26]. More general (non-Lipschitzian) nonlinearities are approximated by (optimal) Lipschitzian nonlinearities on a compact set, if necessary. In this situation, a standard Picard’s iteration scheme (from the proof of the Picard–Lindelöf theorem; see e.g. Hartman [24, Chapt. II, §1, pp. 8–10]) can be used to approximate the desired solution $v\in L^2([0,T]\rightarrow H)$ to problem (2.13), (2.2), i.e., the limit function v in our Theorem 3.4. Also the approximation error and the rate of convergence are estimated by classical arguments.

Before switching entirely to a discussion of typically “analytical” numerical discretizations, we would like to mention also some important studies where the two methods, i.e., Monte Carlo and finite difference/element algorithms, are either combined ([37]) or compared against each other in numerical simulations ([6, Sect. 4] and [37]). A suitable combination of these two methods performed on the “right” problem can lead to significant improvements of the “hybrid” algorithm by Loeper and Pironneau [37] in both, precision and speed. On the other hand, the work in Baustian et al. [6] treats the same subject as does our present work, with an [6, Eq. (6), p. 48] closely related to ours in Eqs. (2.20), (2.2). An interesting comparison of several kinds of Numerical Methods for the linear initial value problem (2.20), (2.2) can be found in that work [6, Sect. 4, pp. 52–54]. In fact, even a quantitative comparison of Monte Carlo and finite difference methods with the exact formula (for the “diffusion (or heat) equation”) is provided there [6, Figs. 2–4, p. 53]. In [6, Sect. 3], Monte Carlo and finite difference/element algorithms are applied separately and then compared against each other in numerical simulations [6, Sect. 4, pp. 52–54]. In general, the comparison can go “both ways”; one has to choose among Monte Carlo, finite difference/element, and “hybrid” algorithms according to the particular problem that should be discretized for a numerical treatment.

In the previous section (Sect. 3) we have presented a method suitable for direct applications of finite difference/element algorithms to of type (2.20), (2.2) with the inhomogeneity (on the right-hand side in (2.20)) given by the term $f(x,t) = G(v(x,t); x,t)$ from Eq. (2.13). Our functional-analytic approach with an arbitrary Lipschitzian nonlinearity suggests an analogous Picard iteration scheme combined with a finite difference (or similar) discretization method for the inhomogeneous linear parabolic initial value problem (2.20), (2.2). Solving this problem is especially easy if all functions $\sigma :\,[0,T]\rightarrow (0,\infty )$ and $q_S, \gamma _S:\,[0,T]\rightarrow {\mathbb {R}}$ in Hypotheses (BS1) and (BS2), respectively, are constants independent from time $t\in [0,T],$ i.e., $\sigma \in (0,\infty )$ and $q_S, \gamma _S\in {\mathbb {R}}.$ Our choice of the constant $r_G\in {\mathbb {R}}_+$ guarantees that Hypothesis (G3) is valid: For almost every pair $(x,t)\in {\mathbb {R}}^1\times (0,T),$ the function $G(\,\cdot \, ; x,t):\,v\mapsto G(v; x,t):\,{\mathbb {R}}\rightarrow {\mathbb {R}}$ is monotone increasing. Under these Hypotheses, (BS1) and (BS2), supplemented by Hypothesis (f1) in place of Hypotheses (G2) and (G5) used for the special inhomogeneity $f(x,t) = G(v(x,t); x,t)$ on the right-hand side in Eq. (2.13), in the next section (Sect. 5) we give an for the solution v(x, t) to the inhomogeneous linear parabolic initial value problem (2.20), (2.2). This formula is obtained in terms of a linear integral operator whose kernel is expressed through a “modification” of the standard “diffusion (or heat) kernel”, cf. Theorem 5.1 below.

5 Numerical methods: the semilinear parabolic problem

In order to calculate (and compute numerically) the monotone iterations $u_0\ge u_1\ge \cdots \ge u_{j-1}\ge u_j\ge \cdots \ge {}- u_0 ,$ treated in Theorem 3.4, we need to solve the inhomogeneous linear parabolic initial value problem (2.13), (2.2) under the assumption that the inhomogeneity $f(x,t) = G(v(x,t); x,t)$ on the right-hand side in Eq. (2.13) is known. In other words we wish to supplement our main result, Theorem 3.4, by a result on solving the initial value problem (2.20), (2.2). We will obtain a unique weak (or mild) solution $v:\,[0,T]\rightarrow H$ under Hypothesis (f1), that is to say, $f:\,[0,T]\rightarrow H$ is a continuous function, where $f:\,{\mathbb {R}}^1\times [0,T]\rightarrow {\mathbb {R}}:\,(x,t)\mapsto f(x,t)$ satisfies $f(t)\equiv f( \,\cdot \, ,t)\in H$ for every $t\in [0,T].$

Furthermore, our special inhomogeneity choice of $f(x,t) = G(v(x,t); x,t)$ satisfies even the stronger hypothesis (2.21) on the Hölder continuity of the function $f:\,[0,T]\rightarrow H$ with the Hölder exponent $\vartheta _f\in (0,1),$ see also ineq. (2.22). As a consequence, our mild solution $v:\,[0,T]\rightarrow H,$ which is continuous by definition, becomes a unique classical solution. Thus, from now on let us assume that $f:\,[0,T]\rightarrow H$ is a Hölder-continuous function with the Hölder exponent $\vartheta _f\in (0,1)$; cf. ineq. (2.21) and also ineq. (2.22) related to the linear initial value problem (2.20), (2.2).

To solve the initial value problem (2.20), (2.2), we consider the abstract differential equation (2.24) with the initial condition (2.25), i.e., $v(0) = v_0\in H$ at time $t=0.$ A unique mild solution $v\in C([0,T]\rightarrow H)$ to this abstract initial value problem is obtained by the following formula,

$$\begin{aligned} \begin{aligned} v(t) = {\mathfrak {T}}(t,0) v_0 + \int _0^t {\mathfrak {T}}(t,\tau )\, f(\tau ) \,{\mathrm {d}}\tau \quad \hbox {in }H \hbox { for }\, t\in [0,T] , \end{aligned} \end{aligned}$$

(5.1)

where ${\mathfrak {T}} = \{ {\mathfrak {T}}(t,s):\,H\rightarrow H:\,0\le s\le t\le T\}$ stands for the evolutionary family of bounded linear operators ${\mathfrak {T}}(t,s):\,H\rightarrow H\ (0\le s\le t\le T)$ that yield the unique mild solution

$$\begin{aligned} v(t) = {\mathfrak {T}}(t,s) v_s \quad \hbox {in } H \hbox { for } t\in [s,T] \end{aligned}$$

to the homogeneous abstract differential equation (with $f(t)\equiv 0$) on the interval [s, T] for each $s\in [0,T),$

$$\begin{aligned}&\frac{\partial v}{\partial t} - {\mathcal {A}}(t) v + r_G\, v = 0 \quad \hbox {in }H \hbox { for } t\in (s,T) ; \end{aligned}$$

(5.2)

$$\begin{aligned}&\hbox {with}\ v(s)= v_s\in H \quad \hbox { (the initial condition).} \end{aligned}$$

(5.3)

We refer the reader to the monograph by Pazy [40, Chapt. 5, §5.6], Theorem 6.8 on p. 164, for greater details on the evolutionary family ${\mathfrak {T}}$ for the homogeneous abstract initial value problem (5.2), (5.3) and to [40, Chapt. 5, §5.7], Theorem 7.1 on p. 168, for the variation-of-constants formula (5.1) applied to the inhomogeneous abstract initial value problem (2.24), (2.25) posed on the entire time interval [0, T].

In our theorem below, Theorem 5.1, we provide an explicit formula for the evolutionary family ${\mathfrak {T}}.$ Each operator ${\mathfrak {T}}(t,s):\,H\rightarrow H\ (0\le s < t\le T)$ turns out to be a (bounded) integral operator with the kernel ${\mathfrak {K}}(x,y;t,s) > 0$ described explicitly in terms of the standard “diffusion (or heat) kernel”

$$\begin{aligned}&{\mathfrak {G}}(x;t){\mathop {=}\limits ^{\mathrm{{def}}}}\frac{1}{\sqrt{4\pi t}}\, \exp \left( {}- \frac{x^2}{4t} \right) \quad ({} > 0) \quad \hbox {for}\ (x,t)\in {\mathbb {R}}^1\times (0,\infty ) ;\nonumber \\&\hbox {hence}, \hbox {we have }\quad \int _{-\infty }^{+\infty } {\mathfrak {G}}(x;t) \,{\mathrm {d}}x = 1 \quad \hbox {for every}\ t\in (0,\infty ) . \end{aligned}$$

(5.4)

We refer to John [29, Chapt. 7, §1(a), pp. 206–213] for greater details about the “diffusion (or heat) equation”. Next, we make use of the abbreviations

$$\begin{aligned} {\mathscr {S}}(t){\mathop {=}\limits ^{\mathrm{{def}}}}\frac{1}{2}\int _0^t [\sigma (\tau )]^2 \,{\mathrm {d}}\tau \quad \hbox { and }\quad {\mathscr {R}}(t){\mathop {=}\limits ^{\mathrm{{def}}}}\int _0^t \left( q_S(\tau ) - \gamma _S(\tau ) - \frac{1}{2}\, [\sigma (\tau )]^2 \right) \,{\mathrm {d}}\tau \end{aligned}$$

(5.5)

defined for every $t\in [0,T].$ By our Hypothesis (BS1), we have $\sigma _0 = \min _{t\in [0,T]} \sigma (t) > 0$ which entails

$$\begin{aligned} {\mathscr {S}}(t) - {\mathscr {S}}(s) = \frac{1}{2}\int _s^t [\sigma (\tau )]^2 \,{\mathrm {d}}\tau \ge \frac{1}{2}\, \sigma _0^2\, (t-s) > 0 \quad \hbox { whenever }\; 0\le s < t\le T . \end{aligned}$$

Theorem 5.1

(Homogeneous linear problem). Under the Hypotheses (BS1) and (BS2) stated in Sect. 2,  the evolutionary family ${\mathfrak {T}}$ on the Hilbert space H takes the following form : 

$$\begin{aligned}&v(t) = \left[ {\mathfrak {T}}(t,s) v_s\right] (x) = \int _{-\infty }^{+\infty } {\mathfrak {K}}(x,y;t,s)\, v_s(y) \,{\mathrm {d}}y\nonumber \\&\quad \hbox { for all }x\in {\mathbb {R}}^1, 0\le s < t\le T, \hbox { and }\, v_s\in H , \end{aligned}$$

(5.6)

with the kernel ${\mathfrak {K}}(x,y;t,s) > 0$ defined by

$$\begin{aligned}&{\mathfrak {K}}(x,y;t,s) \equiv {\mathfrak {K}}(x-y;t,s) \nonumber \\&\quad = \exp ( {}- r_G (t-s) )\cdot {\mathfrak {G}} \left( x-y + {\mathscr {R}}(t) - {\mathscr {R}}(s) ;\, {\mathscr {S}}(t) - {\mathscr {S}}(s) \right) \nonumber \\&\quad = \frac{ \exp ( {}- r_G (t-s) ) }{ \sqrt{ 4\pi [ {\mathscr {S}}(t) - {\mathscr {S}}(s) ] } }\, \exp \left( {}- \frac{( x-y + {\mathscr {R}}(t) - {\mathscr {R}}(s) )^2}{ 4 [ {\mathscr {S}}(t) - {\mathscr {S}}(s) ] } \right) \nonumber \\&\qquad \quad \hbox {for all }x,y\in {\mathbb {R}}^1 \hbox { and }0\le s < t\le T. \end{aligned}$$

(5.7)

Here,  the continuous function $v:\,[s,T]\rightarrow H$ defined in formula (5.6) for $t\in (s,T],$ with $v(s) = v_s\in H,$ is the unique mild solution of the homogeneous abstract initial value problem (5.2),  (5.3) on the interval [s, T]. In fact,  $v:\,[s,T]\rightarrow H$ is also a classical solution to this problem.

In fact, our existence and uniqueness result in this theorem corresponds to that in Pazy [40, Chapt. 5, §5.6], Theorem 6.8 on p. 164. However, since we wish to calculate the kernel ${\mathfrak {K}}(x,y;t,s)$ in formula (5.7) explicitly, we cannot use the abstract proof of [40, Chapt. 5, Theorem 6.8] directly. Thus, we prefer to calculate the evolutionary family ${\mathfrak {T}}$ directly by taking advantage of well-known results [29, Chapt. 7, §1(a), pp. 206–213] for the “diffusion equation”.

But, before giving the proof of Theorem 5.1, let us state an important corollary concerning the solvability of the inhomogeneous abstract initial value problem (2.24), (2.25).

Corollary 5.2

(Inhomogeneous linear problem). Under the Hypotheses (BS1) and (BS2) stated in Sect. 2,  and $f\in C([0,T]\rightarrow H),$ the initial value problem (2.24),  (2.25) possesses a unique mild solution $v:\,[0,T]\rightarrow H$ that is given by the variation-of-constants formula (5.1). If $f:\,[0,T]\rightarrow H$ obeys also Hypothesis (f1) stated in Sect. 2,  then this solution is also a classical solution to problem (2.24),  (2.25).

This corollary follows directly from a simple combination of our Theorem 5.1 and a result in Pazy [40, Chapt. 5, §5.7], Theorem 7.1 on p. 168.

Proof of Theorem 5.1

  The “elliptic” differential operator ${\mathcal {A}}(t):\,D_{\mathbb {C}}\subset H_{\mathbb {C}}\rightarrow H_{\mathbb {C}}$ defined in formula (2.5) with the (complex) domain ${\mathcal {D}}({\mathcal {A}}(t)) = D_{\mathbb {C}},$ i.e.,

$$\begin{aligned}&({\mathcal {A}}(t) v)(x) = ({\mathcal {A}}_1(t) v)(x) + ({\mathcal {A}}_2(t) v)(x) \\&\quad \equiv \frac{1}{2}\, [\sigma (t)]^2 \,\frac{\partial ^2 v}{\partial x^2} + \left( q_S(t) - \gamma _S(t) - \frac{1}{2}\, [\sigma (t)]^2 \right) \frac{\partial v}{\partial x} \\&\qquad \hbox {for}\ v\in D_{\mathbb {C}} \hbox { and every time } t\in (0,T) , \end{aligned}$$

is a sum of two commuting differential operators (meaning that their resolvents commute), respectively, ${\mathcal {A}}_1(t) ,\, {\mathcal {A}}_2(t) :\,D_{\mathbb {C}}\subset H_{\mathbb {C}}\rightarrow H_{\mathbb {C}} ,$ where

$$\begin{aligned}&\left\{ \begin{aligned}&({\mathcal {A}}_1(t) v)(x) {\mathop {=}\limits ^{\mathrm{{def}}}}\frac{1}{2}\, [\sigma (t)]^2 \,\frac{\partial ^2 v}{\partial x^2} \quad \hbox { for }\; x\in {\mathbb {R}}^1\\&\quad \hbox { is a diffusion operator in }D_{\mathbb {C}}\quad \hbox {and } \end{aligned} \right. \\&\left\{ \begin{aligned}&({\mathcal {A}}_2(t) v)(x) {\mathop {=}\limits ^{\mathrm{{def}}}}\left( q_S(t) - \gamma _S(t) - \frac{1}{2}\, [\sigma (t)]^2 \right) \frac{\partial v}{\partial x} \quad \hbox {for}\ x\in {\mathbb {R}}^1\\&\quad \hbox { is a convection (or drift) operator in }D_{\mathbb {C}}. \end{aligned} \right. \end{aligned}$$

$\square$

In order to construct the sum ${\mathcal {A}}(t) = {\mathcal {A}}_1(t) + {\mathcal {A}}_2(t)$ with the domain $D_{\mathbb {C}}\subset H_{\mathbb {C}}\rightarrow H_{\mathbb {C}}$ in the Hilbert space $H_{\mathbb {C}} = L^2({\mathbb {R}};{\mathfrak {w}}),$ one may apply a rather general perturbation result from Pazy [40, Chapt. 3, §3.2], Theorem 2.1 on p. 80. Consequently, the (unique) classical solution $v\in C([0,T]\rightarrow H)$ to the homogeneous abstract differential equation

$$\begin{aligned} \genfrac{}{}{}1{\partial v}{\partial t} - {\mathcal {A}}(t) v = 0 \quad \hbox {in}\ H \hbox { for } t\in (0,T); \end{aligned}$$

(5.8)

with the initial condition $v(0) = v_0\in H$ in Eq. (2.25), takes the “$C^0$-semigroup” form

$$\begin{aligned}&v(t) = {\mathfrak {T}}_{1,2}(t,0) v_0 \quad \hbox {in }H \hbox { for } t\in [0,T] , \quad \hbox {where} \end{aligned}$$

(5.9)

$$\begin{aligned}&{\mathfrak {T}}_{1,2}(t,s) = \exp \left( [ {\mathscr {S}}(t) - {\mathscr {S}}(s) ]\, \frac{\partial ^2}{\partial x^2} \right) \,\exp \left( [ {\mathscr {R}}(t) - {\mathscr {R}}(s) ]\, \frac{\partial }{\partial x} \right) \nonumber \\&\quad \hbox {in }H,\hbox { for } 0\le s\le t\in [0,T] , \end{aligned}$$

(5.10)

is the evolutionary family ${\mathfrak {T}}_{1,2} = \{ {\mathfrak {T}}_{1,2}(t,s):\,H\rightarrow H:\,0\le s\le t\le T\}$ of bounded linear operators ${\mathfrak {T}}_{1,2}(t,s):\,H\rightarrow H\ (0\le s\le t\le T)$ that yield the unique classical solution $v(t) = {\mathfrak {T}}_{1,2}(t,s) v_s$ in H to the homogeneous abstract differential equation (5.8) (with $f(t)\equiv 0$) on the interval [s, T] for each $s\in [0,T),$ with the initial condition $v(s) = v_s\in H.$ We recall that the auxiliary functions $t \,\mapsto \, {\mathscr {S}}(t) ,\, {\mathscr {R}}(t):\,[0,T]\rightarrow {\mathbb {R}}$ have been defined in Eq. (5.5). We refer the reader to the general results in Pazy [40, Chapt. 5], Theorem 6.8 in §5.6 on p. 164 and Theorem 7.1 in §5.7 on p. 168. It is easy to see that

$$\begin{aligned}&\left\{ \begin{aligned}&\genfrac{}{}{}1{{\mathrm {d}}}{{\mathrm {d}}t}\, \exp \left( {\mathscr {S}}(t)\, \frac{\partial ^2}{\partial x^2} \right) v_0 = {\mathscr {S}}'(t) \exp \left( {\mathscr {S}}(t)\, \frac{\partial ^2}{\partial x^2} \right) \frac{{\mathrm {d}}^2 v_0}{{\mathrm {d}}x^2}\\&\quad = {\mathscr {S}}'(t)\cdot \exp \left( {\mathscr {S}}(t)\, \frac{\partial ^2}{\partial x^2} \right) v_0'' = \frac{1}{2}\, [\sigma (t)]^2\cdot \exp \left( {\mathscr {S}}(t)\, \frac{\partial ^2}{\partial x^2} \right) v_0''\\&\qquad \hbox {holds for all } v_0\in D_{\mathbb {C}} \hbox { and for every } t\in (0,T] , \quad \hbox {and}\quad \end{aligned} \right. \end{aligned}$$

(5.11)

$$\begin{aligned}&\begin{aligned}&\left[ \exp \left( {\mathscr {R}}(t)\, \frac{\partial }{\partial x} \right) v_0 \right] (x) = v_0( x + {\mathscr {R}}(t) ) \quad \hbox {for all }x\in {\mathbb {R}}^1 \hbox { and }t\in [0,T],\\&\qquad \hbox { together with }\qquad \end{aligned} \end{aligned}$$

(5.12)

$$\begin{aligned}&\left\{ \begin{aligned}&\genfrac{}{}{}1{{\mathrm {d}}}{{\mathrm {d}}t}\, \exp \left( {\mathscr {R}}(t)\, \frac{\partial }{\partial x} \right) v_0 = {\mathscr {R}}'(t) \exp \left( {\mathscr {R}}(t)\, \frac{\partial }{\partial x} \right) \frac{{\mathrm {d}} v_0}{{\mathrm {d}}x}\\&\quad = {\mathscr {R}}'(t)\cdot v_0'( x + {\mathscr {R}}(t) ) = \left( q_S(t) - \gamma _S(t) - \frac{1}{2}\, [\sigma (t)]^2 \right) \cdot v_0'( x + {\mathscr {R}}(t) )\\&\qquad \hbox {for all}\ v_0\in D_{\mathbb {C}} \hbox { and for every } t\in [0,T] . \end{aligned} \right. \end{aligned}$$

(5.13)

Finally, to construct the full linear operator

$$\begin{aligned} {\mathcal {A}}(t) - r_G\,\bullet :\,D_{\mathbb {C}}\subset H_{\mathbb {C}}\rightarrow H_{\mathbb {C}} :\,v \,\longmapsto \, {\mathcal {A}}(t) v - r_G\, v \end{aligned}$$

that appears on the left-hand side of the inhomogeneous abstract differential equation (2.24), let us combine formula (5.9) with Eqs. (5.11) and (5.13) in order to define the evolutionary family ${\mathfrak {T}} = \{ {\mathfrak {T}}(t,s):\,H\rightarrow H:\,0\le s\le t\le T\}$ by the following composition of bounded linear operators

$$\begin{aligned}&{\mathfrak {T}}(t,s) v_s {\mathop {=}\limits ^{\mathrm{{def}}}}\exp ( {}- r_G\, (t-s) )\cdot {\mathfrak {T}}_{1,2}(t,s) v_s \nonumber \\&\quad = \exp ( {}- r_G\, (t-s) )\cdot \exp \left( [ {\mathscr {S}}(t) - {\mathscr {S}}(s) ]\, \frac{\partial ^2}{\partial x^2} \right) \,\exp \left( [ {\mathscr {R}}(t) - {\mathscr {R}}(s) ]\, \frac{\partial }{\partial x} \right) v_s \nonumber \\&\qquad \hbox {in }H \hbox { for }v_s\in H \hbox { and }\, 0\le s\le t\le T . \end{aligned}$$

(5.14)

Applying Eqs. (5.11) and (5.13) to this formula we obtain

$$\begin{aligned}&\genfrac{}{}{}1{{\mathrm {d}}}{{\mathrm {d}}t}\, {\mathfrak {T}}(t,s) v_s = \left( {\mathcal {A}}(t) - r_G\,\bullet \right) {\mathfrak {T}}(t,s) v_s\nonumber \\&\qquad \hbox { in }H_{\mathbb {C}}\hbox { for }v_s\in D_{\mathbb {C}}\hbox { and }\, 0\le s < t\le T . \end{aligned}$$

(5.15)

Let us recall that $f:\,[0,T]\rightarrow H$ is a continuous function. Given any initial value $v_s\in H,$ we may apply [40, Chapt. 5, §5.7], Theorem 7.1 on p. 168, to conclude that the function

$$\begin{aligned} v(t) = {\mathfrak {T}}(t,s) v_s + \int _s^t {\mathfrak {T}}(t,\tau )\, f(\tau ) \,{\mathrm {d}}\tau \quad \hbox {in }H \hbox { for } t\in [s,T] \end{aligned}$$

(5.16)

is the unique mild solution of the inhomogeneous abstract differential equation (2.24) in the time interval [s, T] with the initial value $v(s) = v_s\in H.$

In our proof of Theorem 3.4 above, we have taken advantage of Hypotheses (G2) and (G1’) in order to conclude that the limit function $v\in L^2([0,T]\rightarrow H)$ obtained in formula (3.25) is a strong (classical) solution to problem (2.13), (2.2). Indeed, these two Hypotheses, (G2) and (G1’), guarantee that all iterates (functions) $u_j:\,[0,T]\rightarrow H$; $j=0,1,2,\ldots ,$ in Theorem 3.4 are uniformly Hölder-continuous functions with some Hölder exponent $\vartheta _v\in (0,1)$ and their monotone sequence $u_0\ge u_1\ge \cdots \ge u_{j-1}\ge u_j\ge \cdots \ge {}- u_0$ is uniformly bounded in the Hölder space $C^{\vartheta _v}( [0,T]\rightarrow H).$ Hence, we have $v\in C^{\vartheta _v}([0,T]\rightarrow H),$ as well. We conclude that our function $f(x,t) = G(v(x,t); x,t)$ from Eq. (2.13) obeys Hypothesis (f1). Consequently, from now on let us assume that $f:\,[0,T]\rightarrow H$ is a Hölder-continuous function with the Hölder exponent $\vartheta _f\in (0,1)$; cf. ineq. (2.21) and also ineq. (2.22) related to the linear initial value problem (2.20), (2.2). Then, given any initial value $v_s\in H,$ it follows again from [40, Chapt. 5, §5.7], Theorem 7.1 on p. 168, that the function $v:\,[s,T]\rightarrow H$ defined above by the variation-of-constants formula in (5.16), is even a classical solution of the inhomogeneous abstract differential equation (2.24) in the time interval [s, T] with the initial value $v(s) = v_s\in H.$ From formulas (5.9), (5.12), and (5.16) we deduce that each mapping ${\mathfrak {T}}(t,s):\,H\rightarrow H\ (0\le s < t\le T)$ is an integral operator with the kernel ${\mathfrak {K}}(x,y;t,s) > 0$; see Eq. (5.6).

Next, we will derive an explicit formula for the kernel ${\mathfrak {K}}(x,y;t,s)$ with a help from the standard “diffusion kernel” ${\mathfrak {G}}(x;t)$ defined in Eq. (5.4). First, the standard “diffusion semigroup” ${\mathcal {T}}_1 = \{ {\mathcal {T}}_1(t):\,H\rightarrow H:\,0\le t < \infty \}$ consists of (bounded) integral operators ${\mathcal {T}}_1(t) = \exp \left( t\, \frac{\partial ^2}{\partial x^2} \right) :\,H\rightarrow H\ (0\le t < \infty )$ with the kernel ${\mathfrak {G}}(x-y;t) > 0,$ i.e.,

$$\begin{aligned}&\left[ {\mathcal {T}}_1(t) v_0\right] (x) = \int _{-\infty }^{+\infty } {\mathfrak {G}}(x-y;t)\, v_0(y) \,{\mathrm {d}}y\\&\quad \hbox { defined for all }x\in {\mathbb {R}}^1, 0< t < \infty ,\hbox { and }\, v_0\in H . \end{aligned}$$

We refer to John [29, Chapt. 7, §1(a), pp. 206–213] for greater details. Hence, in formula (5.14) we have the (bounded) integral operator defined for all $0\le s < t\le T$:

$$\begin{aligned}&\left[ {\mathfrak {T}}_1(t,s) v_s\right] (x) {\mathop {=}\limits ^{\mathrm{{def}}}}\left[ \exp \left( [ {\mathscr {S}}(t) - {\mathscr {S}}(s) ]\, \frac{\partial ^2}{\partial x^2} \right) v_s \right] (x)\\&\quad = \int _{-\infty }^{+\infty } {\mathfrak {G}}_1(x-y;t,s)\, v_s(y) \,{\mathrm {d}}y \quad \hbox { for all }x\in {\mathbb {R}}^1 \hbox { and }\, v_s\in H , \end{aligned}$$

with the kernel ${\mathfrak {G}}_1(x-y;t,s) > 0$ given by, for all $x,y\in {\mathbb {R}}^1$ and $0\le s < t\le T$:

$$\begin{aligned} {\mathfrak {G}}_1(x-y;t,s)&= {\mathfrak {G}}\left( x-y ;\, {\mathscr {S}}(t) - {\mathscr {S}}(s) \right) \\ {}&= \frac{1}{ \sqrt{ 4\pi [ {\mathscr {S}}(t) - {\mathscr {S}}(s) ] } }\, \exp \left( {}- \frac{(x-y)^2}{ 4 [ {\mathscr {S}}(t) - {\mathscr {S}}(s) ] } \right) . \end{aligned}$$

Recall that, by Eq. (5.5), we have ${\mathscr {S}}(t) - {\mathscr {S}}(s) = \frac{1}{2}\int _s^t [\sigma (\tau )]^2 \,{\mathrm {d}}\tau > 0$ whenever $0\le s < t\le T.$

Second, the “shift (or translation) group” ${\mathcal {T}}_2 = \{ {\mathcal {T}}_2(t):\,H\rightarrow H:\,t\in {\mathbb {R}}\}$ consists of (bounded) translation operators ${\mathcal {T}}_2(t) = \exp \left( t\, \frac{\partial }{\partial x} \right) :\,H\rightarrow H\ ( -\infty< t < +\infty )$ given by

$$\begin{aligned} \left[ {\mathcal {T}}_2(t) v_0\right] (x) = v_0(x+t) \quad \hbox { for all }x\in {\mathbb {R}}^1, t\in {\mathbb {R}},\hbox { and }\, v_0\in H . \end{aligned}$$

Hence, on the right-hand side in formula (5.14) we have the translation operator

$$\begin{aligned}&\left[ {\mathfrak {T}}_2(t,s) v_s\right] (x) {\mathop {=}\limits ^{\mathrm{{def}}}}\left[ \exp \left( [ {\mathscr {R}}(t) - {\mathscr {R}}(s) ]\, \frac{\partial }{\partial x} \right) v_s \right] (x) = v_s( x + {\mathscr {R}}(t) - {\mathscr {R}}(s) )\\&\quad \hbox {for all }x\in {\mathbb {R}}^1, 0\le s\le t\le T, \hbox { and }\, v_s\in H . \end{aligned}$$

Consequently, the (bounded linear) composition operator on H,  defined by formula (5.10), i.e.,

$$\begin{aligned}&{\mathfrak {T}}_{1,2}(t,s) {\mathop {=}\limits ^{\mathrm{{def}}}}{\mathfrak {T}}_1(t,s)\, {\mathfrak {T}}_2(t,s) \equiv {\mathfrak {T}}_1(t,s)\circ {\mathfrak {T}}_2(t,s)\\&\quad = \exp \left( [ {\mathscr {S}}(t) - {\mathscr {S}}(s) ]\, \frac{\partial ^2}{\partial x^2} \right) \,\exp \left( [ {\mathscr {R}}(t) - {\mathscr {R}}(s) ]\, \frac{\partial }{\partial x} \right) :\,H\rightarrow H , \end{aligned}$$

for $0\le s < t\le T,$ is an integral operator defined for all $0\le s < t\le T$:

$$\begin{aligned} \left[ {\mathfrak {T}}_{1,2}(t,s) v_s\right] (x)&{} = \int _{-\infty }^{+\infty } {\mathfrak {G}}_{1,2}(x,y;t,s)\, v_s(y) \,{\mathrm {d}}y \quad \hbox { for all }x\in {\mathbb {R}}^1 \hbox { and }\, v_s\in H , \end{aligned}$$

with the kernel ${\mathfrak {G}}_{1,2}(x,y;t,s) > 0$ given by, for all $x,y\in {\mathbb {R}}^1$ and $0\le s < t\le T$:

$$\begin{aligned} {\mathfrak {G}}_{1,2}(x,y;t,s)&\equiv {\mathfrak {G}}_{1,2}(x-y;t,s) = {\mathfrak {G}}_1 \bigl ( x + {\mathscr {R}}(t) - {\mathscr {R}}(s) - y ;\, t,s \bigr )\\&= {\mathfrak {G}} \bigl ( x-y + {\mathscr {R}}(t) - {\mathscr {R}}(s) ;\, {\mathscr {S}}(t) - {\mathscr {S}}(s) \bigr )\\&= \frac{1}{ \sqrt{ 4\pi [ {\mathscr {S}}(t) - {\mathscr {S}}(s) ] } }\, \exp \left( {}- \frac{( x-y + {\mathscr {R}}(t) - {\mathscr {R}}(s) )^2}{ 4 [ {\mathscr {S}}(t) - {\mathscr {S}}(s) ] } \right) . \end{aligned}$$

Finally, comparing ${\mathcal {A}}(t):\,D_{\mathbb {C}}\subset H_{\mathbb {C}}\rightarrow H_{\mathbb {C}}$ with the full linear operator

$$\begin{aligned} {\mathcal {A}}(t) - r_G\,\bullet :\,D_{\mathbb {C}}\subset H_{\mathbb {C}}\rightarrow H_{\mathbb {C}} :\,v \,\longmapsto \, {\mathcal {A}}(t) v - r_G\, v , \end{aligned}$$

we arrive at ${\mathfrak {T}}(t,s) = \exp ( {}- r_G\, (t-s) )\cdot {\mathfrak {T}}_{1,2}(t,s)$ in H for all $0\le s < t\le T.$ Hence, Eq. (5.14) is valid and the desired formulas (5.6) and (5.7) follow.

Our proof of Theorem 5.1 is finished. $\square$

Remark 5.3

(An alternative proof of Theorem 5.1). In order to derive formula (5.7) for the kernel ${\mathfrak {K}}(x,y;t,s) > 0$ of the integral operator ${\mathfrak {T}}(t,s):\,H\rightarrow H$ in Eq. (5.6), one can replace the calculations in our proof of Theorem 5.1 above by applying first the Fourier transformation to the homogeneous parabolic problem (5.2), (5.3) to calculate the Fourier transform of the kernel ${\mathfrak {K}}(x,y;t,s)\equiv {\mathfrak {K}}(x-y;t,s)$ with respect to the variable (difference) $x-y\in {\mathbb {R}}^1,$ followed by a standard application of the inverse Fourier transformation to obtain the desired kernel ${\mathfrak {K}}(x,y;t,s)$ for all $x,y\in {\mathbb {R}}^1$ and $0\le s< t < \infty .$ This procedure is analogous with the derivation of the “diffusion kernel” ${\mathfrak {G}}(x;t)$ (defined here in Eq. (5.4)) by Fourier transformation in John [29, Chapt. 7, §1(a), pp. 208–209]; see Eq. (1.10d) on p. 209. $\square$

Our explicit formula (5.1) verified in Corollary 5.2 (to Theorem 5.1) renders a precise analytic way of calculating the exact solution to problem (2.24), (2.25), which is the abstract form of the Cauchy problem (2.20), (2.2). In practical, real world applications to problems in Mathematical Finance this solution needs to be approximated by suitable Numerical Methods; typically with reasonable speed and within reasonable precision. Whereas the , modelled by the nonlinear “reaction” function $G(\,\cdot \, ; x,t):\,v\mapsto G(v; x,t):\,{\mathbb {R}}\rightarrow {\mathbb {R}},$ defined in Eq. (2.8), are computed in a most standard way through the inhomogeneity $f(x,t) = G(v(x,t); x,t)$ from Eq. (2.13), with the (nonlinear) substitution operator $G(\,\cdot \, ; x,t):\,{\mathbb {R}}\rightarrow {\mathbb {R}},$ the numerical computation of the (unique) solution $v:\,{\mathbb {R}}^1\times (0,T)\rightarrow {\mathbb {R}}$ to the initial value Cauchy problem (2.20), (2.2) is much more complicated:

Remark 5.4

(i) Eq. (2.20) being a linear “diffusion equation”, the Monte Carlo method is a very natural way for approximating the precise analytic solution to problem (2.20), (2.2) (i.e., (2.24), (2.25)) by numerical simulations produced by this method. We refer to the works by Arregui et al. [3, Sects. 3–4, pp. 18–23], Baustian et al. [6, §3.6, p. 52] and to Plotnikov [41, Sect. 1, pp. 121–125] for greater details. The first two references, [3, 6], treat exactly the problem of “Option pricing under some Value Adjustment” () in Mathematical Finance; under “Credit Value Adjustment” (), for instance. The third one, [41], treats linear integral equations of type (5.16) which envolve the evolutionary family of bounded linear integral operators ${\mathfrak {T}}(t,s):\,H\rightarrow H\ (0\le s < t\le T)$ on the Hilbert space H. In fact, a numerical method for solving the full, nonlinear integral equation (5.1) for the unknown function $v\in C([0,T]\rightarrow H)$ substituted into the (subsequently unknown) nonlinearity

$$\begin{aligned} f(\tau )\equiv f( \,\cdot \, ,\tau ):\,{\mathbb {R}}^1\rightarrow {\mathbb {R}}:\,x \,\mapsto \, f(x,\tau ) = G(v(x,\tau ); x,\tau ) \quad \hbox {with } f(\tau )\in H \end{aligned}$$

for every $\tau \in [0,T],$ is provided in [41]. This work is based in an iteration method for a nonlinear integral equation similar to ours, cf. [41, Eq. (1.1), p. 121].

(ii) In contrast with the probabilistic Monte Carlo methods for solving the Cauchy problem (2.20), (2.2), analytic methods based on a finite difference (or finite element) scheme provide a highly competitive alternative to Monte Carlo in a series of works, such as the monograph by Achdou and Pironneau [1] and the articles by Arregui et al. [3, §3.4, pp. 20–21], [4, Sect. 4, pp. 734–737], Baustian et al. [6, §3.5, pp. 50–51], Koleva [30, Sect. 3, pp. 367–368], and Koleva and Vulkov [31, Sects. 3–4, pp. 510–515].

(iii) Last but not least, a “hybrid” algorithm mixing Monte Carlo with finite difference/element methods in quest for optimization on both, precision and speed, is presented in Loeper and Pironneau [37]. $\square$