# On Merton’s Optimal Portfolio Problem
under Sporadic Bankruptcy

Yaacov Kopeliovich and Michael Pokojovy  University of Connecticut, Finance
Department, 2100 Hillside Road, Storrs, CT 06269 and Old Dominion University,
Department of Mathematics and Statistics, 4300 Elkhorn Ave, Norfolk, VA 23529
<EMAIL_ADDRESS>
###### Abstract.

Consider a stock market following a geometric Brownian motion and a riskless
asset continuously compounded at a constant rate. Assuming the stock can go
bankrupt, i.e., lose all of its value, at some exogenous random time
(independent of the stock price) modeled as the first arrival time of a
homogeneous Poisson process, we study the Merton’s optimal portfolio problem
consisting of maximizing the expected logarithmic utility of the total wealth
at a preselected finite maturity time. First, we present a heuristic
derivation based on a new type of Hamilton-Jacobi-Bellman equation. Then, we
formally reduce the problem to a classical controlled Markovian diffusion with
a new type of terminal and running costs. A new version of Merton’s ratio is
rigorously derived using Bellman’s dynamic programming principle and validated
with a suitable type of verification theorem. A real-world example comparing
the latter ratio to the classical Merton’s ratio is given.
Keywords: Merton’s portfolio problem; Markov-switching processes; absorbing
processes; Hamilton-Jacobi-Bellman equation; logarithmic utility.

Helpful discussions with Professor Ramon van Handel on the verification
theorem for classical Merton problem with logarithmic utility are greatly
appreciated.

##  1\. Introduction

In this paper, we investigate the dynamic optimal allocation problem in a
situation where security prices are subject to exogenous events like
bankruptcy. The optimal portfolio problem was first considered by R. Merton
[25] for multiple stocks that follow a multi-dimensional geometric Brownian
motion [27]. For a single stock, the price $S_{t}$ at time $t$ is modeled as a
Markovian diffusion process that solves the stochastic differential equation

(1) $\mathrm{d}S_{t}=\mu S_{t}\mathrm{d}t+\sigma S_{t}\mathrm{d}B_{t}$

with a standard Wiener process $(B_{t})_{t\geq 0}$, while the risk-free asset
follows the ordinary differential equation

(2) $\mathrm{d}R_{t}=rR_{t}$

with a constant risk-free rate $r$. Assuming an isoelastic utility function

(3) $U(W)=\begin{cases}\frac{W^{1-\gamma}-1}{1-\gamma},&\gamma\geq 0,\\\
\log(W),&\gamma=1\end{cases}$

for some $\gamma\geq 0$, he derived [25] his, now classical, Merton ratio
$\pi$ of the wealth to invest in stock

(4) $\pi=\frac{\mu-r}{\gamma\sigma^{2}},$

while the remaining fraction, $(1-\pi)$, of investor’s wealth at each time $t$
is invested in the risk-free asset.

The theoretical importance of Merton ratio is that it gives investors a simple
and easily interpretable, yet rigorous rule how to allocate capital into the
risky asset as to maximize the expected utility. The practical importance of
this ratio is highlighted in a recent The Economist article [30]. We note that
a similar formula exists that provides a guideline how to invest in multiple
stocks. Indeed, assume the price vector
$\boldsymbol{S}_{t}=(S_{t}^{1},S_{t}^{2},\dots,S_{t}^{p})^{T}$ of $p$ stocks
follows a multivariate geometric Brownian motion

(5)
$\mathrm{d}\boldsymbol{S}_{t}=\operatorname{diag}(\boldsymbol{S})\big{(}\boldsymbol{\mu}\,\mathrm{d}t+\boldsymbol{\Sigma}^{1/2}\mathrm{d}\boldsymbol{B}_{t}\big{)}$

with a $p$-variate standard Brownian motion $(\boldsymbol{B}_{t})_{t\geq 0}$,
where $\boldsymbol{\mu}\in\mathbb{R}^{p}$ is the vector of expected returns
and $\boldsymbol{\Sigma}\in\mathbb{R}^{p\times p}$ is a symmetric, positive
definite covariance matrix between the $p$-variate innovations. According to
[25], the optimal weights
$\boldsymbol{\pi}=(\pi^{1},\pi^{2},\dots,\pi^{p})^{T}$ are given by

(6)
$\boldsymbol{\pi}=\frac{1}{\gamma}\boldsymbol{\Sigma}^{-1}\big{(}\boldsymbol{\mu}-r\mathds{1}_{p}\big{)}\text{
with }\mathds{1}_{p}=(1,1,\dots,1)\in\mathbb{R}^{p}.$

Since the groundbreaking work of Merton [25], many authors investigated
various variants of Merton optimal portfolio problem. For example, Wachter
[31] considered the optimal portfolio problem for risky assets following a
mean reverting processes and obtained an exact solution for the optimal
allocation when markets are complete. Merton himself [26] applied his dynamic
portfolio optimization approach to extend the usual capital asset pricing
model (CAPM) to an in-temporal CAPM. The latter generalizes the usual CAPM
model developed by Sharpe based on Markowitz portfolio problem. Efforts to
incorporate other types of risky assets (such as fixed income securities) have
been attempted as well. For instance, Brennan et al. [6] considered a Merton-
type problem for fixed income securities and stocks but without any bankruptcy
assumptions on the bonds. Karatzas et al. [16] considered the optimal
portfolio problem for a general utility function and incorporated bankruptcy
events. The bankruptcy event considered in the latter paper was not modeled as
an exogenous event but was rather assumed to arise from the self-exciting
(autoregressive) dynamics, i.e., the authors asked the question when it will
be optimal for the security to default given an infinite time horizon and
solved this problem of optimal bankruptcy. For applications in corporate
finance, see [21] and [22].

The goal of this paper is to revisit Merton’s original problem and investigate
the optimal allocation between stock and cash as to optimize the expected
logarithmic utility at a finite maturity $T$ assuming, however, the stock can
sporadically jump into bankruptcy at an exogenous random time
$\tau\sim\mathrm{Exp}(\lambda)$, independent of the stock price, which is
modeled as an exponential random variable with some parameter $\lambda>0$
corresponding to default intensity. The amended price dynamics reads then as

(7) $\mathrm{d}S_{t}=\left\\{\begin{array}[]{cl}\mu S_{t-}\mathrm{d}t+\sigma
S_{t-}\mathrm{d}B_{t}&\text{for }t<\tau\\\ 0\,\mathrm{d}t,&\text{for
}t>\tau\end{array}\right.\text{ and }\Delta S_{\tau}=-S_{\tau-}$

with the jump operator (viz. [23, Equation (1)])

(8) $\Delta S_{t}:=S_{t}-S_{t-}$

where $(B_{t})_{t\geq 0}$ is a standard Wiener process. The process in
Equation (7) is a special two-state case of a more general Markov-switching
jump diffusion process considered in [4, 32]. Indeed, the random bankruptcy
time $\tau$ can be viewed as the first arrival time of a Poisson process or
the jump time of a two-state continuous time Markov process with two states:
pre-bankruptcy and post-bankruptcy, where the latter state is absorbing. The
average time to bankruptcy $\mathbb{E}[\tau]$ is $\frac{1}{\lambda}$, while
the conditional probability of banktruptcy is memoryless

(9) $\mathbb{P}\big{\\{}t\leq\tau<t+\Delta t\,|\,\tau\geq
t\big{\\}}=\lambda\Delta t+o(\Delta t)\text{ as }\Delta t\searrow 0\text{ for
any }t\geq 0$

as it is known to be the case for Poisson processes. A stochastic dynamics
similar to Equation (7), albeit with non-random jump times, has recently been
studied in [23]. Assuming the logarithmic utility function in this paper, we
obtain the Merton ratio

(10)
$\pi_{t}=\left\\{\begin{array}[]{cl}\frac{\mu-r+\sigma^{2}-\sqrt{\left(\sigma^{2}-(\mu-r)\right)^{2}+4\lambda\sigma^{2}}}{2\sigma^{2}},&t<\tau,\\\
0,&t\geq 0.\end{array}\right.$

The expression in Equation (10) is reminiscent of the multivariate optimal
allocation vector obtained in [1] in the context of portfolio choice in
markets with contagion.

To some extent, our work has connection to portfolio optimization under
Brownian motion perturbed with jump processes. For example, [9] investigated
optimal portfolios under Levý innovations. Their optimization approach
followed the classical Markowitz paradigm of mean variance. The authors
adopted a framework of risk constraints, which is related to our problem as we
consider stocks that have a risk of sporadic bankruptcy. An important
distinction that jumps under a Levý process have an emerging property while
stocks (contrary to bonds) do not emerge from bankruptcy.

A dynamic portfolio optimization problem for a stochastic model with jumps was
considered in [10]. Like other works, the processes adopted were emerging
processes with regime switching regulated by a Markov control process, which
is separate from our case where the process enters an absorbing state. Another
significant work in this direction is [1], in which the authors consider asset
allocation when the stock follows a geometric Brownian motion and a jump
process. Even though the stochastic processes considered are complicated, the
authors were able to obtain an exact solution to the Merton allocation problem
for logarithmic utility.

The optimization problem we consider in this work is inspired by similar
problems on the bond side. To the best of our knowledge, this type of problem
was first studied in [5]. The problem considered there was an optimization for
stock cash and a bond with some default probability. The authors formulated
the underlying optimization problem and presented a solution but omitted the
final derivation of the Hamilton-Jacobi-Bellman (HJB) equation or a stochastic
integral approach we employ in our paper. The problem was fully solved in [7,
8], where an optimal allocation problem for a bond with bankruptcy and a stock
was considered. The authors set up an HJB equation with two regimes: the post-
default and the pre-default regime. Our problem is separate from this scenario
due to the fact that we consider stocks with exogenous bankruptcy occurrence.
To the best of our knownledge, this type of problems with exogenous
bankruptcies for stocks have not been previously considered. Our solution
presented in Equation (10) is interesting on its own account as it mandates
that one should never allocates the complete amount of investment into a stock
with with positive bankruptcy probability or borrow money to invest in such
stocks. It is also worthwhile to note that bankruptcy is an absorbing process
and, as far as we know for optimal portfolios, such processes have not been
considered before. The probability of bankruptcy can be significant for stocks
with low credit quality and, thus, our solution is not a mere theoretical
exercise but has a practically important bearing as we will show in the
sequel. It should be reiterated that the process we consider is different from
conventional models based on Poisson jump processes since the stock cannot
recover from bankruptcy so that the latter is an absorbing state. On the other
hand, we note that it is possible to consider similar models (cf. [4, 32])
that would incorporate semi-Markovian switching (viz. [24, Chapter 7]) between
economic regimes that can be either absorbing or not absorbing. In this sense,
our present work provides a blueprint for investigating optimal portfolio
allocation problems for this type of models.

The rest of this paper is structured as follows. In Section 2, we introduce a
probabilistic framework and state the problem. In Section 3, a simple
heuristic calculation is presented to derive the Merton ratio from Equation
(10). Without being rigorous, it is new and interesting on its own merit since
it employs a non-conventional form of the Itô’s lemma to formally derive the
optimality conditions following the idea originally presented in [19]. Similar
to [4, 7], a pair of coupled HJB equations is derived to find the optimal
weights. In Section 4, a rigorous solution to the optimal portfolio problem is
obtained by reducing the expected terminal logarithmic utility objective for a
stock with bankruptcy to a combination of exponentially weighted expected
terminal and running logarithmic utilities without bankruptcy. The optimal
Merton ratio (viz. Equation (10)) obtained with both approaches agree. In
Section 5, we then illustrate our approach by compute the optimal allocation
for a real-world stock with a non-trivial bankruptcy probability and confirm a
material difference from the allocation produced with the usual Merton ratio.
Summary and conclusions are given in Section 6. Proofs of auxiliary results
are relegated to Appendix A.

##  2\. Probabilistic Setup and Preliminaries

For time $t\geq 0$, let $S_{t}$ and $R_{t}$ denote the price of a defaultable
stock following Equation (1) driven by a standard Brownian motion
$(B_{t})_{t\geq 0}$ and and a riskfree asset satisfying Equation (2), where
$\mu$, $\sigma>0$ and $r>0$ are the drift, volatility and risk-free rate,
respectively. Assuming $W_{t}$ is the total wealth at time $t$ consisting of
$\pi_{t}W_{t}$ dollars in stock and $(1-\pi_{t})X_{t}$ dollars in the riskless
asset for some $\pi_{t}$, the dynamics of $X_{t}$ is given by equation

(11)
$\begin{split}\mathrm{d}W_{t}&=\left\\{\begin{array}[]{cl}\big{(}\mu\pi_{t-}+r(1-\pi_{t-})\big{)}W_{t-}\mathrm{d}t+\sigma\pi_{t-}W_{t-}\mathrm{d}B_{t}&\text{for
}0\leq t<\tau,\\\ r(1-\pi_{t})W_{t}&\text{for }\tau\leq t\leq
T,\end{array}\right.\\\ \Delta W_{\tau}&=-\pi_{t}W_{\tau-},\\\
W_{0}&=w\end{split}$

where $w>0$ is a (given) initial wealth. Recall that
$\tau\sim\mathrm{Exp}(\tau)$ is assumed stochastically independent of
$(B_{t})_{t\geq 0}$. To emphasize the dependence of $W_{t}$ on $\pi_{t}$ and
$\tau$, we will often write $W_{t}^{\pi,\tau}$.

To put Equation (11) into a formal framework similar to [23], consider the
semi-Markov process

(12) $\xi_{t}:=\mathds{1}_{\\{t<\tau\\}}$

with the non-absorbing pre-default state $0$ and the absorbing post-default
state $1$. Introducing the drift and volatility functions

(13)
$\begin{split}a(t,W,\pi,\xi)&=\big{(}\mu\pi+r(1-\pi)\big{)}(1-\xi)+(1-\pi)\xi,\\\
b(t,W,\pi,\xi)&=\sigma\pi(1-\xi),\end{split}$

Equation (11) can be expressed as

(14)
$\begin{split}\mathrm{d}W_{t}&=a(t-,W_{t-},\pi_{t-},\xi_{t-})\mathrm{d}t+b(t-,W_{t-},\pi_{t-},\xi_{t-})\mathrm{d}B_{t}\text{
for }0\leq t\leq T,\\\ \Delta W_{\tau}&=-\pi_{\tau-}W_{\tau-}\text{ if
}\tau\leq T,\\\ W_{0}&=w.\end{split}$

It should be emphasized that, unlike [23], the switching time in Equation (14)
is random.

Consider the filtered probability space $(\Omega,(\mathcal{F}_{t})_{t\geq
0},\mathbb{P})$ with the minimal $\sigma$-algebra

(15) $\mathcal{F}_{t}=\sigma(B_{t},\xi_{t})\text{ for }t\geq 0.$

Arguing similar to [4], for any adapted process $(\pi_{t})_{t\geq 0}$ with

$\int_{0}^{T}\pi_{t}^{2}\mathrm{d}t<\infty\;\mathbb{P}\text{-a.s.},$

Equation (14) possesses a unique strong solution (being a cádlág process)
denoted by $(W_{t}^{\pi,\tau})_{t\geq 0}$ on any time interval $[0,T]$. The
set of admissible controls for Equation (14) can be defined as

(16) $\begin{split}\mathcal{A}_{w}:=\big{\\{}\pi\,|\,&(\pi)_{t\in[0,T]}\text{
is adapted with
}\int_{0}^{T}\pi_{t}^{2}\mathrm{d}s<\infty\;\mathbb{P}\text{-a.s.},\\\
&(W_{t})_{t\in[0,T]}\text{ is a unique strong solution such that
}W_{0}^{\pi,\tau}=w,\\\
&W_{t}^{\pi,\tau}(\omega)>0\;\mathbb{P}\otimes\nu\text{ a.e.~{}in
}[0,T]\times\Omega,\;\mathbb{E}\big{[}\mathop{\operatorname{ess\,sup}}_{t\in[0,T]}U(W_{t}^{\pi,\tau})_{-}\big{]}<\infty\big{\\}}\end{split}$

where $x_{-}=\max\\{0,-x\\}$ denotes the negative part of $x$ and $\nu$ is the
Borel measure on $[0,\infty)$.

Employing the logarithmic utility function

(17) $U(W)=\log(W),$

consider the terminal expected utility functional

(18) $J(\pi)=\mathbb{E}_{\mathcal{F}_{T}}\big{[}U(W^{\pi,\tau}_{T})\big{]}.$

Thus, the optimal portfolio problem reads as

(19)
$J(\pi)\to\min_{\pi\in\mathcal{A}_{w}},\quad\pi^{\ast}=\mathop{\operatorname{arg\,max}}_{\pi\in\mathcal{A}_{w}}J(W^{\pi,\tau},\pi),\quad
W_{t}^{\ast}=W_{t}^{\pi^{\ast},\tau}.$

To solve the optimal control problem, we first want to explicitly characterize
the forward map $\pi\mapsto W^{\pi,\tau}$ mapping each admissible control
$\pi\in\mathcal{A}_{w}$ to the unique strong solution
$(W^{\pi,\tau}_{t})_{t\geq 0}$ of Equation (11). To this end, let
$(X_{t})_{t\geq 0}$ denote the unique strong solution to the classical wealth
SDE without bankruptcy

(20)
$\mathrm{d}X_{t}=\big{(}\mu\pi_{t}+r(1-\pi_{t})\big{)}X_{t}\mathrm{d}t+\sigma\pi_{t}X_{t}\mathrm{d}B_{t}\text{
for }t>0,\quad X_{0}=w$

explicitly given [18, Chapter 4] via

(21)
$X_{t}=\Phi_{t}w-\sigma\int_{0}^{t}\Phi_{t}\Phi_{s}^{-1}\pi_{s}\mathrm{d}s$

with

$\Phi_{t}=\exp\Big{(}\int_{0}^{t}\Big{(}\big{(}\mu\pi_{s}+r(1-\pi_{s})\big{)}-\frac{1}{2}\sigma^{2}\pi_{s}^{2}\Big{)}\mathrm{d}s+\sigma\int_{0}^{t}\pi_{s}\mathrm{d}B_{s}\Big{)}.$

With this notation, the wealth process $(W_{t})_{t\in[0,T]}$ solving Equation
(11) can be expressed as

(22)
$W_{t}=X_{t}\mathds{1}_{\\{t<\tau\\}}+\exp\Big{(}\int_{\tau}^{t}r(1-\pi_{s})\mathrm{d}s\Big{)}(1-\pi_{\tau-})X_{\tau}\mathds{1}_{\\{t\geq\tau\\}}.$

Indeed, if the bankruptcy has not happened at time $t$ or before, $W_{t}$ is
same as $X_{t}$. Otherwise, all capital allocated in the stock at time $\tau-$
(viz. $\pi_{\tau-}X_{\tau-}$) will be lost, while the remaining capital
allocated in the riskless asset (viz. $(1-\pi_{\tau-})X_{\tau-}$) will
continue accruing from time $\tau$ to time $t$.

As can be seen from Equation (22), the positivity of $W_{t}>0$
$\mathbb{P}\otimes\nu$-a.e. required in the definition of $\mathcal{A}_{w}$ in
Equation (23) necessitates that

$\pi_{t}(\omega)<1\text{ for }\mathbb{P}\otimes\nu\text{ for a.e.
}(t,\omega)\in\Omega\times[0,T].$

Indeed, letting $I=\\{t\in[0,T]\,|\,\pi_{t}\geq 1\\}$, we obtain

$\begin{split}\mathbb{P}\\{W_{t}\leq
0\\}&=\mathbb{P}\Big{\\{}\exp\Big{(}\int_{\tau}^{t}r(1-\pi_{s})\mathrm{d}s\Big{)}(1-\pi_{\tau-})X_{\tau}\mathds{1}_{\\{T\geq\tau\\}}\leq
0\Big{\\}}\\\
&=\mathbb{P}\Big{\\{}(1-\pi_{\tau-})X_{\tau}\mathds{1}_{\\{\tau\in I\\}}\leq
0\Big{\\}}=\mathbb{P}\\{\tau\in I\\},\end{split}$

where we used the fact that $X_{t}>0$ $\mathbb{P}$-a.s. for a.e. $t\in[0,T]$.
The latter probability vanishes if and only if $\nu(I)=0$.

Thus, the admissible set in Equation (23) can be equivalently expressed as

(23) $\begin{split}\mathcal{A}_{w}:=\big{\\{}\pi\,|\,&(\pi)_{t\in[0,T]}\text{
is adapted with
}\int_{0}^{T}\pi_{t}^{2}\mathrm{d}t<\infty\;\mathbb{P}\text{-a.s.},\;\pi_{t}(\omega)<1\;\mathbb{P}\otimes\nu\text{
a.e.~{}in }\Omega\times[0,T],\\\ &(W_{t})_{t\in[0,T]}\text{ is a unique strong
solution with
}W_{0}^{\pi,\tau}=w,\;\mathbb{E}\big{[}\mathop{\operatorname{ess\,sup}}_{t\in[0,T]}U(W_{t}^{\pi,\tau})_{-}\big{]}<\infty\big{\\}}\end{split}$

Analyzing Equation (22), we easily observe that any wealth fraction $\pi_{t}$
invested into stock after bankruptcy (i.e., $t\geq\tau$) in lieu of the
riskfree asset can only reduce the terminal wealth, $W_{T}$, and, thus, the
expected utility $J(\cdot)$. Therefore, without loss of generality, we can
assume $\pi_{t}=0$ for $t\geq\tau$ leading to state-switching controls

(24)
$\pi_{t}=\pi^{\ast}_{\mathrm{pre}}\mathds{1}_{\\{t<\tau\\}}+\pi^{\ast}_{\mathrm{post}}\mathds{1}_{\\{t\geq\tau\\}}=\left\\{\begin{array}[]{cl}\pi_{\mathrm{pre},t},&t<\tau,\\\
\pi_{\mathrm{post},t},&t\geq\tau\end{array}\right.\quad\text{ with
}\pi_{\mathrm{post},t}\equiv 0.$

This furnishes a simplified version of Equation (22) given by

(25)
$W_{t}=X_{t}\mathds{1}_{\\{t<\tau\\}}+e^{r(t-\tau)}(1-\pi_{\tau-})X_{\tau}\mathds{1}_{\\{t\geq\tau\\}}$

without violating any admissibility conditions.

##  3\. Heuristic Derivation

We begin with a heuristic calculation of the optimal allocation strategy. To
this end, we define the admissible set

$\begin{split}\mathcal{A}_{t,w}:=\big{\\{}\pi\,|\,&(\pi)_{s\in[t,T]}\text{ is
adapted with
}\int_{t}^{T}\pi_{s}^{2}\mathrm{d}s<\infty\;\mathbb{P}\text{-a.s.},\;\pi_{s}(\omega)<1\text{
for }\mathbb{P}\otimes\nu\text{ a.e.~{}in }\Omega\times[t,T],\\\
&(W_{s})_{s\in[t,T]}\text{ is a unique strong solution with
}W_{t}^{\pi,\tau}=w,\;\mathbb{E}\big{[}\mathop{\operatorname{ess\,sup}}_{s\in[t,T]}U(W_{s}^{\pi,\tau})_{-}\big{]}<\infty\big{\\}}\end{split}$

where $(W_{s}^{\pi,\tau})_{s\in[t,T]}$ solves the jump SDE in Equation (11)
over $[t,T]$. Further, consider the value function

(26)
$V(t,w)=\sup_{\pi\in\mathcal{A}_{t,w}}\mathbb{E}_{\mathcal{F}_{T}}\big{[}U(W^{\pi,\tau}_{T})\,|\,W_{t}^{\pi,\tau}=w\big{]}.$

In spirit of [4], we also introduce two auxiliary “conditional” value
functions

(27)
$\begin{split}V^{\mathrm{pre}}(t,w)&=\sup_{\pi\in\mathcal{A}_{t,w}}\mathbb{E}_{\mathcal{F}_{T}}\big{[}U(W^{\pi,\tau}_{T})\,|\,t<\tau\big{]},\\\
V^{\mathrm{post}}(t,w)&=\sup_{\pi\in\mathcal{A}_{t,w}}\mathbb{E}_{\mathcal{F}_{T}}\big{[}U(W^{\pi,\tau}_{T})\,|\,t=\tau\big{]}=\mathbb{E}_{\mathcal{F}_{T}}\big{[}U(W^{\pi=0,\tau}_{T})\,|\,t=\tau\big{]}\end{split}$

corresponding to the case no bankruptcy has occurred at time $t$ or before and
the case the bankruptcy occurs at time $t$, respectively. It should be noted
that this framework allows for default-specific optimal control processes
$\pi_{t}$, i.e.,

Using the law of total probability, the total value function in Equation (26)
can be expressed as

(28)
$\begin{split}V(t,w)&=V^{\mathrm{pre}}(t,w)\,\mathbb{P}\\{\tau<t\\}+V^{\mathrm{post}}(t,w)\,\mathbb{P}\\{\tau\geq
t\\}\\\ &=(1-e^{-\lambda t})V^{\mathrm{pre}}(t,w)+e^{-\lambda
t}V^{\mathrm{post}}(t,w).\end{split}$

Arguing similar to Equation (25), we can explicitly compute the “post”
conditional value function via

(29)
$V^{\mathrm{post}}(t,w)=U\big{(}e^{r(T-t)}(1-\pi_{t-})w\big{)}=r(T-t)+\log(1-\pi_{t-})+\log(w),$

recalling that $U(\cdot)$ is the logarithmic utility function (viz. Equation
(17)).

As for the “pre” conditional value function, on the strength of Equation (9),
the conditional probability of no crash occurring in $[t,t+\mathrm{d}t]$
assuming no crash has occurred until time $t$ is
$(1-\lambda\mathrm{d}t)+o(\mathrm{d}t)$, while the conditional probability of
crash is $\lambda\mathrm{d}t+o(\mathrm{d}t)$. Thus, using an argumentation
similar to [4, Theorem 3.3], we can heuristically derive an infinitesimal
recursive formula for $V^{\mathrm{pre}}(t,w)$ reading as

(30)
$\begin{split}\mathrm{d}V^{\mathrm{pre}}(t,w)&=(1-\lambda\mathrm{d}t)\mathrm{d}V^{\mathrm{pre}}(t,w)+\lambda\mathrm{d}t\big{(}V^{\mathrm{post}}(t,w)-V^{\mathrm{pre}}(t,w)\big{)}+o(\mathrm{d}t)\\\
&=(1-\lambda\mathrm{d}t)\mathrm{d}V^{\mathrm{pre}}(t,w)+\lambda\mathrm{d}t\big{(}r(T-t)+\log(1-\pi_{t-})+\log(w)-V^{\mathrm{pre}}(t,w)\big{)}+o(\mathrm{d}t)\end{split}$

Let $\pi_{t}\in\mathcal{A}_{t,w}$. For a wealth process
$(W_{s}^{\pi,\tau})_{s\in[t,T]}$ solving the jump SDE in Equation (11) subject
to the initial condition $W_{t}^{\pi,\tau}=w$, we evaluate the differential
$\mathrm{d}V^{\mathrm{pre}}(s,W_{s})$ for $t\leq s\leq T$. According to
Equation (20), the wealth process satisfies

(31)
$\begin{split}\mathrm{d}W_{s}&=W_{s}\left(\pi_{s}\mu\mathrm{d}s+\pi_{s}\sigma\mathrm{d}B_{s}+(1-\pi_{s})r\mathrm{d}s\right),\\\
(\mathrm{d}W_{s})^{2}&=W_{s}^{2}\pi_{s}^{2}\sigma^{2}\mathrm{d}s.\end{split}$

Assuming $V^{\mathrm{pre}}(\cdot)$ is a smooth function, Itô’s rule (cf. [4,
Theorem 3.5]) furnishes

(32)
$\begin{split}\mathrm{d}V^{\mathrm{pre}}(t,w)&=\partial_{w}V^{\mathrm{pre}}(t,w)(w\left(\pi_{t}\mu\mathrm{d}t+(1-\pi_{t})r\mathrm{d}t\right)+\frac{1}{2}\partial_{w}^{2}V^{\mathrm{pre}}(t,w)w^{2}\pi_{t}^{2}\sigma^{2}\mathrm{d}t+\partial_{t}V^{\mathrm{pre}}(t,w)\mathrm{d}t\\\
&+\lambda\mathrm{d}t\left(r(T-t)+\log(1-\pi_{t})+\log(w)-V^{\mathrm{pre}}(t,w)\right),\end{split}$

where we assume the continuity of $\pi_{s}$ (and, therefore, that of
$\pi_{\mathrm{pre},s}$) for $t\leq s<\tau$ implying $\pi_{s-}=\pi_{s}$.
Dividing by $\mathrm{d}t$ and passing to the limit $\mathrm{d}t\to 0$, we
obtain an HJM-type equation for the “pre” conditional value function

(33)
$\partial_{t}V^{\mathrm{pre}}+\mathcal{H}(t,V^{\mathrm{pre}},\pi_{t},\partial_{w}V^{\mathrm{pre}},\partial_{w}^{2}V^{\mathrm{pre}})=0$

where

(34)
$\begin{split}\mathcal{H}(t,w,\pi,V^{\mathrm{pre}},\partial_{w}V^{\mathrm{pre}},\partial_{ww}V^{\mathrm{pre}})=\sup_{\pi_{t}<1}&\Big{(}\partial_{w}V^{\mathrm{pre}}(w\left(\pi_{t}\mu+(1-\pi_{t})r\right)+\tfrac{1}{2}\partial_{w}^{2}V^{\mathrm{pre}}w^{2}\pi_{t}^{2}\sigma^{2}\\\
&+\lambda\left(r(T-t)+\log(1-\pi_{t})+\log(w)-V^{\mathrm{pre}}\right)\Big{)}.\end{split}$

In contrast to the usual HJB equation, we can observe that our HJM operator
$\mathcal{H}(\cdot)$ additionally depends on the derivatives of the value
function but also the value function itself.

Using the ansatz

(35) $V(t,w)=\log(w)+f(t),$

the first-order Fermat optimality condition for $g(\cdot)$ reads as

(36) $(\mu-r)-\pi_{t}\sigma^{2}-\frac{\lambda}{1-\pi_{t}}=0.$

Recalling $\pi_{t}^{\ast}=0$ for $t\geq\tau$ (viz. Equation (24)), the optimal
allocation (cf. Equation (10)) is given as:

$\pi_{t}^{\ast}=\pi^{\ast}_{\mathrm{pre}}\mathds{1}_{\\{t<\tau\\}}$

with

$\pi_{\mathrm{pre}}^{\ast}=\frac{\mu-r+\sigma^{2}-\sqrt{\left(\sigma^{2}-(\mu-r)\right)^{2}+4\lambda\sigma^{2}}}{2\sigma^{2}}.$

Plugging $\pi_{t}^{\ast}$ back into Equation (34), the latter can be
explicitly solved using the ansatz from Equation (35) yielding the value
function $V=V(t,x)$. One would naturally expect the latter to satisfy

$\max_{\pi\in\mathcal{A}_{t,w}}J(\pi)=V(t,x)=J(\pi^{\ast}).$

Unfortunately, the optimal control problem in Equation (34) is not directly
covered by the verification theorem for Markov-switching SDEs established in
[4]. While latter results can likely be extended to our situation, an simpler
strategy to avoid laborious proofs and calculations is to transform the
problem in Equation (34) to a conventional controlled Markovian diffusion and
invoking one of several rigorous verification theorems known in the literature
for this class of optimal control problems.

##  4\. Reduction to Controlled Markovian Diffusion

In this section, we present an alternative derivation that reduces the non-
conventional optimal control problem discussed in Section 2 to the classical
case. For the sake of simplicity, with some abuse of notation, we will write
$\pi_{t}$ to denote $\pi_{\mathrm{pre},t}$ for the most part of this section
until Equation (45), where $\pi_{t}$ reclaims its original meaning (cf.
Equation (24)).

Our approach is inspired by [17]. Assuming the conditions of Fubini &
Tonelli’s theorem are satisfied, we can write

$\displaystyle\mathbb{E}_{\mathbb{P}_{W}\otimes\mathbb{P}_{\tau}}\big{[}u(Y_{T}^{\pi})\big{]}$
$\displaystyle=\mathbb{E}_{\mathbb{P}_{W}}\Big{[}\mathbb{E}_{\tau}\big{[}u(Y_{T}^{\pi})\big{]}\Big{]}=\mathbb{E}_{\mathbb{P}_{W}}\Big{[}\int_{0}^{\infty}u(Y_{T}^{\pi})\varphi(t)\mathrm{d}t\Big{]}$
$\displaystyle=\mathbb{E}_{\mathbb{P}_{W}}\Big{[}\int_{0}^{\infty}\Big{(}u(X_{T})\mathds{1}_{\\{T<t\\}}+u\big{(}e^{r(T-t)}(1-\pi_{t})X_{t}\big{)}\mathds{1}_{\\{T\geq
t\\}}\Big{)}\varphi(t)\mathrm{d}t\Big{]}$
$\displaystyle=\mathbb{E}_{\mathbb{P}_{W}}\Big{[}\int_{T}^{\infty}u(X_{T})\varphi(t)\mathrm{d}t\Big{]}+\mathbb{E}_{\mathbb{P}_{W}}\Big{[}\int_{0}^{T}u\big{(}e^{r(T-t)}(1-\pi_{t})X_{t}\big{)}\varphi(t)\mathrm{d}t\Big{]}$
$\displaystyle=\mathbb{E}_{\mathbb{P}_{W}}\Big{[}\Big{(}\int_{T}^{\infty}\varphi(t)\mathrm{d}t\Big{)}u(X_{T})\Big{]}+\mathbb{E}_{\mathbb{P}_{W}}\Big{[}\int_{0}^{T}u\big{(}e^{r(T-t)}(1-\pi_{t})X_{t}\big{)}\varphi(t)\mathrm{d}t\Big{]}$
$\displaystyle=\mathbb{P}_{\mathbb{P}_{\tau}}\\{\tau\geq
T\\}\mathbb{E}_{\mathbb{P}_{W}}\big{[}u(X_{T})\big{]}+\mathbb{E}_{\mathbb{P}_{W}}\Big{[}\int_{0}^{T}u\big{(}e^{r(T-t)}(1-\pi_{t})X_{t}\big{)}\varphi(t)\mathrm{d}t\Big{]}$

where $\mathbb{P}_{\mathbb{P}_{\tau}}\\{\tau\geq
T\\}=\int_{T}^{\infty}\varphi(t)\mathrm{d}t$ solely depends on $T$ with
$\varphi(\cdot)$ denoting the probability density function of $\tau$. As
before, we choose $u(x)=\log(x)$ to be the logarithmic utility and assume
$\tau\sim\mathrm{Exp}(\lambda)$ is exponentially distributed with some known
intensity $\lambda>0$, i.e., $\varphi(t)=\lambda\exp(-\lambda
t)\mathds{1}_{\\{t\geq 0\\}}$. Then the expected utility can be further
simplified as

(37)
$\displaystyle\begin{split}\mathbb{E}_{\mathbb{P}_{W}\otimes\mathbb{P}_{\tau}}\big{[}u(Y_{T}^{\pi})\big{]}&=\mathbb{E}_{\mathbb{P}_{W}}\Big{[}\int_{0}^{T}\lambda
e^{-\lambda t}\log\big{(}(1-\pi_{t})X_{t}\big{)}\mathrm{d}t\Big{]}\\\
&+\mathbb{E}_{\mathbb{P}_{W}}\big{[}e^{-\lambda
T}\log(X_{T})\big{]}+r\lambda\int_{0}^{T}(T-t)e^{-\lambda t}\mathrm{d}t\\\
&=\mathbb{E}_{\mathbb{P}_{W}}\Big{[}\int_{0}^{T}\lambda e^{-\lambda
t}\log\big{(}(1-\pi_{t})X_{t}\big{)}\mathrm{d}t\Big{]}\\\
&+\mathbb{E}_{\mathbb{P}_{W}}\Big{[}e^{-\lambda
T}\log(X_{T})+r\Big{(}T-\frac{1-e^{-\lambda
T}}{\lambda}\Big{)}\Big{]}\end{split}$

where $\tau$ is completely eliminated from the latter equation.

Therefore, we were able to reduce the original control problem with sporadic
bankruptcy in Equations (11), (17), (18) and (19) to a classical stochastic
control problem in Equations (20) and (37) with running and terminal costs for
wealth $W_{t}$ so that the usual Bellman’s principle for Markovian diffusions
can now be applied. Throughout the rest of the section, we adopt the notation
of [13, Chapter IV]. Note that maximization of expected utility is
equivalently expressed as minimimization of negative utility within this
framework. Thus, we seek to minimize the functional

(38) $J(\pi)=\mathbb{E}\Big{[}-\int_{0}^{T}\lambda e^{-\lambda
t}\log\big{(}(1-\pi_{t})X_{t}\big{)}\mathrm{d}t\Big{]}+\mathbb{E}\Big{[}-e^{-\lambda
T}\log(X_{T})-r\Big{(}T-\frac{1-e^{-\lambda T}}{\lambda}\Big{)}\Big{]}.$

Since the terminal cost function in Equation (38) cannot be continuously
extended to $x=0$, classical results (viz. [12, 13]) cannot be directly
applied to our problem. Instead, we will employ the verification theorem
proved in the recent work [2]. Other known results, e.g., [14, 15, 29], do not
appear to be directly applicable in our situation.

The set of admissible controls is defined as

(39)
$\begin{split}\mathcal{A}_{t,x}:=\big{\\{}\pi\,|\,&(\pi)_{s\in[t,T]}\text{ is
adapted with
}\int_{t}^{T}\pi_{s}^{2}\mathrm{d}s<\infty\;\mathbb{P}\text{-a.s.},\\\
&(X_{s})_{s\in[t,T]}\text{ is a unique strong solution to
SDE~{}\eqref{EQUATION:SDE_WEALTH} with }X_{t}^{\pi}=x,\\\
&\mathbb{E}\big{[}\mathop{\operatorname{ess\,sup}}_{s\in[t,T]}\big{(}\log(X_{s}^{\pi})\big{)}_{-}\big{]}<\infty\big{\\}}\end{split}$

where, unlike Equation (15), the filtration is now solely generated by the
Brownian motion $(B_{t})_{t}$. Arguing similar to previous sections, the last
condition in the definition of $\mathcal{A}_{t,x}$ in Equation (39)
necessitates

$\pi_{t}(\omega)<1\;\mathbb{P}\otimes\nu\text{ a.e.~{}in }\Omega\times[0,T].$

Introducing the Hamilton-Jacobi-Bellman (HJB) operator

$\displaystyle(t,V,\pi)\mapsto\mathcal{H}(t,\pi,\partial_{x}V,\partial_{xx}V)$

with

(40)
$\displaystyle\begin{split}\mathcal{H}(t,x,p,A)&=\inf_{\pi\in\mathbb{R}}\Big{(}f(t,x,\pi)\cdot
p+\frac{1}{2}\sigma^{2}(t,x,\pi)\cdot A-L(t,x,\pi)\Big{)}\\\
&=\inf_{\pi<1}\Big{(}f(t,x,\pi)\cdot p+\frac{1}{2}\sigma^{2}(t,x,\pi)\cdot
A-L(t,x,\pi)\Big{)}\end{split}$

where $f(t,x,\pi)=\big{(}\mu\pi+r(1-\pi)\big{)}x$, $\sigma(t,x,\pi)=\sigma\pi
x$, $L(t,x,\pi)=\lambda e^{-\lambda t}\log\big{(}(1-\pi)x\big{)}$, consider
the terminal value problem for the HJB partial differential equation (PDE)

(41)
$\displaystyle\partial_{t}V(t,x)+\mathcal{H}\big{(}t,x,\partial_{x}V(t,x),\partial_{xx}V(t,x)\big{)}$
$\displaystyle=0$ $\displaystyle\text{for }(t,x)\in(0,T)\times(0,\infty),$
(42) $\displaystyle V(T,x)$ $\displaystyle=\Psi(T,x)$ $\displaystyle\text{for
}x\in(0,\infty)$

with the terminal value $\Psi(T,x)=-e^{-\lambda
T}\log(x)-r\Big{(}T-\frac{1-e^{-\lambda T}}{\lambda}\Big{)}$. On the strength
of Lemma 3 in the Appendix, the HJB operator can be expressed as a fully
nonlinear elliptic operator

(43)
$\mathcal{H}(t,x,p,A)=\big{(}\mu\pi^{\ast}+r(1-\pi)\big{)}xp+\frac{1}{2}\sigma^{2}(\pi^{\ast})^{2}x^{2}-\lambda
e^{-\lambda t}\log\big{(}(1-\pi^{\ast})x\big{)}\text{ for }A>0$

with $\pi^{\ast}=\pi^{\ast}(t,x,p,A)$ given in Equation (50). For $A\leq 0$,
the HJB operator is undefined since the underlying minimization problem is
ill-posed.

To solve the HJB Equations (41)–(42), we use the the ansatz

(44) $V(t,x)=-e^{-\lambda t}\log(x)+G(t)$

with $\partial_{xx}V(t,x)=e^{-\lambda t}(1/x^{2})>0$ for $t\geq 0$ and $x>0$.
Plugging into Equation (50) yields the (constant) pre-default Merton ratio:

(45)
$\pi^{\ast}_{\mathrm{pre}}\equiv\frac{(\mu-r)+\sigma^{2}-\sqrt{(\sigma^{2}-(\mu-r))^{2}+4\lambda\sigma^{2}}}{2\sigma^{2}}.$

Like the classical Merton ratio, the Merton ratio in Equation (45) neither
depends on $x$, nor $t$. See Figure 1.

Figure 1.  Pre-default Merton ratio $\pi^{\ast}_{\mathrm{pre}}$ vs. $\lambda$

Moreover, for $\mu-r\leq\sigma^{2}$, $\pi^{\ast}_{\mathrm{pre}}(\lambda)$
converges to the classical expression $\frac{\mu-r}{\sigma^{2}}$ as
$\lambda\searrow 0$. It is also interesting to observe that our Merton ratio
in Equation (45) always suggests to put less money into the stock than the
usual Merton ratio would recommend. In particular, it also encourages to short
the stock more aggressively, especially for larger $\lambda$, which makes
practical sense in view of the chance of default. For $\mu-r<\sigma^{2}$,
$\pi^{\ast}_{\mathrm{pre}}(\lambda)\nearrow 1$ as $\lambda\searrow 0$. Unlike
classical situation, this means the investor should never invest all of his
capital into stock or borrow from the bank to purchase stock. Indeed, since
the risk of the stock price to default sporadically over any non-trivial time
interval is positive, a Merton ratio $\pi\geq 1$ is not admissible in the
sense of Equation (39) as it renders the terminal wealth $Y_{T}$ in Equation
(22) zero or negative with a positive probability, therefore, rendering the
expected utility equal $-\infty$. This is a major qualitative difference to
the classical scenario.

######  Proposition 1.

The value function is given by

(46) $V(t,x)=-e^{-\lambda
t}\log(x)-\frac{C^{\ast}(\lambda)}{\lambda}\big{(}e^{-\lambda t}-e^{-\lambda
T}\big{)}-r\Big{(}T-\frac{1-e^{-\lambda T}}{\lambda}\Big{)}\quad\text{ for
}t\geq 0,x>0$

with $\pi^{\ast}_{\mathrm{pre}}$ from Equation (45) and

$C^{\ast}(\lambda)=\big{(}\mu\pi^{\ast}_{\mathrm{pre}}+r(1-\pi^{\ast}_{\mathrm{pre}})\big{)}-\frac{1}{2}\sigma^{2}(\pi^{\ast}_{\mathrm{pre}})^{2}+\lambda\log(1-\pi^{\ast}_{\mathrm{pre}}).$

###### Proof.

Plugging

$\partial_{t}V(t,x)=\lambda e^{-\lambda
t}\log(x)+\dot{G}(t),\quad\partial_{x}V(t,x)=e^{-\lambda
t}(1/x),\quad\partial_{xx}V(t,x)=-e^{-\lambda t}(1/x^{2})$

into HJB Equation (41), we obtain

$\lambda e^{-\lambda
t}\log(x)+\dot{G}(t)-\big{(}\mu\pi^{\ast}_{\mathrm{pre}}+r(1-\pi^{\ast}_{\mathrm{pre}})\big{)}e^{-\lambda
t}+\frac{1}{2}\sigma^{2}(\pi^{\ast}_{\mathrm{pre}})^{2}e^{-\lambda t}-\lambda
e^{-\lambda t}\log(1-\pi^{\ast}_{\mathrm{pre}})-\lambda e^{-\lambda
t}\log(x)=0$

or, equivalently,

(47) $\dot{G}(t)=Ce^{-\lambda t}\text{ for }t\geq 0$

with

(48)
$C=\big{(}\mu\pi^{\ast}_{\mathrm{pre}}+r(1-\pi^{\ast}_{\mathrm{pre}})\big{)}-\frac{1}{2}\sigma^{2}(\pi^{\ast}_{\mathrm{pre}})^{2}+\lambda\log(1-\pi^{\ast}_{\mathrm{pre}}).$

Integrating Equation (47), we obtain $G(t)=-\frac{C}{\lambda}e^{-\lambda t}+c$
for some constant $c$. On the strength of Equation (42), we have the terminal
condition

$-\frac{C}{\lambda}e^{-\lambda T}+c=-r\Big{(}T-\frac{1-e^{-\lambda
T}}{\lambda}\Big{)},$

whence $c=\frac{C}{\lambda}e^{-\lambda T}-r\Big{(}T-\frac{1-e^{-\lambda
T}}{\lambda}\Big{)}$ and

$G(t)=\frac{C}{\lambda}\big{(}e^{-\lambda T}-e^{-\lambda
t}\big{)}-r\Big{(}T-\frac{1-e^{-\lambda T}}{\lambda}\Big{)}$

for $t\geq 0$ and $x>0$. ∎

Interestingly, the following observations can be made. For $\mu-r<\sigma^{2}$,
using L’Hôpital’s rule, $V(t,x;\lambda)$ can be shown to converge to the
classical Merton value function, i.e.,

$\lim_{\lambda\searrow
0}V(t,x;\lambda)=-\log(x)-\Big{(}r+\frac{(\mu-r)^{2}}{\sigma^{2}}\Big{)}(T-t).$

In contrast, for $\mu-r>\sigma^{2}$,

$\lim_{\lambda\searrow
0}V(t,x;\lambda)=-\log(x)-\Big{(}\mu-\frac{1}{2}\sigma^{2}\Big{)}(T-t).$

since $\lim_{\lambda\searrow 0}\pi^{\ast}_{\mathrm{pre}}(\lambda)=1$ and
$\lim_{\lambda\searrow 0}C^{\ast}(\lambda)=\mu-\frac{1}{2}\sigma^{2}$.

Observing

$V(t,x),\partial_{t}V(t,x),\partial_{xx}V(t,x)\in
C^{0}([0,T]\times(0,\infty),\mathbb{R})\text{ with
}\partial_{xx}V(t,x)=e^{-\lambda t}x^{-2}>0\text{ and
}\pi^{\ast}_{t}\in\mathcal{A}_{tx},$

the conditions of Verification Theorem [2, Theorem 3] are met, furnishing the
following optimality result.

######  Proposition 2.

Let the wealth process be given by Equation (22).

  * •

The Merton ratio

$\pi^{\ast}_{t}=\pi^{\ast}_{\mathrm{pre}}\mathds{1}_{\\{t<\tau\\}}$

with $\pi^{\ast}_{\mathrm{pre}}$ given in Equation (45) minimizes the negative
expected utility $J(\pi)$ in Equation (38) over all admissible controls
$\pi\in\mathcal{A}_{tx}$ and satisfies

$J(\pi^{\ast})=\inf_{\pi\in\mathcal{A}_{tx}}J(\pi)=V(t,x).$

  * •

The optimal wealth process, $W^{\ast}_{t}$, is given by Equations (20)–(22)
with $\pi^{\ast}_{t}$ as above.

##  5\. Example: Bombardier Stock Data Analysis

To illustrate our proposed methodology, we consider the B-rated stock of
Bombardier Inc. (BBD-B.TO), a Canadian aerospace and transportation company,
quoted by Toronto Stock Exchange. The annualized rating default probability is
estimated as $\hat{\lambda}=2.99\%$ [20, Table 23]. For simplicity, we use
$\hat{\lambda}=3\%$ in our calculations. Downloading a three years’ worth of
daily closing prices from 2912021 to 2912024 from Yahoo! Finance [11], we
compute the log-differenced returns

$X_{i}=\log(S_{t_{i+1}})-\log(S_{t_{i}})\text{ for }i=1,2,\dots,n$

with $n=751$ furnishing the the annualized drift and volatility estimates (cf.
[3, 28])

(49)
$\hat{\mu}=\frac{n_{\mathrm{trade}}}{n}\sum_{i=1}^{n}X_{i},\quad\hat{\sigma}=\sqrt{n_{\mathrm{trade}}}\Big{(}\frac{1}{n-1}\sum_{i=1}X_{i}^{2}-\frac{n}{n-1}\Big{(}\frac{1}{n}\sum_{i=1}^{n}X_{i}\Big{)}^{2}\Big{)}^{1/2}$

where $n_{\mathrm{trade}}=252$ stands for the number of trading days. See
columns 1 and 2 in Table 1.

$\hat{\mu}$ | $\hat{\sigma}$ | $\hat{\lambda}$ | $\hat{\pi}_{\mathrm{classic}}^{\ast}$ | $\hat{\pi}_{\mathrm{pre}}^{\ast}$ | $\hat{\pi}_{\mathrm{post}}^{\ast}$
---|---|---|---|---|---
$40\%$ | $60\%$ | $3\%$ | $101.12\%$ | $74.21\%$ | $0\%$
Table 1.  Estimated model parameters $\hat{\mu}$, $\hat{\sigma}$,
$\hat{\lambda}$ (annualized) and resulting optimal allocations
$\hat{\pi}_{\mathrm{classic}}^{\ast}$ (both pre- and post-default) according
to Merton and $\hat{\pi}_{\mathrm{pre}}^{\ast}$ and
$\hat{\pi}_{\mathrm{post}}^{\ast}$ derived in this paper.

Using Equations , the classical Merton’s ratio
$\hat{\pi}_{\mathrm{classic}}^{\ast}$ recommends to invest $101.12\%$ of
capital in stock while only $74.32\%$ of capital is allocated in stock pre-
default under our new ratio $\hat{\pi}_{\mathrm{pre}}^{\ast}$. In addition to
offering a much less conservative investment strategy, classical Merton ratio
is not even admissible (cf. Equation (39)) in the presence of a sporadic
bankruptcy as it renders the expected terminal utility equal negative
infinity. Therefore, the difference between our framework and the classical
framework of Merton is not merely quantitative, but qualitative.

##  6\. Summary and Conclusions

In this work, we presented an explicit solution to the optimal wealth
allocation problem between a stock and a cash bond assuming the stock
bankruptcy is an exogenous event with a constant intensity $\lambda$. To the
best of our knowledge, neither this type of problems, nor the solution we
developed for the $\log$-utility have previously appeared in the literature.
The stochastic integral derivation for this type of processes appears to be
new as well. Our solution has an interesting feature that one should never
allocates all money in a stock that can go bankrupt to hedge against an entire
loss of one’s money. This seems to correlate with the investor behavior in
real market. The framework we developed is not limited to bankruptcy events
only. In fact, stochastic control problems with regime or security changes can
be solve using a similar method. These examples include convertible bonds
(i.e., bonds with a conversion feature), bonds with embedded options assuming
the probability of embedded options are governed by an exogenous probability,
mortgage bonds, etc. Our framework can also be applied in the context of
switching between economic regimes caused by exogenous events in the market.
We intend to investigate these and other scenarios in subsequent work.

## Supplementary Materials

The following materials are available at
https://github.com/mpokojovy/MertonAbsorbing:

  1. (1)

Python code to produce Figure 1.

  2. (2)

R code implementing the results of Section 5.

##  Appendix A Auxiliary Results

######  Lemma 3.

Let $g(t,x,\pi,p,A)=f(t,x,\pi)\cdot p+\frac{1}{2}\sigma^{2}(t,x,\pi)\cdot
A-L(t,x,\pi)$ with

$\displaystyle
f(t,x,\pi)=\big{(}\mu\pi+r(1-\pi)\big{)}x,\quad\sigma(t,x,\pi)=\sigma\pi
x,\quad L(t,x,\pi)=\lambda e^{-\lambda t}\log\big{(}(1-\pi)x\big{)}.$

Assuming $x>0$ and $A>0$,

$\displaystyle\mathop{\operatorname{arg\,min}}_{\pi<1}g(t,x,\pi,p,A)$
$\displaystyle=\pi^{\ast}(t,x,p,A),$ $\displaystyle\inf_{\pi<1}g(t,x,\pi,p,A)$
$\displaystyle=g\big{(}t,x,\pi^{\ast}(t,x,p,A),p,A\big{)}$

with

(50)
$\pi^{\ast}(t,x,p,A)=\frac{-(\mu-r)xp+\sigma^{2}x^{2}A-\sqrt{((\mu-r)xp+\sigma^{2}x^{2}A)^{2}+4\lambda
e^{-\lambda t}(\sigma^{2}x^{2}A)}}{2(\sigma^{2}x^{2}A)}.$

###### Proof.

The function

$\displaystyle g(t,x,\pi,p,A)=\big{(}\mu\pi+r(1-\pi)\big{)}x\cdot
p+\frac{1}{2}\sigma^{2}\pi^{2}x^{2}\cdot A-\lambda e^{-\lambda
t}\log\big{(}(1-\pi)x\big{)}$

is smooth with respect to all of its arguments $t\geq 0$, $x>0$, $\pi<1$,
$p\in\mathbb{R}$ and $A>0$. Since $\lim_{\pi\nearrow 0}g(t,x,\pi,p,A)=\infty$,
no minimum can be attained in a vicinity of $\pi=1$.

Computing the derivative

$\displaystyle\partial_{\pi}g(t,x,\pi,p,A)=(\mu-r)xp+\sigma^{2}\pi
x^{2}A+\frac{\lambda e^{-\lambda t}}{1-\pi},$

the Fermat first-order optimality equation $\partial_{\pi}g(t,x,\pi,p,A)=0$
reads as

$\sigma^{2}\pi x^{2}A+(\mu-r)xp+\frac{\lambda e^{-\lambda t}}{(1-\pi)}=0$

or, multiplying with $(\pi-1)$,

(51) $P(\pi;\lambda
t)\equiv\big{(}\sigma^{2}x^{2}A\big{)}\pi^{2}+\big{(}(\mu-r)xp-\sigma^{2}x^{2}A\big{)}\pi-\big{(}(\mu-r)xp+\lambda
e^{-\lambda t}\big{)}=0.$

Solving the former quadratic equation, we obtain two roots given by

$\displaystyle\pi_{1,2}$
$\displaystyle=\frac{-\big{(}(\mu-r)xp-\sigma^{2}x^{2}A\big{)}\pm\sqrt{\big{(}(\mu-r)xp-\sigma^{2}x^{2}A\big{)}^{2}+4\big{(}\sigma^{2}x^{2}A\big{)}\big{(}(\mu-r)xp+\lambda
e^{-\lambda t}\big{)}}}{2\big{(}\sigma^{2}x^{2}A\big{)}}$
$\displaystyle=\frac{-\big{(}(\mu-r)xp-\sigma^{2}x^{2}A\big{)}\pm\sqrt{\big{(}(\mu-r)xp+\sigma^{2}x^{2}A\big{)}^{2}+\lambda
e^{-\lambda
t}\big{(}\sigma^{2}x^{2}A\big{)}}}{2\big{(}\sigma^{2}x^{2}A\big{)}}.$

In view of $A>0$, both roots are real and satisfy $\pi_{1}\leq\pi_{2}$.
Moreover, since $P(\pi;\lambda,t)>P(\pi;0,t)$ for $\lambda>0$ and $t\geq 0$,
we observe

$\pi_{1}(\lambda t)<\pi_{1}(0)\quad\text{ and }\quad\pi_{2}(\lambda
t)>\pi_{2}(0)\quad\text{ for }\lambda>0\text{ and }t\geq 0.$

Consider the following two cases.

  * •

Assuming $(\mu-r)xp\leq\sigma^{2}Ax^{2}$, we obtain

$\pi_{1}(\lambda
t)<\pi_{1}(0)\equiv-\frac{\mu-r}{\sigma^{2}}\Big{(}\frac{p}{Ax}\Big{)}\leq
1\quad\text{ and }\quad\pi_{2}(\lambda t)>\pi_{2}(0)\equiv 1\text{ for
}\lambda>0\text{ and }t\geq 0$

and, thus, $\partial_{\pi}g(t,x,\pi,p,A)\leq 0$ for $\pi<\pi_{1}(\lambda t)$
and $\partial_{\pi}g(t,x,\pi,p,A)\geq 0$ for $\pi_{1}(\lambda t)<\pi<1$.

  * •

Otherwise, $(\mu-r)xp>\sigma^{2}Ax^{2}$, which implies

$\pi_{1}(\lambda t)<\pi_{1}(0)\equiv 1\quad\text{ and }\quad\pi_{2}(\lambda
t)\geq\pi_{2}(0)\equiv-\frac{\mu-r}{\sigma^{2}}\Big{(}\frac{p}{Ax}\Big{)}>1\text{
for }\lambda>0\text{ and }t\geq 0$

and, thus, again $\partial_{\pi}g(t,x,\pi,p,A)\leq 0$ for $\pi<\pi_{1}(\lambda
t)$ and $\partial_{\pi}g(t,x,\pi,p,A)\geq 0$ for $\pi_{1}(\lambda t)<\pi<1$.

Therefore, in either scenario, the global minimum of $g=g(\pi)$ is attained at
$\pi=\pi_{1}(\lambda t)$ for $\lambda>0$, $t\geq 0$ and $x>0$ assuming $A>0$.
∎

## References

  * [1] Y. Aït-Sahalia and T. R. Hurd.  Portfolio choice in markets with contagion.  Journal of Financial Econometrics, 14, 2016.
  * [2] S. Albosaily and S. Pergamenchtchikov.  Optimal investment and consumption for multidimensional spread financial markets with logarithmic utility.  Stats, 4:1012–1026, 2021.
  * [3] A. T. Anum and M. Pokojovy.  A hybrid method for density power divergence minimization with application to robust univariate location and scale estimation.  Communications in Statistics – Theory and Methods, 2023.
  * [4] N. Azevedo, D. Pinheiro, and G.-W. Weber.  Dynamic programming for a Markov-switching jump-diffusion.  Journal of Computational and Applied Mathematics, 267:1–19, 2014\.
  * [5] T. R. Bielecki and I. Jang.  Portfolio optimization with a defaultable security.  Asia Pacific Financial Markets, 13, 2006.
  * [6] M. Brennan, E. Schwartz, and R. Lagnado.  Strategic asset allocation.  Journal of Economic Dynamics and Control, Volume 21, 1997.
  * [7] A. Capponi and J. Figueroa Lopez.  Dynamic portfolio optimization with a defaultable security and regime switching.  Mathematical Finance, 24, 2011.
  * [8] A. Capponi and M. Larsson.  Default and systemic risk in equilibrium.  Mathematical Finance, 25, 2015.
  * [9] S. Emmer and C. Kluppelberg.  Optimal portfolios when stock prices follow an exponential Levý process.  Finance and Stochastics, 8, 2004.
  * [10] W. Fei.  Optimal control of uncertain stochastic systems with Markovian switching and its application to portfolio decisions.  Cybernetics and Systems – An International Journal, 45, 2014.
  * [11] Yahoo! Finance.  Bombardier Inc. (BBD-B.TO) Historical Prices.  https://finance.yahoo.com/quote/BBD-B.TO/history, February 2024\.
  * [12] W. H. Fleming and R. W. Rishel.  Deterministic and Stochastic Optimal Control, volume 25 of Stochastic Modelling and Applied Probability.  Springer, New York, NY, 1975.
  * [13] W. H. Fleming and H. M. Soner.  Controlled Markov Processes and Viscosity Solutions, volume 25.  Springer Science & Business Media, New York, NY, 2006.
  * [14] J.-P. Fouque, R. Sircar, and T. Zariphopoulou.  Portfolio optimization and stochastic volatility asymptotics.  Mathematical Finance, 27(3):704–745, 2017.
  * [15] T. Goll and J. Kallsen.  Optimal portfolios for logarithmic utility.  Stochastic Processes and their Applications, 89:31–48, 2000.
  * [16] I. Karatzas, J. Lehoczky, S. Sethi, and S. Shreve.  Explicit solution of a general consumption/investment problem.  Mathematics of Operations Research, 11(2):261–294, 1986.
  * [17] A. Klinger.  Control with stochastic stopping time.  International Journal of Control, 11:541–549, 1970.
  * [18] P. E. Kloeden and E. Platen.  Numerical Solution of Stochastic Differential Equations, volume 23 of Stochastic Modelling and Applied Probability.  Springer, Berlin, Heidelberg, 1992.
  * [19] Y. Kopeliovich.  Optimal control problems for stochastic processes with absorbing regime.  Journal of Stochastic Analysis, 4, 2023.
  * [20] N. W. Kraemer, J. Palmer, N. M. Richhariya, S. Iyer, L. Fernandes, A. Meher, and S. Pranshu.  Default, Transition, and Recovery: 2021 Annual Global Corporate Default and Rating Transition Study.  https://www.spglobal.com/ratings/en/research/articles/220413-default-transition-and-recovery-2021-annual-global-corporate-default-and-rating-transition-study-12336975, April 2022.
  * [21] H. Leland.  Corporate debt value, bond covenants, and optimal capital structure.  Journal of Finance, 49, 1994.
  * [22] H. Leland, R. Goldstein, and N. Ju.  An EBIT-based model of optimal capital structure.  Journal of Business, 74, 2001.
  * [23] T. Lukashiv, Y. Litvinchuk, I.V. Malyk, A. Golebiewska, and P.V. Nazarov.  Stabilization of stochastic dynamical systems of a random structure with Markov switches and Poisson perturbations.  Mathematics, 11, 2023.
  * [24] J. Medhi.  Stochastic Processes.  New Age International Publishers, New Delhi, 2nd edition, 1984.
  * [25] R. Merton.  Lifetime portfolio selection under uncertainty: the continuous-time case.  The Review of Economics and Statistics, 51, 1969.
  * [26] R. Merton.  An intertemporal capital asset pricing model.  Econometrica, 41, 1973.
  * [27] E. Platen and R. Rendek.  Exact scenario simulation for selected multi-dimensional stochastic processes.  Communications on Stochastic Analysis, 3(3):443–465, 2009.
  * [28] M. Pokojovy and A. T. Anum.  A fast initial response approach to sequential financial surveillance.  In H. Han and E. Baker, editors, The Recent Advances in Transdisciplinary Data Science, pages 19–33, Cham, 2022. Springer Nature Switzerland.
  * [29] C. S. Pun and H. Y. Wong.  Robust investment?reinsurance optimization with multiscale stochastic volatility.  Insurance: Mathematics and Economics, 62:245?256, 2015.
  * [30] The Economist. Buttonwood (Online Column).  How to avoid a common investment mistake.  https://www.economist.com/finance-and-economics/2023/09/21/how-to-avoid-a-common-investment-mistake, September 2023.
  * [31] J. Wachter.  Portfolio and consumption decisions under mean-reverting returns: An exact solution for complete markets.  Journal of Financial And Quantitative Analysis, 37, 2002.
  * [32] J. Wei, Y. Shen, and Q. Zhao.  Portfolio selection with regime-switching and state-dependent preferences.  Journal of Computational and Applied Mathematics, 365:112361, 2020\.


# Radiotherapy Dosimetry: A Review on Open-Source Optimizer

Paul Dubois  TheraPanacea  MICS, CentraleSupélec  Institut du Cancer de
Montpellier

###### Abstract

Radiotherapy dosimetry plays a crucial role in optimizing treatment plans for
cancer patients. In this study, we investigate the performance of a dozen
standard state-of-the-art open-source optimizers for radiotherapy dosimetry.
Our evaluation includes the use of TGG119 benchmark cases as well as one real
case obtained from the Institute du Cancer de Montpellier (ICM). Among the
tested optimizers, Newton CG demonstrates the fastest convergence in terms of
the number of iterations. However, when considering the computation time per
iteration, LBFGS emerges as the most efficient optimizer. These findings shed
light on the performance of open-source optimizers for radiotherapy dosimetry,
aiding practitioners in selecting suitable optimization tools for efficient
treatment planning.

##  1 Introduction

Radiotherapy, a widely utilized intervention for cancer treatment, employs
ionizing radiation to eliminate malignant cells. Intensity-modulated radiation
therapy (IMRT) has emerged as a notable technique within radiotherapy, aiming
to deliver high radiation doses to tumors while minimizing exposure to healthy
surrounding tissues [4]. Traditional IMRT strategies typically employ a set
number of beams, often 5, 7, or 9, originating from various angles around the
patient, commonly distributed evenly [3]. Each beam’s intensity is modulated
to optimize the delivery of radiation doses to the tumor while reducing
exposure to healthy tissues. This approach surpasses the effectiveness of the
3D-conformal radiotherapy (3D-CRT) technique [10] [16] [21]. To facilitate
precise and efficient radiation delivery, a computer-controlled device called
the multi-leaf collimator (MLC) is utilized to shape the radiation beam
according to the contours of the tumor.

The effectiveness of a radiotherapy treatment plan relies on the optimization
procedure, which involves a series of steps aimed at ensuring the optimal
delivery of radiation in accordance with the prescribed guidelines of medical
practitioners. Typically, computer software is employed to facilitate the
optimization process, taking into account various factors such as patient
anatomy, the size and location of the tumor and organs, and the radiation
objectives defined by medical professionals.

While there has been comparison between commercial software[19] [13], the aim
of this paper is to investigate the best optimizer for this task among the
open source optimizers.

#### Pre-dose-optimization

The initial stage of the optimization process entails the creation of a
virtual representation of the patient’s anatomical structure using advanced
medical imaging modalities, such as computed tomography (CT) or magnetic
resonance imaging (MRI) scans. This model is subsequently utilized to
accurately determine the size and location of the tumor, as well as to
delineate the surrounding healthy tissues that necessitate protection from
radiation exposure. Following this, the radiation dose required for effective
treatment is established, typically based on dose-volume objectives defined by
physicians (e.g., ensuring that 95% of the planning target volume receives a
minimum dose of 75 Gy). Determining the appropriate dose takes into
consideration factors such as tumor characteristics, location, size, as well
as the patient’s medical history and overall health status. These essential
steps in the optimization process are carried out by medical professionals
with expertise in radiotherapy treatment planning.

#### Radiotherapy doses

The subsequent step entails the computation of the radiation dose distribution
within the patient’s volumetric anatomy. This is achieved by simulating a
particular configuration of the multi-leaf collimator (MLC) on the patient’s
body, utilizing the available medical imaging data. The resulting computed
dose represents a mapping from the three-dimensional volume of the patient’s
anatomy to a scalar value measured in Grays (Gy), which denotes the absorbed
radiation energy. In practical implementation, a discrete representation of
the dose distribution is utilized, wherein the dose is calculated for each
individual voxel comprising the patient’s anatomical structure.

#### Dose-Volume Histograms

Medical professionals have meticulously identified and delineated the
pertinent anatomical structures within the patient’s anatomy. This allows the
computation of dose-volume histograms (DVHs) for each structure, predicated on
a specified dose distribution. The dose-volume objectives are subsequently
represented as specific points on the DVH curve, which correspond to the
desired minimum or maximum dose constraints. These objectives delineate the
desired thresholds that should be upheld, with points on the DVH curve either
located above (for minimum dose constraints) or below (for maximum dose
constraints) the prescribed thresholds.

#### Dose Evaluation

Physicians employ multiple criteria to assess the quality of a radiation dose
administered during treatment. Initially, they scrutinize the three-
dimensional distribution of the dose across the patient’s anatomy, focusing on
the spatial allocation among different anatomical structures, as well as
identifying the presence, number, and locations of regions with excessive
radiation (referred to as ”hot spots”). Subsequently, physicians conduct a
thorough analysis of the dose-volume histograms (DVHs) to evaluate the degree
of compliance with predefined DVH objectives. This crucial evaluation step
aims to safeguard the adjacent healthy tissues from unnecessary radiation
exposure. By optimizing the treatment plan and meticulously assessing the
quality of the dose distribution, physicians strive to ensure the attainment
of the most favorable outcome for the patient.

##  2 Methods

To ensure the precision and effectiveness of radiation therapy, a robust dose
optimization process are essential.

###  2.1 Radiotherapy Dosimetry

Radiotherapy dose optimization can be conceptualized as an inverse problem,
whereby the objective is to determine the most suitable radiation dose
distribution that aligns with the desired treatment outcome [22]. In other
terms, the challenge lies in identifying the radiation intensity or fluence
maps that deliver the prescribed dose to the tumor while minimizing exposure
to healthy tissues.

Mathematical optimization algorithms are employed to address the inverse
problem of radiotherapy dose optimization. These algorithms aim to identify
the optimal solution by minimizing a predefined objective function that
encompasses treatment goals and constraints. Typically, the objective function
includes terms that penalize both underdosing and overdosing of the tumor and
overdosing of surrounding healthy tissues. It may also incorporate terms that
account for the complexity or deliverability of the treatment plan.

Efficiently solving this optimization problem often involves designing the
objective function to be convex, thereby providing a well-defined target for
the optimization process. Gradient-based methods, Newtonian algorithms, or
quasi-Newtonian algorithms are commonly employed for this purpose. We aim at
benchmarking state-of-the-art open-source optimization algorithms for the
specific task of radiotherapy dosimetry.

###  2.2 Data

In this research endeavor, our focus was to evaluate the various open-source
optimizers. We used the widely recognized TG-119 [14] cases as a benchmark for
evaluating radiation therapy plans optimization. The TG-119 dataset provides
specific dose goals, which we incorporated into our proposed cost function.

We also used one real case of prostate cancer treatment from ICM. For this
case, doctors had provided specific dose goals, that we again incorporated
into our proposed cost function.

The TGG 119 multiple PTVs is a theoretical case, unlikely to happen in real
life. However, the three other cases represent a comprehensive set of what
dosimetrists could encounter on a daily basis.

The simulation of the beams was done using TheraPanacea dose engine, which
uses collapse cone convolution techniques, and is conformal to other simulator
available on the market.

###  2.3 Objective function

The cost function is formulated as a weighted sum of multiple objectives, with
each objective corresponding to a specific dose goal. The formulation is as
follows:

$f(\mathbf{d})=\sum_{o\in\mathcal{O}}w_{o}f_{o}(\mathbf{d})$

where:

  * •

$\mathbf{d}$ represents the dose distribution at the voxel level, and
$\mathbf{d}[s]$ denotes the dose on voxels within the structure $s$

  * •

$\mathcal{O}$ denotes the set of objectives corresponding to dose volume goals

  * •

$w_{o}$ signifies the weight assigned to the objective $o\in\mathcal{O}$

  * •

$o_{s}$, $o_{d}$, and $o_{v}$ refer to the structure, dose, and volume goals
of the objective $o\in\mathcal{O}$

The objective function $f_{o}(\mathbf{d})$ is computed based on the specific
type of dose volume constraint: If $o$ represents a maximum dose volume
constraint111e.g.: top 20% of the volume should receive at most 30 Gy,
$f_{o}(\mathbf{d})$ is calculated as
$\sum_{d\in\mathbf{d}[o_{s}]}(d-o_{d})_{+}^{2}$; if $o$ represents a minimum
dose volume constraint222e.g.: 95% of the volume should receive at least 70
Gy, then $f_{o}(\mathbf{d})$ is calculated as
$\sum_{d\in\mathbf{d}[o_{s}]}(o_{d}-d)_{+}^{2}$. The formulation involves a
squared over/under-dose penalty function.

In addition to the above, we introduced a regularization term that penalizes
variations in bixel values between neighboring regions, also employing a
squared penalty.

The optimization process involves finding the optimal bixel values
($\mathbf{b}$) by solving $\mathbf{d}=\mathbf{L}\mathbf{b}$, where
$\mathbf{L}$ is a precomputed dose-influence matrix mapping bixels to voxels.
Notably, since negative energy rays are physically infeasible, we ensured that
each bixel value is non-negative ($b\geq 0\quad\forall b\in\mathbf{b}$). To
achieve this, we computed $\mathbf{d}=\mathbf{L}|{\mathbf{b}}|$, where
$|{\mathbf{b}}|$ denotes the element-wise absolute value of $\mathbf{b}$.

By construction, the objective function is convex. Consequently, minimizing
the objective function with a given set of weights should invariably converge
to the same radiotherapy plan. To generate different treatment doses for the
same patient case, dosimetrists can play withe the weights of each sub-
objective function. This is outside the scope of this small review article, so
we decided to just set the weights of all constraints equal to one.

###  2.4 Open-source Optimizers

We tried to have a comprehensive test of available open-source optimizers,
here is a short description of the ones tested:

#### (Stochastic) Gradient Descent

Is an optimization algorithm that iteratively updates the model parameters in
the direction of the negative gradient of the objective function. In our case,
it is not stochastic, since it calculates the gradient using the current
solution333Our objective function has all its inputs as parameters, so there
is no notion of stochasticity. [8].

#### Conjugate Gradient

Is an iterative optimization algorithm commonly used to solve systems of
linear equations or quadratic optimization problems. It iteratively computes
conjugate directions and updates the solution along these directions, aiming
to minimize the objective function [6]. Conjugate Gradient is often applied in
scenarios where the Hessian matrix is unavailable or expensive to compute.

#### Newton

Newton’s method is an iterative optimization algorithm that uses the second-
order derivative (Hessian matrix) to find the minimum of a function. It
updates the current estimate by taking into account both the first-order
derivative (gradient) and the second-order derivative [15].

#### SLSQP

(Sequential Least Squares Programming) is a sequential quadratic programming
algorithm used for constrained optimization. It iteratively solves a sequence
of quadratic programming subproblems to find the optimal solution subject to
constraints [2].

#### RMSprop

(Root Mean Square Propagation) is an optimization algorithm that addresses the
problem of diminishing learning rates in traditional gradient descent methods.
It divides the learning rate by the root mean square of the past gradients,
which helps to stabilize and speed up convergence [7].

#### BFGS-based

##### Pure BFGS

(Broyden-Fletcher-Goldfarb-Shanno) is a quasi-Newton method that approximates
the Hessian matrix using updates based on gradient information. It performs a
line search to determine the step size that minimizes the objective function
along the search direction [5].

##### L-BFGS

(Limited-memory BFGS) is a variation of BFGS that uses a limited-memory
approach to approximate the Hessian matrix. It stores a limited number of past
gradient and parameter values to compute an approximate inverse Hessian matrix
efficiently [11].

#### Adam-based

##### Pure Adam

(Adaptive Moment Estimation) is an optimization algorithm that combines ideas
from both adaptive learning rates and momentum methods. It computes adaptive
learning rates for each parameter based on estimates of the first and second
moments of the gradients [9].

##### RAdam

(Rectified Adam) is a variant of the Adam optimizer that introduces a
rectification term to stabilize the adaptive learning rate. It aims to address
some convergence issues that can occur in Adam by dynamically adjusting the
variance of the adaptive learning rate [12].

##### NAdam

(Nesterov Adam) combines the Nesterov accelerated gradient method with the
Adam optimizer. It incorporates Nesterov momentum into the Adam update rule to
improve convergence and provide better generalization [20].

##### AdamDelta

Is another variant of the Adam optimizer that replaces the second moment
estimates (variance) with a delta parameter. It eliminates the need for
storing and updating the moving average of the squared gradients, which can be
beneficial in memory-constrained settings [23].

##### Adamax

Is an extension of the Adam optimizer that uses the infinity norm (max norm)
of the gradients instead of the L2 norm. It is designed to handle sparse
gradients more effectively and can be particularly useful in deep learning
models [1].

#### Rprop

(Resilient Backpropagation) is an optimization algorithm specifically designed
for neural networks. It updates the weights based on the sign of the gradient,
adjusting the step size adaptively. Rprop performs weight updates
independently for each weight parameter [18].

#### Other optimizers variations

In addition, we also tested AdamW, Adagrad and ASGD. However, AdamW & Adagrad
behaved similarly to Adam, and ASGD behaved similarly to SGD. For readability
purposes, we did not include them in the results.

##  3 Results

Figure 1: TGG 119: Multiple PTVs

Figure 2: TGG 119: Head and Neck

Figure 3: TGG 119: Prostate

Figure 4: ICM: (Typical) Prostate

#### Newton’s method

Based on the iterations-wise graph analysis, Newton’s method emerges as the
most optimal, consistently achieving a stable converged state within a mere 10
steps across all four examined cases. However, Newton’s method steps are very
expensive to compute, since it uses second order derivative (the Hessian) that
is expensive to compute (in terms of calculation time).

It is widely recognized that Newton’s method excels in optimizing convex
functions [17]. Our objective function is built to be convex. Hence, it make
sense that this optimization algorithm is particularly effective.

#### Best Algorithms

Other than Newton’s method, three algorithms have similar performances: Adam,
RMSprop and LBFGS. Adam (and RMSprop, to a lesser extent) appear to have more
”wavy” cost curves, while LBFGS cost decreases in a more stable way. These
observations are true both iteration and time-wise.

TGG 119 Multiple PTVs (figure 1) is the smallest problem, and th real ICM
prostate case (fig. 4) is the largest problem (in terms of
patient/organs/structure volume size); TGG 119 fake head & neck (fig. 2) and
TGG 119 fake prostate (fig. 3) have similar sizes. Notably, there is an
observable trend indicating that as the problem size increases, LBFGS
outperforms both RMSprop and Adam optimization algorithms.

#### LBFGS vs BFGS

It would be expected that BFGS performs better than LBFGS in terms of
iterations, but not in terms of time (since LBFGS is a fast approximation of
the BFGS technique). However, we observe that LBFGS outperforms BFGS even on
the iterations-wise graph. This suggests that the limited memory approximation
made are biased towards suitable directions in these type of problems.

##  4 Discussion

If it was possible to make Newton’s method faster, than we would advise to use
Newton’s optimization algorithm. However, to the best of our knowledge,
computing the Hessian remains long, not only in our implementation.

Hence, we advise to use the LBFGS algorithm for the problem of dose
optimization i radiotherapy; it appears to be the fastest to converge, and
converged steadily on the four cases tested.

## References

  * [1] Somenath Bera and Vimal K Shrivastava.  Analysis of various optimizers on deep convolutional neural network model in the application of hyperspectral remote sensing image classification.  International Journal of Remote Sensing, 41(7):2664–2683, 2020\.
  * [2] JF Bonnans, Jean Charles Gilbert, C Lemaréchal, and Claudia Sagastizabal.  Numerical Optimization – Theoretical and Practical Aspects.  01 2006.
  * [3] Thomas Bortfeld.  Imrt: a review and preview.  Physics in Medicine and Biology, 51(13):R363, jun 2006.
  * [4] Gary A. Ezzell, James M. Galvin, Daniel Low, Jatinder R. Palta, Isaac Rosen, Michael B. Sharpe, Ping Xia, Ying Xiao, Lei Xing, and Cedric X. Yu.  Guidance document on delivery, treatment planning, and clinical implementation of imrt: Report of the imrt subcommittee of the aapm radiation therapy committee.  Medical Physics, 30(8):2089–2115, 2003.
  * [5] Roger Fletcher and Chengian Xu.  Hybrid methods for nonlinear least squares.  IMA Journal of Numerical Analysis, 7(3):371–389, 1987.
  * [6] Magnus R. Hestenes and Eduard Stiefel.  Methods of conjugate gradients for solving linear systems.  Journal of research of the National Bureau of Standards, 49:409–435, 1952.
  * [7] Geoffrey Hinton, Nitish Srivastava, and Kevin Swersky.  Neural networks for machine learning lecture 6a overview of mini-batch gradient descent.  Cited on, 14(8):2, 2012.
  * [8] J.B. Hiriart-Urruty and C. Lemaréchal.  Fundamentals of Convex Analysis.  Grundlehren Text Editions. Springer Berlin Heidelberg, 2012.
  * [9] Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization, 2017.
  * [10] Thomas P. Kole, Osarhieme Aghayere, Jason Kwah, Ellen D. Yorke, and Karyn A. Goodman.  Comparison of heart and coronary artery doses associated with intensity-modulated radiotherapy versus three-dimensional conformal radiotherapy for distal esophageal cancer.  International Journal of Radiation Oncology*Biology*Physics, 83(5):1580–1586, 2012.
  * [11] Dong C. Liu and Jorge Nocedal.  On the limited memory bfgs method for large scale optimization.  Mathematical Programming, 45(1):503–528, Aug 1989.
  * [12] Liyuan Liu, Haoming Jiang, Pengcheng He, Weizhu Chen, Xiaodong Liu, Jianfeng Gao, and Jiawei Han.  On the variance of the adaptive learning rate and beyond, 2021.
  * [13] Erick Mora-Ramirez, Lore Santoro, Emmanuelle Cassol, Juan C Ocampo-Ramos, Naomi Clayton, Gunjan Kayal, Soufiane Chouaf, Dorian Trauchessec, Jean-Pierre Pouget, Pierre-Olivier Kotzki, Emmanuel Deshayes, and Manuel Bardiès.  Comparison of commercial dosimetric software platforms in patients treated with 177 Lu‐DOTATATE for peptide receptor radionuclide therapy.  Med. Phys., 47(9):4602–4615, September 2020.
  * [14] Dinesh Kumar Mynampati, Ravindra Yaparpalvi, Linda Hong, Hsiang-Chi Kuo, and Dennis Mah.  Application of aapm tg 119 to volumetric arc therapy (vmat).  Journal of Applied Clinical Medical Physics, 13(5):108–116, 2012\.
  * [15] Jorge Nocedal and Stephen J. Wright.  Numerical optimization.  Springer series in operations research and financial engineering. Springer, New York, NY, 2nd ed. edition, 2006.  MAB0014.001: M 11.0268.
  * [16] David Palma, Emily Vollans, Kerry James, Sandy Nakano, Vitali Moiseenko, Richard Shaffer, Michael McKenzie, James Morris, and Karl Otto.  Volumetric modulated arc therapy for delivery of prostate radiotherapy: Comparison with intensity-modulated radiotherapy and three-dimensional conformal radiotherapy.  International Journal of Radiation Oncology*Biology*Physics, 72(4):996–1001, 2008.
  * [17] Barnabas Poczos and Ryan Tibshirani.  Convex Optimization.
  * [18] Martin A. Riedmiller and Heinrich Braun.  Rprop - a fast adaptive learning algorithm.  1992\.
  * [19] Bin Tang, Xiangyu Liu, Xianliang Wang, Shengwei Kang, Pei Wang, Jie Li, and Lucia Clara Orlandini.  Dosimetric comparison of graphical optimization and inverse planning simulated annealing for brachytherapy of cervical cancer.  J Contemp Brachytherapy, 11(4):379–383, August 2019.
  * [20] Ange Tato and Roger Nkambou.  Improving adam optimizer.  2018\.
  * [21] Lynn J. Verhey.  Comparison of three-dimensional conformal radiation therapy and intensity-modulated radiation therapy systems.  Seminars in Radiation Oncology, 9(1):78–98, 1999.  Radiation Therapy Treatment Optimization.
  * [22] S Webb.  The physical basis of imrt and inverse planning.  The British Journal of Radiology, 76(910):678–689, 2003.  PMID: 14512327.
  * [23] Matthew D. Zeiler.  Adadelta: An adaptive learning rate method, 2012.