# Valuation of the Convertible Bonds under Penalty TF model using Finite
Element Method

Rakhymzhan Kazbek  Nazarbayev University, Department of Mathematics, School of
Sciences and Humanities, 53 Kabanbay Batyr Ave, Nur-Sultan 010000, Kazakhstan
Yogi Erlangga  Zayed University, Department of Mathematics, College of Natural
and Health Sciences, Abu Dhabi Campus, P.O. Box 144534, United Arab Emirates
Corresponding author  Yerlan Amanbek  Nazarbayev University, Department of
Mathematics, School of Sciences and Humanities, 53 Kabanbay Batyr Ave, Nur-
Sultan 010000, Kazakhstan  Dongming Wei  Nazarbayev University, Department of
Mathematics, School of Sciences and Humanities, 53 Kabanbay Batyr Ave, Nur-
Sultan 010000, Kazakhstan

###### Abstract

In this paper, the TF system of two-coupled Black-Scholes equations for
pricing the convertible bonds is solved numerically by using the P1 and P2
finite elements with the inequality constraints approximated by the penalty
method. The corresponding finite element ODE system is numerically solved by
using a modified Crank-Nicolson scheme, in which the non-linear system is
solved at each time step by the Newton-Raphson method for non-smooth
functions. Moreover, the corresponding Greeks are also calculated by taking
advantage of the P1-P2 finite element approximation functions. Numerical
solutions by the finite element method compare favorably with the solutions by
the finite difference method in literature.

Keywords— TF model, Pricing convertible bonds, Finite element method, Finite
difference method, Financial derivatives, Greeks.

##  1 Introduction

Convertible bond (CB) is a type of financial derivatives used in financial
risk management and in trading by investors and the issuers [18, 27, 26]. It
has become a popular choice of corporations as a stable investment vehicle.
Pricing financial derivatives for CBs is however a complicated problem. Tree
methods (e.g., in [23]), Monte-Carlo simulation (e.g., [4]), and partial
differential equations are among the widely used techniques in the literature.

One of the pioneering work using PDE modelling for CB pricing is due to
Ignersoll [19], who developed a method for the determination of the optimal
conversion and call policies for convertible securities using Black-Scholes
[8] methodology. Brennan and Schwartz [9] extended Ignersoll’s model to
incorporate callability and dividend payments. They later included stochastic
interest rate, resulting in a two-factor model for CB pricing with callability
and conversion strategies [10]. In the late 1990s, Tsiveriotis and Fernandes
[28] proposed an innovative two-factor model for pricing CBs under
callability, puttability and conversion provisions with hypothetical cash-only
convertible bond (COCB) part, known as the TF model. Another model for pricing
CBs with credit risk and default strategy was developed by Ayache et al. [6]
based on a system of triple partial differential equations, called the AFV
model.

Because of their complexity, the TF and AFV models have to be solved
numerically. The review paper of Al Saedi and Tularam [1] provides a recent
account on methodologies for option pricing under the Black-Scholes equation,
such as the R3C scheme [5], cubic spline wavelets and multiwavelet bases
methods [11], finite difference methods (FDM) [14], finite-volume methods
(FVM) [22], and finite-element methods (FEM) [17, 7, 13, 20]. Due to their
simplicity, FDM have become a popular numerical method in computational
finance. In computational bond pricing, the authors of [9, 10, 28, 6] use FDM
to solve their proposed model.

In this paper, we discuss and expanded the existing literature on numerical
solutions of the TF model for CB pricing with credit risk and without
dividends by implementing finite-element methods. The TF model for convertible
bond pricing is based on the system of partial differential equations (PDEs)

$\displaystyle\displaystyle\frac{\partial U}{\partial
t}+\frac{\sigma^{2}}{2}S^{2}\frac{\partial^{2}U}{\partial
S^{2}}+r_{g}S\frac{\partial U}{\partial S}-r(U-V)-(r+r_{c})V=0,$ (1)
$\displaystyle\displaystyle\frac{\partial V}{\partial
t}+\frac{\sigma^{2}}{2}S^{2}\frac{\partial^{2}V}{\partial
S^{2}}+r_{g}S\frac{\partial V}{\partial S}-(r+r_{c})V=0,$ (2)

for the time $t\in(0,T)$ and the underlying stock price $S\in(0,\infty)$, with
$U$ being the value of the CB, $V$ the value of hypothetical COCB, $r$ the
risk-free rate, $r_{g}$ the growth rate, which can be counted as risk-free
rate $r$ (see [18]), $r_{c}$ the credit spread reflecting payoff default risk,
$\sigma$ be volatility, and $T$ the maturity time.

The terminal condition at the maturity time $T$ means that once the CB is
expired, no one can call or put it back. Holder of the CB will get as much as
possible depending on the conversion ratio $k$ and stock price. There is
however a minimum based on the face value $F$ and the coupon payment $K$,
yielding the terminal condition:

$\displaystyle U(S,T)=\begin{cases}F+K,&\text{if }F+K\geq kS,\\\
kS,&\text{otherwise},\end{cases}$ (3)

and

$\displaystyle V(S,T)=\begin{cases}F+K,&\text{if }F+K\geq kS,\\\
0,&\text{otherwise}.\end{cases}$ (4)

Throughout its lifetime, however, the CB can be converted to the underlying
stock at the value $kS$, and issuer should pay the principal value $F+K$ to
the holder, if the issuer has not converted until maturity time. These rights
lead to three conditions that constrain the CB price:

  1. 1.

Upside constraint due to conversion of bonds: for $t\in[0,T]$,

$\displaystyle U(S,t)\geq kS,$ (5) $\displaystyle V(S,t)=0\text{ if
}U(S,t)\leq kS;$ (6)

  2. 2.

Upside constraint due to callability with the call price $B_{call}$: for
$t\in[T_{call},T]$, with $T_{call}$ the earliest time the bond issuer is
allowed to call back the bond,

$\displaystyle U(S,t)\leq\max(B_{call}(t),kS),$ (7) $\displaystyle\quad\quad
V(S,t)=0\text{ if }U(S,t)\geq B_{call}(t);$ (8)

  3. 3.

Downside constraint due to putability with the put price $B_{put}$: for
$t\in[T_{put},T]$, with $T_{put}$ the earliest time the investor is allowed to
put the bond back,

$\displaystyle U(S,t)\geq B_{put}(t),$ (9) $\displaystyle
V(S,t)=B_{put}(t)\text{ if }U(S,t)\leq B_{put}(t).$ (10)

Following [6], the call and put price in the callability and puttability
constraints include the effect of future’s coupon payment and the underlying
interest as follows. Let $\mathcal{T}_{\text{coupon}}=\\{t_{i}\\}$, the set of
the coupon payment time, with $0<t_{i-1}<t_{i}\leq T$. Then

$\displaystyle B_{put,call}(t)=B_{put,call}^{cl}+AccI(t),$ (11)

where

$\displaystyle AccI(t)=K_{i}\frac{t-t_{i-1}}{t_{i}-t_{i-1}},$ (12)

the accrued interest at any time $t$ between the time of the last coupon
payment $t_{i-1}$ and the time of the next (pending) coupon payment $t_{i}$.
Note that the constraints (5) and (9) can be combined to
$U(S,t)\geq\max\\{B_{put},kS\\}$. In this way, the constraints for $U$ can be
rewritten as

$\displaystyle U(S,t)$ $\displaystyle\leq\max(B_{call}(t),kS),$ (13)
$\displaystyle U(S,t)$ $\displaystyle\geq\max(B_{put}(t),kS),$ (14)

with

$\displaystyle B_{call}(t)$
$\displaystyle=\begin{cases}B_{call}^{cl}+AccI(t),&\text{if
}t\in[T_{call},T],\\\ +\infty,&\text{otherwise},\end{cases}$ (15)
$\displaystyle B_{put}(t)$
$\displaystyle=\begin{cases}B_{put}^{cl}+AccI(t),&\text{if
}t\in[T_{put},T],\\\ 0,&\text{otherwise}.\end{cases}$ (16)

Two boundary conditions need to be supplemented to the PDEs (1) and (2). At
$S=0$, the PDEs are reduced to

$\displaystyle\displaystyle\frac{\partial U(0,t)}{\partial t}$
$\displaystyle=rU(0,t)+r_{c}V(0,t),$ (17)

and

$\displaystyle\displaystyle\frac{\partial V(0,t)}{\partial t}$
$\displaystyle=(r+r_{c})V(0,t),$ (18)

with the terminal conditions $U(0,T)=V(0,T)=F+K$ (see Eqs. (3) and (4)). The
other boundary condition is associated with the situation when the the stock
price $S$ increases unboundedly, under which the CB is converted into stock.
Therefore,

$\displaystyle\lim_{S\to\infty}\begin{cases}U(S,t)=kS,\\\
V(S,t)=0.\end{cases}$ (19)

After spatial discretization using, for instance, FDM or FEM, the initial-
boundary value problems (IBVP) with constraints given by (1)–(19) can be
solved by using a time-integration method, such as Crank-Nicolson, in
combination with projected SSOR (PSSOR) method to tackle the nonlinearity. In
our paper, however, we shall consider a formulation of the above-stated IBVP
in a penalty PDE to explicitly include some constraints in the PDEs, and apply
FEM and Crank-Nicolson method on the resulting penalty PDE.

The paper is organized as follows. In Section 2, transformation of the TF
model with penalty in a suitable form for numerical method is presented to
implement our numerical schemes. Finite element formulation for the penalty TF
model is established in Section 3. We detail the time-integration method for
solving the resulting differential algebraic equation in Section 4. Section 5
presents numerical results and discussions on the convergence of the FEM
method. We finish up the paper by drawing some conclusion and remarks in
Section 6.

##  2 The transformed TF model

The standard procedure for solving the Black-Scholes-type PDEs requires
transformation of the terminal-boundary value problem to an initial-boundary
value problem, on which many numerical methods can be devised. Let
$\displaystyle\tau=T-t$ and $\displaystyle
x=\ln\left(S/S_{\text{int}}\right)$, where $S_{\text{int}}$ is the stock price
at the initial time $t=0$. With this change of variables, it can be shown that
the PDEs (1)–(2) are transformed to

$\displaystyle\displaystyle\frac{\partial U}{\partial\tau}$
$\displaystyle=\frac{\sigma^{2}}{2}\frac{\partial^{2}U}{\partial
x^{2}}+\left(r-\frac{\sigma^{2}}{2}\right)\frac{\partial U}{\partial
x}-r(U-V)-(r+r_{c})V,$ (20) $\displaystyle\displaystyle\frac{\partial
V}{\partial\tau}$
$\displaystyle=\frac{\sigma^{2}}{2}\frac{\partial^{2}V}{\partial
x^{2}}+\left(r-\frac{\sigma^{2}}{2}\right)\frac{\partial V}{\partial
x}-(r+r_{c})V,$ (21)

with $(\tau,x)\in(0,T)\times(-\infty,\infty)$. The terminal conditions at
$t=T$ now become initial conditions at $\tau=0$, given by

$\displaystyle U(x,0)=\begin{cases}F+K,&\text{if }F+K\geq
kS_{\text{int}}e^{x},\\\ kS_{\text{int}}e^{x},&\text{otherwise},\end{cases}$
(22)

and

$\displaystyle V(x,0)=\begin{cases}F+K,&\text{if }F+K\geq
kS_{\text{int}}e^{x},\\\ 0,&\text{otherwise}.\end{cases}$ (23)

Furthermore, the constraints for $U$ are transformed into

$\displaystyle U(x,\tau)$
$\displaystyle\geq\max\\{B_{put}(\tau),kS_{\text{int}}e^{x}\\}=:U_{put}^{\star}(\tau),$
(24) $\displaystyle U(x,\tau)$
$\displaystyle\leq\max\\{B_{call}(\tau),kS_{\text{int}}e^{x}\\}=:U_{call}^{\star}(\tau),$
(25)

where

$\displaystyle B_{call}(\tau)$
$\displaystyle=\begin{cases}B_{call}^{cl}+AccI(\tau),&\text{if
}t\in[0,\tau_{call}],\\\ +\infty,&\text{otherwise},\end{cases}$ (26)

and

$\displaystyle B_{put}(\tau)$
$\displaystyle=\begin{cases}B_{put}^{cl}+AccI(\tau),&\text{if
}t\in[0,\tau_{put}],\\\ 0,&\text{otherwise},\end{cases}$ (27)

with $\tau_{put,call}=T-T_{put,call}$ and, for
$\tau_{i}\in\mathcal{T}_{\tau,\text{coupon}}=\mathcal{T}_{\text{coupon}}$,

$\displaystyle\displaystyle
B_{put,call}(\tau)=B^{cl}_{put,call}+K_{i}\frac{\tau-\tau_{i-1}}{\tau_{i}-\tau_{i-1}}.$
(28)

The constraints for $V$ now read

  1. 1.

for $\tau\in[0,T]$,

$\displaystyle V(x,\tau)=0,\text{ if }U(\tau)\leq kS_{\text{int}}e^{x},$ (29)

  2. 2.

for $\tau\in[0,\tau_{call}]$,

$\displaystyle V(x,\tau)=0,\text{ if }U(\tau)\geq B_{call}(\tau)$ (30)

  3. 3.

for $\tau\in[0,\tau_{put}]$,

$\displaystyle V(x,\tau)=B_{put},\text{ if }U(\tau)\leq B_{put}(\tau).$ (31)

For the boundary conditions, we note here that $x$ is not defined at $S=0$. As
in the actual numerical computation we set $S$ as close as possible to $0$, we
assume that (17) also holds at the proximity of $S=0$, which corresponds to
$x_{\min}\to-\infty$ in the $x$-space. This results in the boundary conditions
at $x_{\min}$

$\displaystyle\frac{\partial U(x_{\min},\tau)}{\partial\tau}$
$\displaystyle=-rU(x_{\min},\tau)-r_{c}V(x_{\min},\tau),$ (32)
$\displaystyle\frac{\partial V(x_{\min},\tau)}{\partial t}$
$\displaystyle=-(r+r_{c})V(x_{\min},\tau).$ (33)

Transformation of the boundary conditions at $S\to+\infty$ is straightforward:
at $x_{\max}\to+\infty$,

$\displaystyle U(x_{\max},\tau)$ $\displaystyle=kS_{int}e^{x_{\max}},$ (34)
$\displaystyle\displaystyle V(x_{\max},\tau)$ $\displaystyle=0.$ (35)

As stated in Section 1, our focus in this paper is on the the penalty TF
model. We therefore need to reformulate the model into a PDE model with
penalty terms associated with some constraints. In particular, as in practice
we are mainly concerned with the CB price $U$, not $V$, we shall reformulate
the CB PDE (20) with the associated constraints (24) and (25) into a penalty
PDE. To this end, note that the linear complementarity problem (LCP) for (20)
with constraints (24) and (25) is given by [6]

$\left(\begin{array}[]{l}\mathcal{L}U-r_{c}V=0\\\ U\geq\max(B_{p},\kappa
S_{\text{int}}e^{x})\\\ U\leq\max(B_{c},\kappa
S_{\text{int}}e^{x})\end{array}\right)\vee\left(\begin{array}[]{l}\mathcal{L}U-r_{c}V\leq
0\\\ U=\max(B_{p},\kappa S_{\text{int}}e^{x})\\\ U\leq\max(B_{c},\kappa
S_{\text{int}}e^{x})\end{array}\right)\vee\left(\begin{array}[]{l}\mathcal{L}U-r_{c}V\geq
0\\\ U\geq\max(B_{p},\kappa S_{\text{int}}e^{x})\\\ U=\max(B_{c},\kappa
S_{\text{int}}e^{x})\end{array}\right),$ (36)

where
$\mathcal{L}=-\displaystyle\frac{\partial}{\partial\tau}+\frac{\sigma^{2}}{2}\frac{\partial^{2}}{\partial
x^{2}}+\left(r-\frac{\sigma^{2}}{2}\right)\frac{\partial}{\partial x}-r$.

The penalty PDE for the bond valuation can be constructed from the LCP (36):

$\displaystyle\frac{\partial
U}{\partial\tau}=\frac{\sigma^{2}}{2}\frac{\partial^{2}U}{\partial
x^{2}}+\left(r-\frac{\sigma^{2}}{2}\right)\frac{\partial U}{\partial x}-rU-
r_{c}V+\rho\max(U-U_{call}^{\star},0)+\rho\max(U_{put}^{\star}-U,0),$

where $\rho>0$ is the penalty parameter, which typically is set very large. By
rewriting $\max(U-U_{call}^{\star},0)=\alpha_{call}(U-U_{call}^{\star})$ and
$\max(U_{put}^{\star}-U,0)=\alpha_{put}(U-U_{put}^{\star})$, with

$\displaystyle\alpha_{call}=\begin{cases}1,&\text{if }U-U_{call}^{\star}\geq
0,\\\
0,&\text{otherwise},\end{cases}\quad\text{and}\quad\alpha_{put}=\begin{cases}1,&\text{if
}U_{put}^{\star}-U\geq 0,\\\ 0,&\text{otherwise},\end{cases}$ (37)

the CB PDE can be reformulated into

$\displaystyle\frac{\partial
U}{\partial\tau}=\frac{\sigma^{2}}{2}\frac{\partial^{2}U}{\partial
x^{2}}+(r-\frac{\sigma^{2}}{2})\frac{\partial U}{\partial x}-rU-
r_{c}V+\rho\alpha_{call}(U-U_{call}^{\star})+\rho\alpha_{put}(U_{put}^{\star}-U).$
(38)

##  3 Finite element method

In this section, we construct finite-element methods to approximately
discretize the penalty CB PDE (38) and the COCB PDE (21) in the transformed TF
model. To apply finite element spatial discretization to the PDEs, the domain
$(0,\infty)$ is approximated by the bounded domain
$\Omega=(x_{\min},x_{\max})$ in the following:

Suppose that the solution of the bond PDEs are functions of the following
class:

$U(x,\tau),V(x,\tau)\in\mathcal{H}^{1}=\left\\{f:\Omega\to\mathbb{R}\mid
f,\frac{\partial f}{\partial x}\in L_{2}(\Omega)\text{ and $f$ satisfies
boundary conditions on }\partial\Omega\right\\}.$

Consider two test functions $w,z\in\mathcal{H}^{1}_{0}$, where
$\mathcal{H}^{1}_{0}=\left\\{f:\Omega\to\mathbb{R}\mid f,\frac{\partial
f}{\partial x}\in L_{2}(\Omega)\text{ and }f|_{\partial\Omega}=0\right\\}$.
The weak formulation of the penalty PDE (38) and (21) reads

$\displaystyle\int\displaylimits_{\Omega}w\frac{\partial U}{\partial\tau}$
$\displaystyle=\frac{\sigma^{2}}{2}\int\displaylimits_{\Omega}w\frac{\partial^{2}U}{\partial
x^{2}}+\left(r-\frac{\sigma^{2}}{2}\right)\int\displaylimits_{\Omega}w\frac{\partial
U}{\partial x}-r\int\displaylimits_{\Omega}wU-
r_{c}\int\displaylimits_{\Omega}wV+\rho\int\displaylimits_{\Omega}\alpha_{call}w(U-U^{\star}_{call})+\rho\int\displaylimits_{\Omega}\alpha_{put}w(U^{\star}_{put}-U),$
$\displaystyle\int\displaylimits_{\Omega}z\frac{\partial V}{\partial\tau}$
$\displaystyle=\frac{\sigma^{2}}{2}\int\displaylimits_{\Omega}z\frac{\partial^{2}V}{\partial
x^{2}}+\left(r-\frac{\sigma^{2}}{2}\right)\int\displaylimits_{\Omega}z\frac{\partial
V}{\partial x}-(r+r_{c})\int\displaylimits_{\Omega}zV,$

where the integration is carried out along the $x$-direction ($dx$ is not
indicated to save space). Integration by parts and applying the vanishing
property of the function $w$ and $z$ at the boundary results in the weak
formulation

$\displaystyle\frac{\partial}{\partial\tau}\int\displaylimits_{\Omega}wU$
$\displaystyle=-\frac{\sigma^{2}}{2}\int\displaylimits_{\Omega}\frac{\partial
w}{\partial x}\frac{\partial U}{\partial
x}-\left(r-\frac{\sigma^{2}}{2}\right)\int\displaylimits_{\Omega}\frac{\partial
w}{\partial x}U-r\int\displaylimits_{\Omega}wU-
r_{c}\int\displaylimits_{\Omega}wV+\rho\int\displaylimits_{\Omega}\alpha_{call}w(U-U^{\star}_{call})+\rho\int\displaylimits_{\Omega}\alpha_{put}w(U^{\star}_{put}-U),$
$\displaystyle\frac{\partial}{\partial\tau}\int\displaylimits_{\Omega}zV$
$\displaystyle=-\frac{\sigma^{2}}{2}\int\displaylimits_{\Omega}\frac{\partial
z}{\partial x}\frac{\partial V}{\partial
x}-\left(r-\frac{\sigma^{2}}{2}\right)\int\displaylimits_{\Omega}\frac{\partial
z}{\partial x}V-(r+r_{c})\int\displaylimits_{\Omega}zV.$

To build our finite-element approximation, we consider the finite dimensional
subspace $S^{h}_{0}\subset\mathcal{H}^{1}_{0}$, spanned by the basis
$\\{\psi_{1},\psi_{2},\dots,\psi_{n}\\}$. The finite-element approximation to
the solution $U$ is the function

$U_{h}=\sum_{i=1}^{n}u_{i}\psi_{i}+\sum_{i\in\mathcal{I}_{\partial}}u_{i}\psi_{i}\simeq
U,\quad u_{i}\in\mathbb{R},$

where $\psi_{i\in\mathcal{I}_{\partial}}$ are additional functions needed to
interpolate the given solutions at the boundaries. A similar form of
approximation to $V$ is devised, namely

$V_{h}=\sum_{i=1}^{n}v_{i}\psi_{i}+\sum_{i\in\mathcal{I}_{\partial}}v_{i}\psi_{i}\simeq
V,\quad v_{i}\in\mathbb{R}.$

The use of the above approximations results in the weak formulation in the
finite-dimensional space:

$\displaystyle\frac{\partial}{\partial\tau}\left(\sum_{i=1}^{n}u_{i}\int\displaylimits_{\Omega}w\psi_{i}+\sum_{i\in\mathcal{I}_{\partial}}u_{i}\int\displaylimits_{\Omega}w\psi_{i}\right)$
$\displaystyle=-\frac{\sigma^{2}}{2}\left(\sum_{i=1}^{n}u_{i}\int\displaylimits_{\Omega}\frac{\partial
w}{\partial x}\frac{\partial\psi_{i}}{\partial
x}+\sum_{i\in\mathcal{I}_{\partial}}u_{i}\int\displaylimits_{\Omega}\frac{\partial
w}{\partial x}\frac{\partial\psi_{i}}{\partial x}\right)$ (39)
$\displaystyle-\left(r-\frac{\sigma^{2}}{2}\right)\left(\sum_{i=1}^{n}u_{i}\int\displaylimits_{\Omega}\frac{\partial
w}{\partial
x}\psi_{i}+\sum_{i\in\mathcal{I}_{\partial}}^{n}u_{i}\int\displaylimits_{\Omega}\frac{\partial
w}{\partial x}\psi_{i}\right)$
$\displaystyle-r\left(\sum_{i=1}^{n}u_{i}\int\displaylimits_{\Omega}w\psi_{i}+\sum_{i\in\mathcal{I}_{\partial}}u_{i}\int\displaylimits_{\Omega}w\psi_{i}\right)-r_{c}\left(\sum_{i=1}^{n}v_{i}\int\displaylimits_{\Omega}w\psi_{i}+\sum_{i\in\mathcal{I}_{\partial}}v_{i}\int\displaylimits_{\Omega}w\psi_{i}\right)$
$\displaystyle+\mathcal{P}_{call}+\mathcal{P}_{put},$

where
$\displaystyle\mathcal{P}_{call}=\rho\int\displaylimits_{\Omega}\alpha_{call}w\left(U-U_{call}^{\star}\right)dx$
and
$\displaystyle\mathcal{P}_{put}=\rho\int\displaylimits_{\Omega}\alpha_{put}w\left(U_{put}^{\star}-U\right)dx$,
and

$\displaystyle\frac{\partial}{\partial\tau}\left(\sum_{i=1}^{n}v_{i}\int\displaylimits_{\Omega}z\psi_{i}+\sum_{i\in\mathcal{I}_{\partial}}v_{i}\int\displaylimits_{\Omega}z\psi_{i}\right)$
$\displaystyle=-\frac{\sigma^{2}}{2}\left(\sum_{i=1}^{n}v_{i}\int\displaylimits_{\Omega}\frac{\partial
z}{\partial x}\frac{\partial\psi_{i}}{\partial
x}+\sum_{i\in\mathcal{I}_{\partial}}v_{i}\int\displaylimits_{\Omega}\frac{\partial
z}{\partial x}\frac{\partial\psi_{i}}{\partial x}\right)$ (40)
$\displaystyle-\left(r-\frac{\sigma^{2}}{2}\right)\left(\sum_{i=1}^{n}v_{i}\int\displaylimits_{\Omega}\frac{\partial
z}{\partial
x}\psi_{i}+\sum_{i\in\mathcal{I}_{\partial}}^{n}v_{i}\int\displaylimits_{\Omega}\frac{\partial
z}{\partial x}\psi_{i}\right)$
$\displaystyle-(r+r_{c})\left(\sum_{i=1}^{n}v_{i}\int\displaylimits_{\Omega}z\psi_{i}+\sum_{i\in\mathcal{I}_{\partial}}v_{i}\int\displaylimits_{\Omega}z\psi_{i}\right).$

In the Galerkin method [21], the test function $w$ and $z$ are chosen to
coincide with the basis function $\psi_{i}$. Imposing this condition for
$w,z=\psi_{j}$, $j=1,\dots,n$ results in the system of equations

$\displaystyle\frac{\partial}{\partial\tau}\left(\sum_{i=1}^{n}u_{i}\int\displaylimits_{\Omega}\psi_{j}\psi_{i}+\sum_{i\in\mathcal{I}_{\partial}}u_{i}\int\displaylimits_{\Omega}\psi_{j}\psi_{i}\right)$
$\displaystyle=-\frac{\sigma^{2}}{2}\left(\sum_{i=1}^{n}u_{i}\int\displaylimits_{\Omega}\frac{\partial\psi_{j}}{\partial
x}\frac{\partial\psi_{i}}{\partial
x}+\sum_{i\in\mathcal{I}_{\partial}}u_{i}\int\displaylimits_{\Omega}\frac{\partial\psi_{j}}{\partial
x}\frac{\partial\psi_{i}}{\partial x}\right)$ (41)
$\displaystyle-\left(r-\frac{\sigma^{2}}{2}\right)\left(\sum_{i=1}^{n}u_{i}\int\displaylimits_{\Omega}\frac{\partial\psi_{j}}{\partial
x}\psi_{i}+\sum_{i\in\mathcal{I}_{\partial}}^{n}u_{i}\int\displaylimits_{\Omega}\frac{\partial\psi_{j}}{\partial
x}\psi_{i}\right)$
$\displaystyle-r\left(\sum_{i=1}^{n}u_{i}\int\displaylimits_{\Omega}\psi_{j}\psi_{i}+\sum_{i\in\mathcal{I}_{\partial}}u_{i}\int\displaylimits_{\Omega}\psi_{j}\psi_{i}\right)-r_{c}\left(\sum_{i=1}^{n}v_{i}\int\displaylimits_{\Omega}\psi_{j}\psi_{i}+\sum_{i\in\mathcal{I}_{\partial}}v_{i}\int\displaylimits_{\Omega}\psi_{j}\psi_{i}\right)$
$\displaystyle+\mathcal{P}_{call,j}+\mathcal{P}_{put,j},$

where $\mathcal{P}_{call,j}$ is $\mathcal{P}_{call}$ with $w$ be replaced by
$\psi_{j}$ and similarly for $\mathcal{P}_{put,j}$, and

$\displaystyle\frac{\partial}{\partial\tau}\left(\sum_{i=1}^{n}v_{i}\int\displaylimits_{\Omega}\psi_{j}\psi_{i}+\sum_{i\in\mathcal{I}_{\partial}}v_{i}\int\displaylimits_{\Omega}\psi_{j}\psi_{i}\right)$
$\displaystyle=-\frac{\sigma^{2}}{2}\left(\sum_{i=1}^{n}v_{i}\int\displaylimits_{\Omega}\frac{\partial\psi_{j}}{\partial
x}\frac{\partial\psi_{i}}{\partial
x}+\sum_{i\in\mathcal{I}_{\partial}}v_{i}\int\displaylimits_{\Omega}\frac{\partial\psi_{j}}{\partial
x}\frac{\partial\psi_{i}}{\partial x}\right)$ (42)
$\displaystyle-\left(r-\frac{\sigma^{2}}{2}\right)\left(\sum_{i=1}^{n}v_{i}\int\displaylimits_{\Omega}\frac{\partial\psi_{j}}{\partial
x}\psi_{i}+\sum_{i\in\mathcal{I}_{\partial}}^{n}v_{i}\int\displaylimits_{\Omega}\frac{\partial\psi_{j}}{\partial
x}\psi_{i}\right)$
$\displaystyle-(r+r_{c})\left(\sum_{i=1}^{n}v_{i}\int\displaylimits_{\Omega}\psi_{j}\psi_{i}+\sum_{i\in\mathcal{I}_{\partial}}v_{i}\int\displaylimits_{\Omega}\psi_{j}\psi_{i}\right).$

In practice, systems of equations (41) and (42) are constructed via an
assembly process using (local) element matrices, whose structures depend on
the choice of the basis functions $\psi_{i}$. The common choice for the basis
functions is a class of functions satisfying the nodal condition

$\displaystyle\psi_{i}(x_{j})=\begin{cases}1,&i=j,\\\
0,&\text{otherwise},\end{cases}$ (43)

where $x_{j}$ is the nodal point. This choice leads to global systems of
linear equations with sparse and banded coefficient matrices.

###  3.1 Linear polynomial bases

Consider partition of the spatial domain $\Omega$ into $n_{E}$ non-overlapping
elements $\Omega_{j}=[x_{j-1},x_{j}]$, with $|\Omega_{j}|=x_{j}-x_{j-1}=h$,
$x_{j}$, $j=0,\dots,n_{E}$, the nodal points, $x_{0}=x_{\min}$ and
$x_{n_{E}}=x_{\max}$. In the basic element $\Omega_{j}=[x_{j-1},x_{j}]$, we
define two linear interpolation basis functions

$\displaystyle\psi_{j-1}(x)$
$\displaystyle=(x-x_{j})/(x_{j-1}-x_{j})=-(x-x_{j})/h,$
$\displaystyle\psi_{j}(x)$
$\displaystyle=(x-x_{j-1})/(x_{j}-x_{j-1})=(x-x_{j-1})/h,$

resulting in the P1 finite element (P1-FEM). Evaluating the integrals in (41)
and (42) using the above-stated basis functions over the element $\Omega_{j}$
results in the following local (element) matrices:

  * •

For the $\displaystyle\int\psi_{j}\psi_{i}dx$ term, the element matrix reads

$M_{j}=\begin{bmatrix}\displaystyle\int\displaylimits_{\Omega_{j}}\psi_{j-1}\psi_{j-1}dx&\displaystyle\int\displaylimits_{\Omega_{j}}\psi_{j-1}\psi_{j}dx\\\
\displaystyle\int\displaylimits_{\Omega_{j}}\psi_{j}\psi_{j-1}dx&\displaystyle\int\displaylimits_{\Omega_{j}}\psi_{j}\psi_{j}dx\end{bmatrix}=\frac{h}{6}\begin{bmatrix}2&1\\\
1&2\end{bmatrix}.$

  * •

For the $-\displaystyle\int\psi_{j,x}\psi_{i,x}dx$ term, the element matrix
reads

$K_{j}=-\begin{bmatrix}\displaystyle\int\displaylimits_{\Omega_{j}}\psi_{j-1,x}\psi_{j-1,x}dx&\displaystyle\int\displaylimits_{\Omega_{j}}\psi_{j-1,x}\psi_{j,x}dx\\\
\displaystyle\int\displaylimits_{\Omega_{j}}\psi_{j,x}\psi_{j-1,x}dx&\displaystyle\int\displaylimits_{\Omega_{j}}\psi_{j,x}\psi_{j,x}dx\end{bmatrix}=-\frac{1}{h}\begin{bmatrix}1&-1\\\
-1&1\end{bmatrix}.$

  * •

For the $\displaystyle\int\psi_{j}\psi_{i,x}dx$ term, the element matrix reads

$N_{j}=\begin{bmatrix}\displaystyle\int\displaylimits_{\Omega_{j}}\psi_{j-1}\psi_{j-1,x}dx&\displaystyle\int\displaylimits_{\Omega_{j}}\psi_{j-1}\psi_{j,x}dx\\\
\displaystyle\int\displaylimits_{\Omega_{j}}\psi_{j}\psi_{j-1,x}dx&\displaystyle\int\displaylimits_{\Omega_{j}}\psi_{j}\psi_{j,x}dx\end{bmatrix}=\frac{1}{2}\begin{bmatrix}-1&-1\\\
1&1\end{bmatrix}.$

###  3.2 Quadratic polynomial bases

In this approach, we add a midpoint $x_{j-\frac{1}{2}}=(x_{j-1}+x_{j})/2$ in
the basic element $\Omega_{j}$, giving three nodal points: $x_{j-1}$,
$x_{j-1/2}$, and $x_{j}$, and define three quadratic interpolation polynomials
satisfying the nodal condition (43):

$\displaystyle\psi_{j-1}(x)$
$\displaystyle=\frac{(x-x_{j-\frac{1}{2}})(x-x_{j})}{(x_{j-1}-x_{j-\frac{1}{2}})(x_{j-1}-x_{j})}=2(x-x_{j-\frac{1}{2}})(x-x_{j})/h^{2},$
(44) $\displaystyle\psi_{j-\frac{1}{2}}(x)$
$\displaystyle=\frac{(x-x_{j-1})(x-x_{j})}{(x_{j-\frac{1}{2}}-x_{j-1})(x_{j-\frac{1}{2}}-x_{j})}=-4(x-x_{j-1})(x-x_{j})/h^{2},$
(45) $\displaystyle\psi_{j}(x)$
$\displaystyle=\frac{(x-x_{j-1})(x-x_{j-\frac{1}{2}})}{(x_{j}-x_{j-1})(x_{j}-x_{j-\frac{1}{2}})}=2(x-x_{j-1})(x-x_{j-\frac{1}{2}})/h^{2},$
(46)

resulting in P2-FEM. The local element matrices are as follows:

  * •

the $\displaystyle\int\psi_{j}\psi_{i}dx$ term:

$M_{j}=\begin{bmatrix}\displaystyle\int\displaylimits_{\Omega_{j}}\psi_{j-1}\psi_{j-1}dx&\displaystyle\int\displaylimits_{\Omega_{j}}\psi_{j-1}\psi_{j-\frac{1}{2}}dx&\displaystyle\int\displaylimits_{\Omega_{j}}\psi_{j-1}\psi_{j}dx\\\
\displaystyle\int\displaylimits_{\Omega_{j}}\psi_{j-\frac{1}{2}}\psi_{j-1}dx&\displaystyle\int\displaylimits_{\Omega_{j}}\psi_{j-\frac{1}{2}}\psi_{j-\frac{1}{2}}dx&\displaystyle\int\displaylimits_{\Omega_{j}}\psi_{j-\frac{1}{2}}\psi_{j}dx\\\
\displaystyle\int\displaylimits_{\Omega_{j}}\psi_{j}\psi_{j-1}dx&\displaystyle\int\displaylimits_{\Omega_{j}}\psi_{j}\psi_{j-\frac{1}{2}}dx&\displaystyle\int\displaylimits_{\Omega_{j}}\psi_{j}\psi_{j}dx\end{bmatrix}=\frac{h}{30}\begin{bmatrix}4&2&-1\\\
2&16&2\\\ -1&2&4\end{bmatrix}.$

  * •

the $\displaystyle-\int\psi_{j,x}\psi_{i,x}dx$ term:

$K_{j}=-\begin{bmatrix}\displaystyle\int\displaylimits_{\Omega_{j}}\psi_{j-1,x}\psi_{j-1,x}dx&\displaystyle\int\displaylimits_{\Omega_{j}}\psi_{j-1,x}\psi_{j-\frac{1}{2},x}dx&\displaystyle\int\displaylimits_{\Omega_{j}}\psi_{j-1,x}\psi_{j,x}dx\\\
\displaystyle\int\displaylimits_{\Omega_{j}}\psi_{j-\frac{1}{2},x}\psi_{j-1,x}dx&\displaystyle\int\displaylimits_{\Omega_{j}}\psi_{j-\frac{1}{2},x}\psi_{j-\frac{1}{2},x}dx&\displaystyle\int\displaylimits_{\Omega_{j}}\psi_{j-\frac{1}{2},x}\psi_{j,x}dx\\\
\displaystyle\int\displaylimits_{\Omega_{j}}\psi_{j,x}\psi_{j-1,x}dx&\displaystyle\int\displaylimits_{\Omega_{j}}\psi_{j,x}\psi_{j-\frac{1}{2},x}dx&\displaystyle\int\displaylimits_{\Omega_{j}}\psi_{j,x}\psi_{j,x}dx\end{bmatrix}=-\frac{1}{3h}\begin{bmatrix}7&-8&1\\\
-8&16&-8\\\ 1&-8&7\end{bmatrix}.$

  * •

the $\displaystyle\int\psi_{j}\psi_{i,x}dx$ term:

$N_{j}=\begin{bmatrix}\displaystyle\int\displaylimits_{\Omega_{j}}\psi_{j-1}\psi_{j-1,x}dx&\displaystyle\int\displaylimits_{\Omega_{j}}\psi_{j-1}\psi_{j-\frac{1}{2},x}dx&\displaystyle\int\displaylimits_{\Omega_{j}}\psi_{j-1}\psi_{j,x}dx\\\
\displaystyle\int\displaylimits_{\Omega_{j}}\psi_{j-\frac{1}{2}}\psi_{j-1,x}dx&\displaystyle\int\displaylimits_{\Omega_{j}}\psi_{j-\frac{1}{2}}\psi_{j-\frac{1}{2},x}dx&\displaystyle\int\displaylimits_{\Omega_{j}}\psi_{j-\frac{1}{2}}\psi_{j,x}dx\\\
\displaystyle\int\displaylimits_{\Omega_{j}}\psi_{j}\psi_{j-1,x}dx&\displaystyle\int\displaylimits_{\Omega_{j}}\psi_{j}\psi_{j-\frac{1}{2},x}dx&\displaystyle\int\displaylimits_{\Omega_{j}}\psi_{j}\psi_{j,x}dx\end{bmatrix}=\frac{1}{6}\begin{bmatrix}-3&-4&1\\\
4&0&-4\\\ -1&4&3\end{bmatrix}.$

###  3.3 Treating the constraints by Penalty method

We now turn to the two nonlinear penalty terms in (38) and construct a finite-
element approximation to them. We in particular apply group finite element
[15] to deal with the nonlinearity. We shall discuss the construction for
$\mathcal{P}_{call,j}$; construction of finite element approximation for
$\mathcal{P}_{put,j}$ is done in the same way.

We assume that the term $\zeta_{call}:=\alpha_{call}(U-U^{\star}_{call})$ is
approximated by

$\zeta_{call}=\sum_{i=1}^{n}\zeta_{i}\psi_{i}+\sum_{\mathcal{I}_{\partial}}\zeta_{i}\psi_{i},$

where
$\zeta_{call,i}=\alpha_{call}(x_{i})(U(x_{i})-U^{\star}_{call}(x_{i}))=:\alpha_{call,i}(u_{i}-u^{\star}_{call,i})$.
Therefore, for $w=\psi_{j}$, $j=1,\dots,n$, we have

$\displaystyle\displaystyle\mathcal{P}_{call,j}$
$\displaystyle=\rho\int\displaylimits_{\Omega}\psi_{j}\left(\sum_{i=1}^{n}\zeta_{i}\psi_{i}+\sum_{\mathcal{I}_{\partial}}\zeta_{i}\psi_{i}\right)dx=\rho\left(\sum_{i=1}^{n}\zeta_{i}\int\displaylimits_{\Omega}\psi_{j}\psi_{i}dx+\sum_{\mathcal{I}_{\partial}}\zeta_{i}\int\displaylimits_{\Omega}\psi_{j}\psi_{i}dx\right)$
$\displaystyle=\rho\left(\sum_{i=1}^{n}\alpha_{call,i}(u_{i}-u^{\star}_{call,i})\int\displaylimits_{\Omega}\psi_{j}\psi_{i}dx+\sum_{\mathcal{I}_{\partial}}\alpha_{call,i}(u_{i}-u^{\star}_{call,i})\int\displaylimits_{\Omega}\psi_{j}\psi_{i}dx\right).$
(47)

Each integral in the above-equation is evaluated element-wise, resulting in
the local element matrices $M_{j}$, $K_{j}$, and $N_{j}$, given in Sections
3.1 and 3.2.

By using the same argument,

$\displaystyle\displaystyle\mathcal{P}_{put,j}=\rho\left(\sum_{i=1}^{n}\alpha_{put,i}(u_{i}-u^{\star}_{put,i})\int\displaylimits_{\Omega}\psi_{j}\psi_{i}dx+\sum_{\mathcal{I}_{\partial}}\alpha_{put,i}(u_{i}-u^{\star}_{put,i})\int\displaylimits_{\Omega}\psi_{j}\psi_{i}dx\right).$
(48)

##  4 Time integration scheme

The global finite-element system obtained from assembling the local finite-
element matrices can be represented by the differential algebraic equations
(DAEs):

$\displaystyle\frac{\partial}{\partial\tau}(M\boldsymbol{u}+\hat{\boldsymbol{b}}_{M,u})$
$\displaystyle=-\frac{\sigma^{2}}{2}K\boldsymbol{u}-\left(r-\frac{\sigma^{2}}{2}\right)N\boldsymbol{u}-rM\boldsymbol{u}-r_{c}M\boldsymbol{v}-\boldsymbol{\beta}_{1}(\boldsymbol{u},\boldsymbol{v})$
$\displaystyle+\rho MP_{put}(\boldsymbol{u}_{put}^{\star}-\boldsymbol{u})+\rho
MP_{call}(\boldsymbol{u}-\boldsymbol{u}_{call}^{\star})+\rho\boldsymbol{b}_{put}+\rho\boldsymbol{b}_{call}:=F_{1}(\boldsymbol{u},\boldsymbol{v}),$
(49)
$\displaystyle\frac{\partial}{\partial\tau}(M\boldsymbol{v}+\hat{\boldsymbol{b}}_{M,v})$
$\displaystyle=-\frac{\sigma^{2}}{2}K\boldsymbol{v}-\left(r-\frac{\sigma^{2}}{2}\right)N\boldsymbol{v}-(r+r_{c})M\boldsymbol{v}-\boldsymbol{\beta}_{2}(\boldsymbol{v}):=F_{2}(\boldsymbol{v}),$
(50)

where $P_{put}=\text{diag}(\alpha_{put,j})$,
$P_{call}=\text{diag}(\alpha_{call,j})$, and

$\displaystyle\displaystyle\boldsymbol{\beta_{1}}(\boldsymbol{u},\boldsymbol{v})$
$\displaystyle=\frac{\sigma^{2}}{2}\boldsymbol{b}_{K,u}+\left(r-\frac{\sigma^{2}}{2}\right)\boldsymbol{b}_{N,u}+r\boldsymbol{b}_{M,u}+r_{c}\boldsymbol{b}_{M,v},$
(51) $\displaystyle\displaystyle\boldsymbol{\beta}_{2}(\boldsymbol{v})$
$\displaystyle=\frac{\sigma^{2}}{2}\boldsymbol{b}_{K,v}+\left(r-\frac{\sigma^{2}}{2}\right)\boldsymbol{b}_{N,v}+(r+r_{c})\boldsymbol{b}_{M,v},$
(52)

are the boundary condition vectors.

Time integration of the DAEs (49) and (50) is carried out by applying the
$\theta$-scheme on both equations, which results in the systems, with
$\theta\in[0,1]$, $\Delta\tau=T/n_{\tau}$, and $n_{\tau}$ the number of time
steps,

$\displaystyle
M\boldsymbol{u}^{m+1}+\hat{\boldsymbol{b}}_{M,u}^{m+1}-M\boldsymbol{u}^{m}-\hat{\boldsymbol{b}}_{M,u}^{m}=\theta\Delta\tau
F_{1}(\boldsymbol{u}^{m+1},\boldsymbol{v}^{m+1})+(1-\theta)\Delta\tau
F_{1}(\boldsymbol{u}^{m},\boldsymbol{v}^{m}),$ $\displaystyle
M\boldsymbol{v}^{m+1}+\hat{\boldsymbol{b}}_{M,v}^{m+1}-M\boldsymbol{v}^{m}-\hat{\boldsymbol{b}}_{M,v}^{m}=\theta\Delta\tau
F_{2}(\boldsymbol{v}^{m+1})+(1-\theta)\Delta\tau F_{2}(\boldsymbol{v}^{m}),$

or

$\displaystyle A_{11}\boldsymbol{u}^{m+1}$
$\displaystyle+A_{12}\boldsymbol{v}^{m+1}-\rho\theta\Delta\tau
M\left(P_{put}^{m+1}(\boldsymbol{u}_{put}^{\star,m+1}-\boldsymbol{u}^{m+1})+P_{call}^{m+1}(\boldsymbol{u}^{m+1}-\boldsymbol{u}^{\star,m+1}_{call})\right)$
$\displaystyle=\widetilde{A}_{11}\boldsymbol{u}^{m}+\widetilde{A}_{12}\boldsymbol{v}^{m}+\rho(1-\theta)\Delta\tau
M\left(P_{put}^{m}(\boldsymbol{u}_{put}^{\star,m}-\boldsymbol{u}^{m})+P_{call}^{m}(\boldsymbol{u}^{m}-\boldsymbol{u}^{\star,m}_{call})\right)$
$\displaystyle+\theta\Delta\tau\boldsymbol{\beta}_{1}^{m+1}+(1-\theta)\Delta\tau\boldsymbol{\beta}_{1}^{m}+\hat{\boldsymbol{b}}^{m}_{M,u}-\hat{\boldsymbol{b}}^{m+1}_{M,u}+\theta\rho\Delta\tau(\boldsymbol{b}_{put}^{m+1}+\boldsymbol{b}_{call}^{m+1})+(1-\theta)\rho\Delta\tau(\boldsymbol{b}_{put}^{m}+\boldsymbol{b}_{call}^{m}),$
(53) $\displaystyle A_{22}\boldsymbol{v}^{m+1}$
$\displaystyle=\widetilde{A}_{22}\boldsymbol{v}^{m}+\theta\Delta\tau\boldsymbol{\beta}_{2}^{m+1}+(1-\theta)\Delta\tau\boldsymbol{\beta}_{2}^{m}+\hat{\boldsymbol{b}}^{m}_{M,v}-\hat{\boldsymbol{b}}^{m+1}_{M,v},$
(54)

where

$\displaystyle A_{11}$
$\displaystyle=M+\theta\Delta\tau\left(\frac{\sigma^{2}}{2}K+\left(r-\frac{\sigma^{2}}{2}\right)N+rM\right),$
$\displaystyle A_{12}$ $\displaystyle=\theta\Delta\tau r_{c}M,$ $\displaystyle
A_{22}$
$\displaystyle=M+\theta\Delta\tau\left(\frac{\sigma^{2}}{2}K+\left(r-\frac{\sigma^{2}}{2}\right)N+(r+r_{c})M\right)$
$\displaystyle\widetilde{A}_{11}$
$\displaystyle=M-(1-\theta)\Delta\tau\left(\frac{\sigma^{2}}{2}K+\left(r-\frac{\sigma^{2}}{2}\right)N+rM\right),$
$\displaystyle\widetilde{A}_{12}$ $\displaystyle=-(1-\theta)\Delta\tau
r_{c}M,$ $\displaystyle\widetilde{A}_{22}$
$\displaystyle=M-(1-\theta)\Delta\tau\left(\frac{\sigma^{2}}{2}K+\left(r-\frac{\sigma^{2}}{2}\right)N+(r+r_{c})M\right).$

Let the solutions $\boldsymbol{u}^{m}$ and $\boldsymbol{v}^{m}$ be known. The
solutions at the next time level $m+1$ can in principle be computed by first
solving (54) for $\boldsymbol{v}^{m+1}$. The solution $\boldsymbol{u}^{m+1}$
is then computed via (53) using the known $\boldsymbol{u}^{m}$,
$\boldsymbol{v}^{m}$, and $\boldsymbol{v}^{m+1}$. This procedure however
requires knowledge of the solutions at the boundaries at the time level $m+1$.

###  4.1 Boundary solutions

At $x_{\min}$, with $u_{0}(\tau):=U(x_{\min},\tau)$,
$v_{0}(\tau):=V(x_{\min},\tau)$, etc, the boundary conditions with penalty in
$U$ can be written as follows:

$\displaystyle\begin{cases}\displaystyle\frac{\partial
u_{0}(\tau)}{\partial\tau}=-ru_{0}(\tau)-r_{c}v_{0}(\tau)+\rho\max\left(u_{0}(\tau)-u^{\star}_{call,0}(\tau),0\right)+\rho\max\left(u^{\star}_{put,0}(\tau)-u_{0}(\tau),0\right),\\\\[10.0pt]
\displaystyle\frac{\partial v_{0}(\tau)}{\partial
t}=-(r+r_{c})v_{0}(\tau).\end{cases}$ (55)

Note that
$\max\left(u_{0}(\tau)-u^{\star}_{call,0}(\tau),0\right)=p_{call,0}(\tau)(u_{0}(\tau)-u^{\star}_{call,0}(\tau))$,
with

$p_{call,0}=\begin{cases}1,&\text{if }u_{0}(\tau)>u^{\star}_{call,0}(\tau),\\\
0,&\text{otherwise},\end{cases}$

and
$\max\left(u^{\star}_{put,0}(\tau)-u_{0}(\tau),0\right)=p_{put,0}(\tau)(u^{\star}_{put,0}(\tau)-u_{0}(\tau))$,
with

$p_{put,0}=\begin{cases}1,&\text{if }u_{0}(\tau)<u^{\star}_{put,0}(\tau),\\\
0,&\text{otherwise}.\end{cases}$

Application of the $\theta$-scheme on (55) leads to the discrete equations:

$\displaystyle u_{0}^{m+1}+\theta\Delta\tau$
$\displaystyle\left(ru_{0}^{m+1}+r_{c}v_{0}^{m+1}+\rho\left(p_{call,0}^{m+1}(u_{0}^{m+1}-u_{call,0}^{\star,m+1})+p_{put}^{m+1}(u_{put,0}^{\star,m+1}-u_{0}^{m+1})\right)\right)$
$\displaystyle=u_{0}^{m}-(1-\theta)\Delta\tau\left(ru_{0}^{m}+r_{c}v_{0}^{m}+\rho\left(p_{call,0}^{m}(u_{0}^{m}-u_{call,0}^{\star,m})+p_{put}^{m}(u_{put,0}^{\star,m}-u_{0}^{m})\right)\right),$
(56) $\displaystyle\left[1+\theta\Delta\tau(r+r_{c})\right]v_{0}^{m+1}$
$\displaystyle=[1-(1-\theta)\Delta\tau(r+r_{c})]v_{0}^{m}.$ (57)

Let the boundary solution $v_{0}^{m}$ be known. Then $v_{0}^{m+1}$ can be
computed from (57). With $u_{0}^{m}$, $v_{0}^{m}$, and $v_{0}^{m+1}$ now
known, (56) becomes a nonlinear function of $u_{0}^{m+1}$, which can be solved
approximately using Newton’s method. First of all we assume that the penalty
term in (56) is approximated in a fully implicit way at the new time level
$m+1$. This results in the equation

$\displaystyle 0=(1+\theta\Delta\tau
r)u_{0}^{m+1}+\Delta\tau\rho\left(p_{call,0}^{m+1}(u_{0}^{m+1}-u_{call,0}^{\star,m+1})+p_{put}^{m+1}(u_{put,0}^{\star,m+1}-u_{0}^{m+1})\right)-\phi_{0}=:f(u_{0}^{m+1}),$
(58)

where

$\phi_{0}=u_{0}^{m}-(1-\theta)\Delta\tau\left(ru_{0}^{m}+r_{c}v_{0}^{m}\right)-\theta\Delta\tau
r_{c}v_{0}^{m+1}.$

With

$\displaystyle f^{\prime}(u_{0}^{m+1})=1+\theta\Delta\tau
r+\Delta\tau\rho\left(p_{call,0}^{m+1}-p_{put,0}^{m+1}\right),$ (59)

Newton’s method for finding $u_{0}^{m+1}$ satisfying (58) can be written as
follows: with an initial guess $u_{0}^{m+1,0}$, compute
$u_{0}^{m+1,k}=u_{0}^{m+1,k-1}-f(u_{0}^{m+1,k-1})/f^{\prime}(u_{0}^{m+1,k-1})$,
for $k=1,2,\dots$.

In this paper, $u_{0}^{m+1,0}$ is chosen to be solution of the unconstrained
boundary problem given by (17). Since we expect that $u_{0}^{m+1,0}$ computed
in this way is a better approximation than, e.g., $u_{0}^{m}$, we can use this
value as well to evaluate the conditions to constraint $v$.

The complete algorithm for computing $u_{0}^{m+1}$ and $v_{0}^{m+1}$ is as
follows:

Algorithm 1: Computing the boundary solutions

  1. 1.

input $u_{0}^{m}$, $v_{0}^{m}$;

  2. 2.

compute $B_{p}(\tau^{m+1})$ and $B_{c}(\tau^{m+1})$;

  3. 3.

compute $v_{0}^{m+1}$ from (57);

  4. 4.

compute $u_{0}^{m+1}$ from (56) without penalty terms;

  5. 5.

apply constraints on $v_{0}^{m+1}$ using $u_{0}^{m+1}$;

  6. 6.

set $u_{0}^{m+1,0}\leftarrow u_{0}^{m+1}$;

  7. 7.

for $k=1,2,\dots$ until convergence

    1. 7.1.

compute $f(u_{0}^{m+1,k-1})$ using (58);

    2. 7.2.

compute $f^{\prime}(u_{0}^{m+1,k-1})$ using (59);

    3. 7.3.

$u_{0}^{m+1,k}\leftarrow
u_{0}^{m+1,k-1}-f(u_{0}^{m+1,k-1})/f^{\prime}(u_{0}^{m+1,k-1})$;

  8. 8.

apply constraints on $v_{0}^{m+1}$ using $u_{0}^{m+1,k}$;

  9. 9.

if $\tau^{m+1}\in\mathcal{T}_{\text{coupon}}$

    1. 9.1

$u_{0}^{m+1}\leftarrow u_{0}^{m+1}+K$;

    2. 9.2

$v_{0}^{m+1}\leftarrow v_{0}^{m+1}+K$;

At $x_{\max}$ we need to compute $u_{n+1}^{m+1}:=U(x_{\max},\tau^{m+1})$ and
$v_{n+1}^{m+1}:=V(x_{\max},\tau^{m+1})$ via (19), apply the constraints, and
pay the coupon if $\tau^{m+1}\in\mathcal{T}_{\text{coupon}}$.

We remark here that solution of Newton’s method exists if the derivative
$f^{\prime}(u_{0}^{m})\neq 0$. This condition however cannot be theoretically
guaranteed. In this regard, for the tuple
$(p_{call,0}^{m},p_{put,0}^{m})\in\\{0,1\\}^{2}$ and nonnegative $\theta$,
$\Delta\tau$, $\rho$, and $r$,

  1. 1.

if $p_{call,0}=p_{put,0}$, then $f^{\prime}(u_{0}^{m})=1+\theta\Delta\tau
r>0$;

  2. 2.

if $p_{call,0}\neq p_{put,0}$, we have two situations:

    1. (a)

if $p_{call,0}=1$ and $p_{put,0}=0$, then
$f^{\prime}(u_{0}^{m})=1+\Delta\tau(\theta r+\rho)>0$;

    2. (b)

if $p_{call,0}=0$ and $p_{put,0}=1$, then
$f^{\prime}(u_{0}^{m})=1+\Delta\tau(\theta r-\rho)=0$ if
$\Delta\tau=1/(\rho-\theta r)$.

Since the penalty parameter $\rho$ is taken very large (e.g., $\rho=10^{12}$
in this paper), vanishing $f^{\prime}$ is very unlikely to happen under a
reasonable choice of the time step $\Delta\tau$. For instance, using the
simulation parameters in Table 1 in Section 5, the vanishing $f^{\prime}$
situation outlined in Point 2(b) above may occur if $\Delta\tau$ is chosen to
be in the order of $\rho^{-1}=10^{-12}$, which makes the time-integration
process extremely impractical.

###  4.2 Interior solutions

With solutions at the boundaries available at $\tau^{m}$ and $\tau^{m+1}$, all
related boundary vectors in (53) and (54) are known. We are thus now in the
position to compute the solutions $\boldsymbol{u}^{m+1}$ and
$\boldsymbol{v}^{m+1}$. $\boldsymbol{v}^{m+1}$ is readily computed from (54).
With $\boldsymbol{v}^{m+1}$ known, (53) reduces to a nonlinear function of
$\boldsymbol{u}^{m+1}$. As we do for the boundary equations, we assume that
the penalty terms in (54) are approximated in a fully implicit way, resulting
in the equation

$\displaystyle 0$ $\displaystyle=A_{11}\boldsymbol{u}^{m+1}-\rho\Delta\tau
M\left(P_{put}^{m+1}(\boldsymbol{u}_{put}^{\star,m+1}-\boldsymbol{u}^{m+1})+P_{call}^{m+1}(\boldsymbol{u}^{m+1}-\boldsymbol{u}^{\star,m+1}_{call})\right)-\boldsymbol{\phi}:=\boldsymbol{f}(\boldsymbol{u}^{m+1}),$
(60)

where

$\displaystyle\boldsymbol{\phi}=\widetilde{A}_{11}\boldsymbol{u}^{m}+\widetilde{A}_{12}\boldsymbol{v}^{m}-A_{12}\boldsymbol{v}^{m+1}+\theta\Delta\tau\boldsymbol{\beta}_{1}^{m+1}+(1-\theta)\Delta\tau\boldsymbol{\beta}_{1}^{m}+\hat{\boldsymbol{b}}^{m}_{M,u}-\hat{\boldsymbol{b}}^{m+1}_{M,u}+\theta\Delta\tau(\boldsymbol{b}_{put}^{m+1}+\boldsymbol{b}_{call}^{m+1}).$
(61)

The nonlinear equation (60) is solved iteratively using Newton’s method.
Starting from an initial guess of the solution $\boldsymbol{u}^{m+1,0}$, the
solution $\boldsymbol{u}^{m+1}$ is approximated using the iterands

$\boldsymbol{u}^{m+1,k}\leftarrow\boldsymbol{u}^{m+1,k-1}-\left(\nabla\boldsymbol{f}(\boldsymbol{u}^{m+1,k-1})\right)^{-1}\boldsymbol{f}(\boldsymbol{u}^{m+1,k-1}),\quad
k=1,2,\dots$

where $\nabla\boldsymbol{f}(\boldsymbol{u}^{m+1,k-1})$ is the Jacobian of
$\boldsymbol{f}$, given by

$\displaystyle\nabla\boldsymbol{f}(\boldsymbol{u}^{m+1,k-1})=A_{11}+\rho\Delta\tau
M\left(P_{call}^{m+1,k-1}-P_{put}^{m+1,k-1}\right).$ (62)

The initial guess $\boldsymbol{u}^{m+1,0}$ is chosen such that it solves
unconstrained CB PDE, which is equivalent to solving (53) without penalty
terms. We shall also use this unconstrained solution $\boldsymbol{u}^{m+1,0}$
to constrain the initially computed $\boldsymbol{v}^{m+1}$ prior to the start
of Newton’s iterations. The procedure for computing the solutions in the
interior after one $\theta$-scheme step is summarized in the following
algorithm:

Algorithm 2: Computing the interior solutions

  1. 1.

input $\boldsymbol{u}^{m}$, $\boldsymbol{v}^{m}$;

  2. 2.

compute $B_{p}(\tau^{m+1})$ and $B_{c}(\tau^{m+1})$;

  3. 3.

compute $\boldsymbol{v}^{m+1}$ from (54);

  4. 4.

compute $\boldsymbol{u}^{m+1}$ from (53) without penalty terms;

  5. 5.

apply constraints on $\boldsymbol{v}^{m+1}$ using $\boldsymbol{u}^{m+1}$;

  6. 6.

set $\boldsymbol{u}^{m+1,0}\leftarrow\boldsymbol{u}^{m+1}$;

  7. 7.

for $k=1,2,\dots$ until convergence

    1. 7.1.

compute $\boldsymbol{f}(\boldsymbol{u}^{m+1,k-1})$ using (60);

    2. 7.2.

compute $\nabla\boldsymbol{f}(\boldsymbol{u}^{m+1,k-1})$ using (62);

    3. 7.3.

$\boldsymbol{u}^{m+1,k}\leftarrow\boldsymbol{u}^{m+1,k-1}-\left(\nabla\boldsymbol{f}(\boldsymbol{u}^{m+1,k-1})\right)^{-1}\boldsymbol{f}(\boldsymbol{u}^{m+1,k-1})$;

  8. 8.

apply constraints on $\boldsymbol{v}^{m+1}$ using $\boldsymbol{u}^{m+1,k}$;

  9. 9.

if $\tau^{m+1}\in\mathcal{T}_{\text{coupon}}$

    1. 9.1

$\boldsymbol{u}^{m+1}\leftarrow\boldsymbol{u}^{m+1}+K$;

    2. 9.2

$\boldsymbol{v}^{m+1}\leftarrow\boldsymbol{v}^{m+1}+K$;

The existence of a solution of Newton’s method requires nonsingularity of the
Jacobian $\nabla f$. Like in the case with Newton’s method for computing
boundary solutions, there may exist values of parameters $r$, $h$, $\sigma$,
$\theta$, $\Delta\tau$, and $\rho$ such that the Jacobian is singular.
Quantification of such conditions requires analysis, which is beyond the scope
of this paper. Numerical tests using various realistic values of parameters
exhibit no convergence issues with the methods, suggesting nonsingularity of
$\nabla f$.

###  4.3 Summary of the time integration method

By including Algorithm 1 and 2 in the $\theta$-scheme, the complete time
stepping procedure to compute the solutions $\boldsymbol{u}(x,\tau)$ and
$\boldsymbol{v}(x,\tau)$ is shown in Algorithm 3.

Algorithm 3: $\theta$-scheme for time integration with constraints

  1. 1.

input the initial solution at $\tau^{0}=0$: $\boldsymbol{u}^{0}$,
$\boldsymbol{v}^{0}$, $u_{0}^{0}$, $u_{n+1}^{0}$, $v_{0}^{0}$, and
$v_{n+1}^{0}$;

  2. 2.

for $m=0,1,\dots,n_{\tau}-1$

    1. 2.1.

Compute $B_{p}(\tau^{m+1})$ and $B_{c}(\tau^{m+1})$;

    2. 2.2.

Compute $u_{0}^{m+1}$ and $v_{0}^{m+1}$ by performing Step 3–9 of Algorithm 1;

    3. 2.3.

Compute $u_{n+1}^{m+1}=\kappa S_{\text{int}}e^{x_{max}}$ and
$v_{n+1}^{m+1}=0$;

    4. 2.4.

Apply the conversion-callability-puttability constraints on $u_{n+1}^{m+1}$
and $v_{n+1}^{m+1}$.

    5. 2.5.

if $\tau^{m+1}\in\mathcal{T}_{coupon}$

      1. 2.5.1.

$u_{n+1}^{m+1}\leftarrow u_{n+1}^{m+1}+K$;

      2. 2.5.2.

$v_{n+1}^{m+1}\leftarrow v_{n+1}^{m+1}+K$;

    6. 2.6.

Compute $\boldsymbol{u}^{m+1}$ and $\boldsymbol{v}^{m+1}$ by performing Step
3–9 of Algorithm 2;

For increased stability, it is possible to initiate the time integration using
Rannacher’s step [24]. In this case, for $m=0$, the solutions
$\boldsymbol{u}^{1}$, $\boldsymbol{v}^{1}$, $u_{0}^{1}$, $u_{n+1}^{1}$,
$v_{0}^{1}$, and $v_{n+1}^{1}$ are computed using the initial conditions using
Step 2.1–2.6 but with a smaller time step than $\Delta\tau$ (e.g,
$\Delta\tau_{rann}=\Delta\tau/n_{rann}$, where $1<n_{rann}\in\mathbb{N}$).
Since we did not see stability issues when $\theta=1/2$, we did not implement
Rannacher’s step to obtain numerical results presented in Section 5.

##  5 Numerical solution of the TF model

###  5.1 Comparison with FDM

In this section, we present numerical results obtained from the FEM and time
integration method discussed in Sections 3 and 4. The modeling and
computational parameters are summarized in Table 1, taken from [6, 16]. The
numerical solution is computed for $x\in[-6,2]$, so that the left boundary is
sufficiently close to $0$ in the $S$-space. We compare the numerical results
with a second-order finite difference method (FDM), combined with the
$\theta$-scheme for time integration. To have a fair comparison, we set the
number of unknowns in both methods to be equal. For instance, the same meshing
can be used for P1-FEM and FDM. For P2-FEM, due to additional unknowns at the
midpoint of each element, we double the number of gridpoints in the FDM
meshing. Throughout, $n_{E}$ and $n_{t}=n_{\tau}$ denote the number of
elements and time steps, respectively. Thus, for P1-FEM, the number of
unknowns (nodal points) is $n=n_{E}-1$ and, for P2-FEM, $n=2n_{E}-1$. For the
$\theta$-scheme, we set $\theta=\frac{1}{2}$ (Crank-Nicolson).

Table 1: Modeling and computational parameters Parameter |  Value
---|---
Time to maturity $T$ |  5 years
Conversion |  0 to 5 years into $k$ shares
Conversion ratio $k$ |  1.0
Face Value $F$ |  $100
Clean call price $B_{c}$ |  $110, from Year 3 to Year 5 ($t\in(3,5]$)
Clean put price $B_{p}$ |  $105, during Year 3 ($t\in(2,3]$)
Coupon payments $K$ |  $4.0
Coupon dates |  .5, 1.0, 1.5, … ,5.0 (semiannual)
Risk-free interest rate r |  5% or 0.05
Credit risk $r_{c}$ |  2% or 0.02
Volatility $\sigma$ |  20% or 0.2
Underlying stock price at $t=0$, $S_{int}$ |  $100
Penalty parameter $\rho$ |  $10^{12}$
Newton-Raphson’s method tolerance $tol$ |  $10^{-12}$

Figures 1, 2, and 3 show surfaces of the CB price $U$ for $40\leq S\leq 160$,
computed using $n_{E}=$ 100 and 200 and with $n_{t}=200$. We note here that
the steps in the solution surfaces correspond to the coupon payment time.
Qualitatively, the surfaces indicate no noticeable difference between
solutions of FDM and P1- or P2-FEM. That the FDM and FEM are qualitatively not
distinguishable can be seen from Figure 4 and 5, where the solutions at $t=0$
are plotted.

Figure 1: Finite difference solution of the penalty TF model over time
$t\in[0,5]$, with $r=0.05$, $r_{c}=0.02$, $\sigma=0.2$, $F=\$100$, $K=\$4$,
$\rho=10^{12}$; Left figure: $n=100$, $n_{t}=100$; Right figure: $n=200$,
$n_{t}=200$.

Figure 2: P1-FEM solution of the penalty TF model over time $t\in[0,5]$, with
$r=0.05$, $r_{c}=0.02$, $\sigma=0.2$, $F=\$100$, $K=\$4$, $\rho=10^{12}$; Left
figure: $n_{E}=100$, $n_{t}=100$; Right figure: $n_{E}=200$, $n_{t}=200$.

Figure 3: P2-FEM solution of the penalty TF model at $t\in[0,5]$, with
$r=0.05$, $r_{c}=0.02$, $\sigma=0.2$, $F=\$100$, $K=\$4$, $\rho=10^{12}$; Left
figure: $n_{E}=100$, $n_{t}=100$; Right figure: $n_{E}=200$, $n_{t}=200$.

Figure 4: Solution of the TF model at $t=0$, with $r=0.05$,$r_{c}=0.02$,
$\sigma=0.2$, $F=\$100$, $K=\$4$, $\rho=10^{12}$; Left figure: for P1-FEM and
P2-FEM, $n_{E}=100$, $n_{t}=100$, and for FDM $n=200$; Right figure: for
P1-FEM and P2-FEM, $n_{E}=200$, $n_{t}=200$, for FDM, $n=400$.

Figure 5: Solution of the TF model at $t=0$, with $r=0.05$,$r_{c}=0.02$,
$\sigma=0.2$, $F=\$100$, $K=\$4$, $\rho=10^{12}$; Left figure: for P1-FEM and
P2-FEM, $n_{E}=400$, $n_{t}=400$, for FDM $n=800$; Right figure: for P1-FEM
and P2-FEM, $n_{E}=800$, $n_{t}=800$, for FDM $n=1600$.

For convergence(mesh-independence) and comparison of the FE solutions with the
FD solutions obtained using various spatial meshes, we present the bond price
$U$ at the initial time $t=0$ and $S=100$ in Table 2 computed using FDM,
P1-FEM, and P2-FEM, with $n_{E}=n_{t}$ to keep the ratio $\Delta\tau/h$
constant ($\Delta\tau/h=0.625<1$). The standard stability condition for the
convection-diffusion equation discretized by FDM with explicit scheme
($\theta=0$) requires that

$\Delta\tau\leq
C_{\theta=0}\min\left(\frac{h}{|a|},\frac{h^{2}}{2\epsilon}\right),\quad
0<C_{\theta=0}\leq 1,$

where, for the TF model, $a=r-\sigma^{2}/2$ and $\epsilon=\sigma^{2}/2$. With
$n_{E}=n_{t}$, there exists however some constant $n_{E}^{*}$ such that the
above stability condition is violated for $n_{E}>n_{E}^{*}$. Instability
however was not observed in our numerical simulation with $\theta=1/2$.
Results in Table 2 suggest that both P2-FEM and FDM solutions converge to a
value around 124.78 (up to 2 decimal places) as the time-spatial mesh is
refined. The P1-FEM solution exhibits a similar convergence towards 124.78 but
at lower rate. As, can be observed that our FE solutions compared favorably
with the FD solutions across these spatial meshes.

Table 2: Convertible bond price at $t=0$ and $S=100$ computed by P1-FEM, P2-FEM, and FDM, with $n_{E}=n_{t}$ $n_{E}$ | P1-FEM | P2-FEM | FDM
---|---|---|---
100 | 124.422 | 124.846 | 124.991
200 | 124.653 | 124.820 | 124.848
400 | 124.740 | 124.814 | 124.814
600 | 124.754 | 124.781 | 124.792
800 | 124.742 | 124.781 | 124.790
1000 | 124.746 | 124.775 | 124.780
1200 | 124.758 | 124.777 | 124.781

###  5.2 The Greeks

In the literature, because FDM is mostly used to compute the solutions, the
calculation of the Greeks is done approximately using finite differencing [12,
16] at grid points. Finite difference approaches can also in principle be used
to approximate the Greeks, once the solutions are available via FEM. In this
paper, we shall compute some of the Greeks, i.e., the Delta ($\Delta$) and the
Gamma ($\Gamma$), directly from the finite-element approximation functions. We
note, however, that since the approximating P1 FE function is piece-wise
continuous and non-differentiable at the element boundaries, we can in this
case only compute the Greeks at any point inside the elements.

For P1-FEM, the approximating function is piece-wise linear continuous.
Therefore, $\Delta$ can be explicitly computed by differentiating the finite
element solution. In the element $\Omega_{j}$, the FEM solution is given by

$\displaystyle U_{h}(x)=-u_{j-1}(x-x_{j})/h+u_{j}(x-x_{j-1})/h,\quad
x\in\Omega_{j}.$ (63)

Therefore, at an instant time $t$,

$\Delta(S_{j-1/2})=\frac{\partial U_{h}}{\partial
S}\bigg{|}_{S_{j-1/2}}=\frac{S_{\text{int}}}{S_{j-1/2}}\frac{\partial
U_{h}}{\partial
x}\bigg{|}_{x_{j-1/2}}=\frac{S_{\text{int}}}{hS_{j-1/2}}(u_{j}-u_{j-1})$

after using the change of variables, $S_{j-1/2}=S_{\text{int}}e^{x_{j-1/2}}$,
$x_{j-1/2}=(x_{j-1}+x_{j})/2$.

The Greeks corresponding to higher-order derivatives however vanish. In this
case, we have to resort to finite differencing, using the FEM solutions at
nodal points. For instance,

$\frac{\partial^{2}U_{h}}{\partial
x^{2}}\bigg{|}_{x_{j-1}}\simeq(u_{j}-2u_{j-1}+u_{j-2})/h^{2}.$

Hence,

$\displaystyle\Gamma(S_{j-1})=\frac{\partial^{2}U_{h}}{\partial
S^{2}}\bigg{|}_{S_{j-1}}=\frac{S_{\text{int}}^{2}}{S^{2}}\frac{\partial^{2}U_{h}}{\partial
x^{2}}\bigg{|}_{x_{j-1}}=\frac{S_{\text{int}}^{2}}{h^{2}S^{2}}(u_{j}-2u_{j-1}+u_{j-2}).$
(64)

For P2-FEM, the approximating FE function is differentiable, which in the
element $\Omega_{j}$ reads

$\displaystyle
U_{h}(x)=2u_{j-1}(x-x_{j-\frac{1}{2}})(x-x_{j})/h^{2}-4u_{j-1/2}(x-x_{j-1})(x-x_{j})/h^{2}+2u_{j}(x-x_{j-1})(x-x_{j-\frac{1}{2}})/h^{2},$

which, after differentiation at $x_{j-1/2}$, yields

$\displaystyle\frac{\partial U_{h}}{\partial
x}\bigg{|}_{x_{j-1/2}}=(u_{j}-u_{j-1})/h\text{ and
}\frac{\partial^{2}U_{h}}{\partial x^{2}}\bigg{|}_{x_{j-1/2}}$
$\displaystyle=4(u_{j}-2u_{j-1/2}+u_{j-1})/h^{2}.$

This in turn results in $\Delta(S_{j-1/2})$ in (63) and

$\displaystyle\Gamma(S_{j-1/2})=\frac{4S_{\text{int}}^{2}}{h^{2}S^{2}}(u_{j}-2u_{j-1/2}+u_{j-1}).$
(65)

An example of Greeks computation results is shown in Figure 6 for $\Delta$ and
$\Gamma$ during the first two years ($t\in[0,2]$). The corresponding Greeks
are computed using P2-FEM for $\Delta$ and $\Gamma$ at $S_{j-1/2}$.

Visible peaks and spikes in the solution surfaces are due to the constraints
and coupon payment, at which point the solution is not differentiable; Delta
and Gamma, however, are smooth at the locations where the inequality
constraints are not imposed.

Figure 6: Greeks by with the parameters $t=[0,2]$, $r=0.05$, $r_{c}=0.02$,
$\sigma=0.2$, $F=\$100$, $K=\$4$, $\rho=10^{12}$, $n_{E}=100$, $n_{t}=600$ by
P2-FEM; Left figure: Delta $\Delta$; Right figure: Gamma $\Gamma$.

###  5.3 Accuracy of the approximations

As a way to measure the accuracy of our FEM approach, we consider a linear
model problem [25] without constraints (and hence no penalty term) and coupon
payment corresponding to the exact solution

$\displaystyle U(x,\tau)$
$\displaystyle=S_{int}^{2}e^{2x}\sqrt{S_{int}e^{x}}-Fe^{-r\tau}\sqrt{S_{int}e^{x}},$
$\displaystyle V(x,\tau)$
$\displaystyle=S_{int}^{2}e^{2x}\sqrt{S_{int}e^{x}}-Fe^{-r\tau}\sqrt{S_{int}e^{x}}+x^{2}\tau,$

where $\tau\in(0,1)$. This manufactured solution corresponds to the initial
conditions

$\displaystyle\begin{cases}\displaystyle
U(x,0)=S_{int}^{2}e^{2x}\sqrt{S_{int}e^{x}}-F\sqrt{S_{int}e^{x}},\\\
\displaystyle
V(x,0)=S_{int}^{2}e^{2x}\sqrt{S_{int}e^{x}}-F\sqrt{S_{int}e^{x}},\end{cases}$

the boundary conditions

$\displaystyle\begin{cases}\displaystyle
U(0,\tau)=S_{int}^{2}\sqrt{S_{int}}-Fe^{-r\tau}\sqrt{S_{int}},\\\
\displaystyle
V(0,\tau)=S_{int}^{2}\sqrt{S_{int}}-Fe^{-r\tau}\sqrt{S_{int}},\end{cases}\begin{cases}\displaystyle
U(1,\tau)=S_{int}^{2}e^{2}\sqrt{S_{int}e}-Fe^{-r\tau}\sqrt{S_{int}e},\\\
\displaystyle
V(1,\tau)=S_{int}^{2}e^{2}\sqrt{S_{int}e}-Fe^{-r\tau}\sqrt{S_{int}e}+\tau,\end{cases}$

and additional forcing terms

$\displaystyle\displaystyle f_{1}$
$\displaystyle=U_{\tau}-\frac{\sigma^{2}}{2}U_{xx}-(r-\frac{\sigma^{2}}{2})U_{x}+rU+r_{c}V,$
$\displaystyle\displaystyle f_{2}$
$\displaystyle=V_{\tau}-\frac{\sigma^{2}}{2}V_{xx}-(r-\frac{\sigma^{2}}{2})V_{x}+(r+r_{c})V,$

in the CB and COCB PDE, respectively. Applying FEM and the $\theta$-scheme
leads to the numerical procedure

$\displaystyle A_{11}\boldsymbol{u}^{m+1}$
$\displaystyle=\widetilde{A}_{11}\boldsymbol{u}^{m}-A_{12}\boldsymbol{v}^{m+1}+\widetilde{A}_{12}\boldsymbol{v}^{m}+\theta\Delta\tau\boldsymbol{\beta}_{1}^{m+1}+(1-\theta)\Delta\tau\boldsymbol{\beta}_{1}^{m}+\hat{\boldsymbol{b}}^{m}_{M,u}-\hat{\boldsymbol{b}}^{m+1}_{M,u}+\theta\Delta\tau
f_{1}+(1-\theta)\Delta\tau f_{1},$ $\displaystyle A_{22}\boldsymbol{v}^{m+1}$
$\displaystyle=\widetilde{A}_{22}\boldsymbol{v}^{m}+\theta\Delta\tau\boldsymbol{\beta}_{2}^{m+1}+(1-\theta)\Delta\tau\boldsymbol{\beta}_{2}^{m}+\hat{\boldsymbol{b}}^{m}_{M,v}-\hat{\boldsymbol{b}}^{m+1}_{M,v}+\theta\Delta\tau
f_{2}+(1-\theta)\Delta\tau f_{2}.$

Calculated errors are presented Figure 7 and 8 using two measures [3, 2]:

$\displaystyle\|Error\|_{L^{2}}$
$\displaystyle=\|U(x,1)-\boldsymbol{u}^{n_{t}}\|_{L^{2}},$
$\displaystyle\|Error\|_{L^{\infty}(L^{2})}$ $\displaystyle=\max_{1\leq m\leq
n_{t}}\left(\|U(x,m\Delta\tau)-\boldsymbol{u}^{m}\|_{L^{2}}\right),$

where $\boldsymbol{u}^{m}$ is the solution of the model problem at
$\tau=m\Delta\tau$, computed by P1-FEM or P2-FEM. In Figure 7, the errors are
calculated for varying $\Delta\tau$ and a fixed value of $h$. The errors
decreases as $\Delta\tau$ is reduced to 0, with a rate that is proportional to
$\Delta\tau$ (first-order convergence). This first order convergence is
consistent with the result in [6].

Figure 7: Error estimates of the model problem in $\tau\in[0,1]$, with
$r=0.05$,$r_{c}=0.02$, $\sigma=0.2$, $h_{P1}=3\times 10^{-4}$, and
$h_{P2}=10^{-3}$.

In Figure 8, the errors are calculated for varying $h$ and fixed $\Delta\tau$.
The plots suggest convergence of P1-FEM and P2-FEM at the rate proportional to
$h$ and $h^{2}$, respectively (Theoretically, e.g. for P1-FEM the convergence
rate is given by the relation $\|Error\|_{L^{\infty}(L^{2})}\leq
C(h+\Delta\tau)$). This convergence rate is as expected for the linear model
used but cannot, however, be expected when nonlinear (e.g., penalty) terms are
added.

Figure 8: Error estimates of the MMS model in $x\in[0,1]$, with
$r=0.05$,$r_{c}=0.02$, $\sigma=0.2$, and $\Delta\tau=10^{-4}$.

##  6 Conclusion

We presented numerical methods based on P1 and P2 finite element methods for
the TF model in convertible bond pricing, reformulated as a penalty PDE. The
resultant differential algebraic equations were solved using $\theta$-scheme,
with nonlinearity due to the penalty handled by Newton’s method.

We reported results from the numerical simulations using the two finite
element approximations and compared them with finite difference methods.
Comparing results from the two methods, we observed extremely comparable
results. Furthermore, we verified the advantage of using P2-FEM over P1-FEM,
when it comes to accuracy of the solutions using fewer elements.

The finite-element method described in this paper can be implemented on
different models for convertible bond pricing, such as the AFV model [6].
Results on this will be reported in the future.

## References

  * [1] Y. H. AlSaedi and G. A. Tularam.  A review of the recent advances made in the Black-Scholes models and respective solutions methods.  Journal of Mathematics and Statistics, 14(1):29–39, January 2018\.
  * [2] Y. Amanbek, G. Singh, G. Pencheva, and M. F. Wheeler.  Error indicators for incompressible Darcy flow problems using enhanced velocity mixed finite element method.  Computer Methods in Applied Mechanics and Engineering, 363:112884, 2020.
  * [3] Y. Amanbek and M. F. Wheeler.  A priori error analysis for transient problems using enhanced velocity approach in the discrete-time setting.  Journal of Computational and Applied Mathematics, 361:459–471, 2019\.
  * [4] M. Ammann, A. Kind, and C. Wilde.  Simulation-based pricing of convertible bonds.  Journal of Empirical Finance, 15(2):310–331, 2008.
  * [5] J. Ankudinova and M. Ehrhardt.  On the numerical solution of nonlinear Black-Scholes equations.  Computers and Mathematics with Applications, 56(3):799–812, August 2008.
  * [6] E. Ayache, P. A. Forsyth, and K. R. Vetzal.  Valuation of convertible bonds with credit risk.  The Journal of Derivatives, 11(1):9–29, 2003.
  * [7] G. Barone-Adesi, A. Bermudez, and J. Hatgioannides.  Two-factor convertible bonds valuation using the method of characteristics/finite elements.  Journal of Economic Dynamics and Control, 27(10):1801–1831, August 2003.
  * [8] F. Black and M. S. Scholes.  The pricing of options and corporate liabilities.  Journal of Political Economy, 81:637–654, 1973.
  * [9] M. J. Brennan and E. S. Schwartz.  Convertible bonds: Valuation and optimal strategies for call and conversion.  The Journal of Finance, 32:1699–1715, 1977.
  * [10] M. J. Brennan and E. S. Schwartz.  Analyzing convertible bonds.  Journal of Financial and Quantitative Analysis, 15:907–929, 1980\.
  * [11] D. Černá.  Numerical solution of the Black-Scholes equation using cubic spline wavelets.  In G. Venkov et al., editor, AIP Conference Proceedings 1789, page 030001. AIP Publishing, 2016.
  * [12] C. C. Christara and R. Wu.  Penalty and penalty-like methods for nonlinear HJB PDEs.  Applied Mathematics and Computation, 425:127015, 2022.
  * [13] J. de Frutos.  A finite element method for two factor convertible bonds.  In Numerical Methods in Finance, pages 109–128. Springer US, 2005\.
  * [14] E. Dremkova and M. Ehrhardt.  A high-order compact method for nonlinear Black-Scholes option pricing equations of American options.  International Journal of Computer Mathematics, 88(13):2782–2797, 2011.
  * [15] C. A. J. Fletcher.  The group finite element formulation.  Computer Methods in Applied Mechanics and Engineering, 37(2):225–244, 1983.
  * [16] P. A. Forsyth and K. R. Vetzal.  Quadratic convergence for valuing American options using a penalty method.  SIAM Journal on Scientific Computing, 23(6):2095–2122, January 2002.
  * [17] P. A. Forsyth, K. R. Vetzal, and R. Zvan.  A finite element approach to the pricing of discrete lookbacks with stochastic volatility.  Applied Mathematical Finance, 6(2):87–106, 1999.
  * [18] J. C. Hull.  Options, Futures, and Other Derivatives.  Pearson Education India, 2003.
  * [19] J. E. Ingersoll Jr.  A contingent-claims valuation of convertible securities.  Journal of Financial Economics, 4:289–321, May 1977.
  * [20] P. Kovalov and V. Linetsky.  Valuing convertible bonds with stock price, volatility, interest rate, and default risk.  SSRN Electronic Journal, 2008.
  * [21] P. K. Kythe and D. Wei.  An Introduction to Linear and Nonlinear Finite Element Analysis: A Computational Approach, volume 57.  The American Society of Mechanical Engineers, 2004.
  * [22] S. Lin and S.-P. Zhu.  Numerically pricing convertible bonds under stochastic volatility or stochastic interest rate with an ADI-based predictor-corrector scheme.  Computers and Mathematics with Applications, 79(5):1393–1419, 2020\.
  * [23] K. Milanov and O. Kounchev.  Binomial tree model for convertible bond pricing within equity to credit risk framework.  SSRN Electronic Journal, 2012.
  * [24] R. Rannacher.  Finite element solution of diffusion problems with irregular data.  Numerische Mathematik, 43(2):309–327, 1984.
  * [25] P. J. Roache.  Code verification by the method of manufactured solutions.  Journal of Fluids Engineering, 124:4–10, 2001.
  * [26] J. De Spiegeleer, W. Schoutens, and C. Van Hulle.  The Handbook of Hybrid Securities: Convertible Bonds, CoCo Bonds and Bail-In.  The Wiley Finance Series. Wiley, 1st edition, 2014.
  * [27] J. De Spiegeleer, W. Schoutens, and P. Jabre.  The Handbook of Convertible Bonds: Pricing, Strategies and Risk Management.  The Wiley Finance Series. Wiley, 1st edition, 2011.
  * [28] K. Tsiveriotis and C. Fernandes.  Valuing convertible bonds with credit risk.  The Journal of Fixed Income, 8:95–102, September 1998.


* Abdo  _et al._ [2010b] A. A. Abdo, M. Ackermann, M. Ajello, W. B. Atwood, L. Baldini, J. Ballet, G. Barbiellini, D. Bastieri, K. Bechtol, R. Bellazzini, B. Berenji, E. D. Bloom, E. Bonamente, A. W. Borgland, J. Bregeon, A. Brez, M. Brigida, P. Bruel, T. H. Burnett, S. Buson, G. A. Caliandro, R. A. Cameron, P. A. Caraveo, J. M. Casandjian, C. Cecchi, A. Chekhtman, C. C. Cheung, J. Chiang, S. Ciprini, R. Claus, J. Cohen-Tanugi, J. Conrad, A. de Angelis, F. de Palma, S. W. Digel, E. d. C. e. Silva, P. S. Drell, A. Drlica-Wagner, R. Dubois, D. Dumora, C. Farnier, C. Favuzzi, S. J. Fegan, W. B. Focke, P. Fortin, M. Frailis, Y. Fukazawa, P. Fusco, F. Gargano, N. Gehrels, S. Germani, B. Giebels, N. Giglietto, F. Giordano, T. Glanzman, G. Godfrey, I. A. Grenier, J. E. Grove, L. Guillemot, S. Guiriec, M. Gustafsson, A. K. Harding, E. Hays, D. Horan, R. E. Hughes, M. S. Jackson, T. E. Jeltema, G. Jóhannesson, A. S. Johnson, R. P. Johnson, W. N. Johnson, T. Kamae, H. Katagiri, J. Kataoka, M. Kerr, J. Knödlseder, M. Kuss, J. Lande, L. Latronico, M. Lemoine-Goumard, F. Longo, F. Loparco, B. Lott, M. N. Lovellette, P. Lubrano, G. M. Madejski, A. Makeev, M. N. Mazziotta, J. E. McEnery, C. Meurer, P. F. Michelson, W. Mitthumsiri, T. Mizuno, A. A. Moiseev, C. Monte, M. E. Monzani, E. Moretti, A. Morselli, I. V. Moskalenko, S. Murgia, P. L. Nolan, J. P. Norris, E. Nuss, T. Ohsugi, N. Omodei, E. Orlando, J. F. Ormes, D. Paneque, J. H. Panetta, D. Parent, V. Pelassa, M. Pepe, M. Pesce-Rollins, F. Piron, T. A. Porter, S. Profumo, S. Rainò, R. Rando, M. Razzano, A. Reimer, O. Reimer, T. Reposeur, S. Ritz, A. Y. Rodriguez, M. Roth, H. F. W. Sadrozinski, A. Sander, P. M. Saz Parkinson, J. D. Scargle, T. L. Schalk, A. Sellerholm, C. Sgrò, E. J. Siskind, D. A. Smith, P. D. Smith, G. Spandre, P. Spinelli, M. S. Strickman, D. J. Suson, H. Takahashi, T. Takahashi, T. Tanaka, J. B. Thayer, J. G. Thayer, D. J. Thompson, L. Tibaldo, D. F. Torres, A. Tramacere, Y. Uchiyama, T. L. Usher, V. Vasileiou, N. Vilchez, V. Vitale, A. P. Waite, P. Wang, B. L. Winer, K. S. Wood, T. Ylinen, M. Ziegler, J. S. Bullock, M. Kaplinghat, G. D. Martinez, and Fermi LAT Collaboration, Observations of Milky Way Dwarf Spheroidal Galaxies with the Fermi-Large Area Telescope Detector and Constraints on Dark Matter Models, ApJ 712, 147 (2010b), arXiv:1001.4531 [astro-ph.CO] .
  * Ando and Nagai [2012] S. Ando and D. Nagai, Fermi-LAT constraints on dark matter annihilation cross section from observations of the Fornax cluster, J. Cosmology Astropart. Phys 2012, 017 (2012), arXiv:1201.0753 [astro-ph.HE] .
  * Steigman  _et al._ [2012] G. Steigman, B. Dasgupta, and J. F. Beacom, Precise relic WIMP abundance and its impact on searches for dark matter annihilation, Phys. Rev. D 86, 023506 (2012), arXiv:1204.3622 [hep-ph] .
  * Levi  _et al._ [2019] M. Levi, L. E. Allen, A. Raichoor, C. Baltay, S. BenZvi, F. Beutler, A. Bolton, F. J. Castander, C.-H. Chuang, A. Cooper, J.-G. Cuby, A. Dey, D. Eisenstein, X. Fan, B. Flaugher, C. Frenk, A. X. Gonzalez-Morales, O. Graur, J. Guy, S. Habib, K. Honscheid, S. Juneau, J.-P. Kneib, O. Lahav, D. Lang, A. Leauthaud, B. Lusso, A. de la Macorra, M. Manera, P. Martini, S. Mao, J. A. Newman, N. Palanque-Delabrouille, W. J. Percival, C. Allende Prieto, C. M. Rockosi, V. Ruhlmann-Kleider, D. Schlegel, H.-J. Seo, Y.-S. Song, G. Tarle, R. Wechsler, D. Weinberg, C. Yeche, and Y. Zu, The Dark Energy Spectroscopic Instrument (DESI), in  _Bulletin of the American Astronomical Society_, Vol. 51 (2019) p. 57, arXiv:1907.10688 [astro-ph.IM] .
  * LSST Dark Energy Science Collaboration [2012] LSST Dark Energy Science Collaboration, Large Synoptic Survey Telescope: Dark Energy Science Collaboration, arXiv e-prints , arXiv:1211.0310 (2012), arXiv:1211.0310 [astro-ph.CO] .
  * Refregier  _et al._ [2010] A. Refregier, A. Amara, T. D. Kitching, A. Rassat, R. Scaramella, and J. Weller, Euclid Imaging Consortium Science Book, arXiv e-prints , arXiv:1001.0061 (2010), arXiv:1001.0061 [astro-ph.IM] .
  * Spergel  _et al._ [2015] D. Spergel, N. Gehrels, C. Baltay, D. Bennett, J. Breckinridge, M. Donahue, A. Dressler, B. S. Gaudi, T. Greene, O. Guyon, C. Hirata, J. Kalirai, N. J. Kasdin, B. Macintosh, W. Moos, S. Perlmutter, M. Postman, B. Rauscher, J. Rhodes, Y. Wang, D. Weinberg, D. Benford, M. Hudson, W. S. Jeong, Y. Mellier, W. Traub, T. Yamada, P. Capak, J. Colbert, D. Masters, M. Penny, D. Savransky, D. Stern, N. Zimmerman, R. Barry, L. Bartusek, K. Carpenter, E. Cheng, D. Content, F. Dekens, R. Demers, K. Grady, C. Jackson, G. Kuan, J. Kruk, M. Melton, B. Nemati, B. Parvin, I. Poberezhskiy, C. Peddie, J. Ruffa, J. K. Wallace, A. Whipple, E. Wollack, and F. Zhao, Wide-Field InfrarRed Survey Telescope-Astrophysics Focused Telescope Assets WFIRST-AFTA 2015 Report, arXiv e-prints , arXiv:1503.03757 (2015), arXiv:1503.03757 [astro-ph.IM] .
  * Berlind and Weinberg [2002] A. A. Berlind and D. H. Weinberg, The Halo Occupation Distribution: Toward an Empirical Determination of the Relation between Galaxies and Mass, ApJ 575, 587 (2002), arXiv:astro-ph/0109001 [astro-ph] .
  * van den Bosch  _et al._ [2013] F. C. van den Bosch, S. More, M. Cacciato, H. Mo, and X. Yang, Cosmological constraints from a combination of galaxy clustering and lensing - I. Theoretical framework, MNRAS 430, 725 (2013), arXiv:1206.6890 [astro-ph.CO] .

##  Appendix A Halo models for the UGRB and galaxies

###  A.1 The mean dark matter profile

To describe the mean density of dark matter around haloes, we will use the
Navarro-Frenk and White profile (NFW, [59]), which takes the form

$\rho_{\rm
NFW}(r)=\frac{\rho_{0}}{x\left(1+x\right)^{2}}\Theta(r<r_{\Delta}),\hskip
12.0ptx\equiv r/r_{s}.$ (64)

Here the comoving spherical overdensity radius $r_{\Delta}$ is related to the
halo mass via

$M=\frac{4\pi}{3}\Delta\,\bar{\rho}_{c}r_{\Delta}^{3},$ (65)

where $\bar{\rho}_{c}$ is the critical density, and $\Delta$ is the so-called
spherical-overdensity parameter, which we will set to $\Delta=200$. The
characteristic radius $r_{s}$ is related to $r_{\Delta}$ through the
concentration-mass relation $c(M)$

$r_{\Delta}=c(M)r_{s}.$ (66)

We will use the parametrisation of [99] to calculate $c(M)$. The normalisation
$\rho_{0}$ can be found by integrating the density profile over volume, and is
related to the halo mass via

$\rho_{0}=\frac{M}{4\pi\,r_{s}^{3}[\log(1+c)-c/(1+c)]}.$ (67)

The simple form of the NFW profile makes it possible to compute its Fourier
transform analytically:

$\displaystyle\rho_{\rm NFW}(k|M)$ $\displaystyle\equiv
4\pi\int_{0}^{\infty}dr\,r^{2}\,\rho_{\rm NFW}(r|M)\,\frac{\sin kr}{kr}$
$\displaystyle=\frac{M}{\log(1+c)-c/(1+c)}\,\left[\cos q\left({\rm
Ci}((1+c)q)-{\rm Ci}(q)\right)+\sin q\left({\rm Si}((1+c)q)-{\rm
Si}(q)\right)-\frac{\sin cq}{1+cq}\right],$ (68)

with $q\equiv kr_{\Delta}/c$.

###  A.2 The halo occupation distribution

As an alternative model describing the non-linear relation between the
galaxies and matter overdensities, we use the halo occupation distribution
framework (HOD, [109, 110]). In this case, the number density of galaxies in a
halo of mass $M$ is parametrised in terms of the mean number of central and
satellite galaxies ($\bar{N}_{c}(M)$, and $\bar{N}_{s}(M)$ respectively), as
well as the profile describing the distribution of satellites around a given
halo $u_{s}(r|M)$. The mean overdensity of galaxies around a halo of mass $M$
is thus given by

$\langle 1+\delta_{g}({\bf
x}|M)\rangle=\frac{\bar{N}_{c}(M)}{\bar{n}_{g}}\left[\delta^{D}({\bf
x})+\bar{N}_{s}(M)\,u_{s}(|{\bf x}||M)\right],$ (69)

where the mean number density of galaxies $\bar{n}_{g}$ is given by

$\bar{n}_{g}=\int dM\,n(M)\,\bar{N}_{c}(M)\left[1+\bar{N}_{s}(M)\right].$ (70)

In this work, we will use the same parametrisation employed in [92, 78]. The
mean number of centrals and satellites is given by:

$\displaystyle\bar{N}_{c}(M)=\frac{1}{2}\left[1+{\rm
erf}\left(\frac{\log(M/M_{\rm min})}{\sigma_{\log M}}\right)\right],\hskip
12.0pt\bar{N}_{s}(M)=\Theta(M-M_{0})\left(\frac{M-M_{0}}{M_{1}}\right)^{\alpha}.$
(71)

Where $M_{\rm min}$ is the mass at which haloes have on average $0.5$ central
galaxies, $M_{0}$ is the mass at which haloes may start forming satellite
galaxies, and $M_{1}$ is the typical mass for haloes hosting satellites. We
further assume that satellites are distributed around haloes following the
matter density, and hence

$u_{s}(r|M)\equiv\rho_{\rm NFW}(r|M)/M.$ (72)

We will fix the free parameters of the model to the best-fit values for the
2MPZ and WI-SC samples found in [78]. In this case, $M_{0}=M_{\rm min}$,
$\sigma_{\log M}=0$, and $\alpha=1$ in all cases. The remaining free
parameters ($M_{\rm min},M_{1}$) are fitted in each redshift bin, and take the
values found in Table 1.

###  A.3 Halo profile for DM annihilation

As noted in Section II.4.3, we consider two limits regarding the impact of
substructure on DM annihilation. In the most conservative limit, we study the
case with no substructure. In this case, all haloes follow an exact NFW
profile, and $\langle\rho_{\rm DM}^{2}\rangle=\langle\rho_{\rm
DM}\rangle^{2}$. Thus, the halo profile for annihilation without substructure
is

$U^{\rm ann,no-sub}_{\gamma}(r|M)=\frac{\rho_{\rm
NFW}^{2}(r|M)}{\bar{\rho}_{M,0}^{2}}.$ (73)

The simplicity of the NFW functional form makes it possible to calculate the
Fourier transform of the squared profile analytically:

$U^{\rm ann,no-sub}_{\gamma}(k|M)={\cal
V}(M)\left[T_{1}(k|M)+T_{2}(k|M)+T_{3}(k|M)+T_{4}(k|M)\right],$ (74)

with

$\displaystyle T_{1}(k|M)\equiv\frac{1}{6q}\left[q(q^{2}-6)\sin
q+3(q^{2}-2)\cos q\right]\,\left[{\rm Si}((1+c)q)-{\rm Si}(q)\right],$ (75)
$\displaystyle T_{2}(k|M)\equiv\frac{1}{6q}\left[q(q^{2}-6)\cos
q-3(q^{2}-2)\sin q\right]\,\left[{\rm Ci}((1+c)q)-{\rm Ci}(q)\right],$ (76)
$\displaystyle
T_{3}(k|M)\equiv-\frac{c((q^{2}-6)c+2q^{2}-15)+q^{2}-11}{6q(1+c)^{3}}\sin(cq),$
(77) $\displaystyle
T_{4}(k|M)\equiv\frac{(3c+4)\cos(qc)}{6(1+c)^{2}}+\frac{{\rm
Si}(qc)}{q}-\frac{2}{3},$ (78) $\displaystyle{\cal V}(M)\equiv\frac{4\pi
r_{s}^{3}\rho_{0}^{2}}{\bar{\rho}_{M,0}^{2}},$ (79)

where $\rho_{0}$ is given in Eq. 67, and $q\equiv kr_{\Delta}/c$ as before.

To describe the effects of substructure, we extend this model as outlined in
[64]. In this case, the contribution from subhaloes is characterised by a
mass-dependent boost factor $b_{\rm sh}(M)$. The associated halo profile is
then given by

$U^{\rm ann}_{\gamma}(k|M)=U^{\rm ann,no-sub}_{\gamma}(k|M)+b_{\rm
sh}(M)\,V_{\rm host}(M)\,u_{\rm sh}(k|M),$ (80)

where $U^{\rm ann,no-sub}_{\gamma}(k|M)$, given in Eq. 74, is the smooth
contribution from the squared host halo density, and the second term is the
contribution from the squared density of subhaloes. This neglects a term of
the form ${\rm host}\times{\rm subhalo}$, which is only relevant on scales
where the densities of host halo and subhaloes are similar. In Eq. 80, $V_{\rm
host}$ is simply the volume integral of the host profile,

$V_{\rm host}(M)\equiv 4\pi\int dr\,r^{2}\,U^{\rm ann,no-
sub}_{\gamma}(r|M)={\cal V}(M)\left[1-\frac{1}{(1+c)^{3}}\right],$ (81)

and $u_{\rm sh}$ is the squared density profile of subhaloes normalised to
have a unit volume integral. For this, we use the functional form of [64]
which, in Fourier space, is

$\displaystyle u_{\rm
sh}(k|M)=A\left[\int_{0}^{1}dx\frac{x^{2}}{(x^{2}+1/16)^{3/2}}\frac{\sin{\kappa
x}}{\kappa
x}+\frac{64}{17^{3/2}}\frac{\kappa\cos{\kappa}+\eta\sin{\kappa}}{\kappa(\kappa^{2}+\eta^{2})}\right],$
(82)

where $A=0.636$ ensures that $u_{\rm sh}(k\rightarrow 0|M)\rightarrow 1$, and
$\kappa\equiv kr_{\Delta}$.

Figure 12: Boost factor as a function of host halo mass for with the three
models of substructure used in this analysis, based on [63] (solid black),
[62] (dashed blue), and [61] (dotted orange).

The amplitude of the contribution from substructures is controlled by the
boost factor $b_{\rm sh}(M)$. We consider three different models of
substructure. Ordering them by the amplitude of the effective boost provided
to the annihilation signal, these are:

  * •

SC14. We use the fitting function of [61]:

$\log_{10}b_{\rm sh}(M)=\sum_{n=0}^{5}b_{n}\log(M)^{n},$ (83)

with $b_{n}=\\{-0.442,0.0796,-0.0025,4.77\times 10^{-6},4.77\times
10^{-6},-9.69\times 10^{-8}\\}$.

  * •

M16. We use the fitting fucnction of [62]:

$\log_{10}b_{\rm sh}(M)=\sum_{n=0}^{5}b_{n}\log_{10}(M)^{n},$ (84)

with $b_{n}=\\{-0.186,0.144,-8.8\times 10^{-3},1.13\times 10^{-3},-3.7\times
10^{-5},-2\times 10^{-7}\\}$.

  * •

G12. We use the fitting function of [63]:

$b_{\rm sh}=110\left(\frac{M}{10^{12}M_{\odot}}\right)^{0.39}.$ (85)

The boost factors associated with each of these models are shown in Fig. 12.

##  Appendix B Minimum halo mass for annihilation

One of the most important sources of theoretical uncertainties when deriving
constraints from DM annihilation, is the minimum halo mass contributing to the
signal. In this work, we have fixed this to $M_{\rm
min}=10^{-6}\,\,M_{\odot}$, following most other analyses in the literature
(e.g. [25]). Fig. 13 shows, however, that the choice of $M_{\rm min}$ can have
a significant impact on the $\langle\sigma v\rangle$ constraints, and on their
interpretation.

The left panel shows the 95%-level constraints on $\langle\sigma v\rangle$
obtained assuming different values for $M_{\rm min}$, in the range
$[10^{-6}\,M_{\odot},10^{6}\,M_{\odot}]$. The constraints may vary by up to
one order of magnitude depending on this choice. The value of $M_{\rm min}$
also affects the redshift evolution of the signal, and thus can affect the
interpretation of tomographic measurements. The right panel shows the redshift
evolution parameter $\Delta\beta_{\rm ann}$, defined in Eq. 63, as a function
of minimum halo mass. Depending on the model of substructure used,
$\Delta\beta$ may take values spanning the range
$-2.8\lesssim\Delta\beta\lesssim-1.5$.

Figure 13: Left panel: 95% upper bounds on the DM annihilation cross section
for the M16 substructure case. The increasing darkness of green corresponds to
increasing the minimum mass $M_{\rm min}$. The lightest shade corresponds to
our fiducial constraints presented in Section V.2.2. Right panel: the redshift
evolution parameter $\Delta\beta_{\rm ann}$ (see Eq. 63 as a function of
$M_{\rm min}$ for the M16 and G12 substructure models (red and blue bands,
respectively).