# Subleading correction to the Asian options volatility in the Black-Scholes
model

Dan Pirjol  School of Business
Stevens Institute of Technology
Hoboken, NJ 07030<EMAIL_ADDRESS>
(Date: 5 July 2024)

###### Abstract.

The short maturity limit $T\to 0$ for the implied volatility of an Asian
option in the Black-Scholes model is determined by the large deviations
property for the time-average of the geometric Brownian motion. In this note
we derive the subleading $O(T)$ correction to this implied volatility, using
an asymptotic expansion for the Hartman-Watson distribution. The result is
used to compute subleading corrections to Asian options prices in a small
maturity expansion, sharpening the leading order result obtained using large
deviations theory. We demonstrate good numerical agreement with precise
benchmarks for Asian options pricing in the Black-Scholes model.

###### Key words and phrases:

Asymptotic expansions, Asian options, time integral of the geometric Brownian
motion

##  1\. Introduction

Asian options are derivatives with payoff linked to the time average of the
asset price

(1) $A_{T}:=\frac{1}{T}\int_{0}^{T}S_{t}dt\,.$

We are interested in pricing Asian options under the Black-Scholes model where
the asset price follows a geometric Brownian motion

(2) $\frac{dS_{t}}{S_{t}}=(r-q)dt+\sigma dW_{t}\,,$

with initial condition $S_{0}>0$. Asian options pricing under the Black-
Scholes model has been widely studied, using and a large variety of
approaches, both numerical and analytical. See [8] for a survey of methods.
Restricting to analytical approaches we mention here the Geman-Yor approach
[9, 5], the Laguerre polynomial expansion method [3], the PDE expansion method
[22, 7], and the spectral method [12]. We mention also the large- and small-
strike asymptotics of Asian options in the Black-Scholes model obtained in
[11] and [24].

The short maturity asymptotics of Asian option prices has been studied using
probabilistic methods from Large Deviations theory [2, 16, 19], assuming that
$S_{t}$ follows a one-dimensional diffusion

(3) $dS_{t}=\sigma(S_{t})S_{t}dW_{t}+(r-q)S_{t}dt\,,$

under suitable technical conditions on the volatility function
$\sigma(\cdot)$. These results include as a limiting case the Black-Scholes
model corresponding to $\sigma(x)=\sigma$. The short maturity asymptotics of
Asian options under stochastic volatility models has been studied recently
using Malliavin calculus methods in [1].

The leading short maturity asymptotics of the out-of-the-money Asian option
prices is given in Theorem 2 in [16]. We quote the result applied to the
Black-Scholes model. For simplicity of notation we take $q=0$ in the
following.

######  Theorem 1.

Assume that the asset price follows the Black-Scholes model
$S_{t}=S_{0}e^{\sigma W_{t}+(r-\frac{1}{2}\sigma^{2})t}$.

(i) for out-of-the-money call Asian options we have

(4) $\lim_{T\to 0}T\log
C(K,T)=-\frac{1}{\sigma^{2}}J_{BS}(K/S_{0})\,,K>S_{0}\,.$

(ii) for out-of-the-money put Asian options we have

(5) $\lim_{T\to 0}T\log
P(K,T)=-\frac{1}{\sigma^{2}}J_{BS}(K/S_{0})\,,K<S_{0}\,.$

The rate function $J_{BS}(k)$ is given in closed form in Proposition 12 in
[16]

(6)
$J_{BS}(k)=\left\\{\begin{array}[]{cc}\frac{1}{2}\beta^{2}-\beta\tanh\frac{\beta}{2}&\,,k\geq
1\\\ 2\xi(\tan\xi-\xi)&\,,0<k\leq 1\\\ \end{array}\right.$

where $\beta$ is the solution of the equation $\frac{\sinh\beta}{\beta}=k$ and
$\xi$ is the solution in $[0,\frac{\pi}{2}]$ of the equation $\frac{\sin
2\xi}{2\xi}=k$.

The analyticity properties of the function $J_{BS}(z)$ in the complex $z$
plane were studied in [13], see Sec. 4.1. This function has no singularities
along the real positive axis, and the closest singularity to $z=1$ is a pole
at $z=0$. For practical computations it is convenient to use the Taylor
expansion of the rate function in powers of $\log k$. The first few terms of
this expansion are

(7)
$J_{BS}(k)=\frac{3}{2}\log^{2}k-\frac{3}{10}\log^{3}k+\frac{109}{1{\small,}400}\log^{4}k+O(\log^{5}k)\,.$

The radius of convergence of this series is determined by the position of the
singularities in the complex plane of the function $z\mapsto J_{BS}(e^{z})$
and is $|\log k|<3.49295$, see Proposition 4.1(ii) in [13]. Outside of this
region the exact result (6) must be used.

An Asian option with maturity $T$ and strike $K$ can be priced as an European
option on a Black-Scholes asset, with the same maturity and strike, with an
implied volatility $\Sigma_{\rm LN}(K,T)$ chosen such that

(8) $C(K,T)=C_{BS}(K,T;A_{\rm fwd},\Sigma_{\rm LN}(K,T))$

where $C(K,T)$ is the Asian option price, and the forward price $A_{\rm fwd}$
is given by

(9) $A_{\rm
fwd}:=\frac{1}{T}\int_{0}^{T}\mathbb{E}[S_{t}]dt=S_{0}\frac{e^{(r-q)T}-1}{(r-q)T}\,.$

We call the volatility $\Sigma_{\rm LN}(K,T)$ the equivalent log-normal
volatility of the Asian option, following Sec. 4.3 of [17]. The representation
(8) is also useful for computing Asian option sensitivities [18].

The short-maturity asymptotics for the Asian option prices of Proposition 1 is
equivalent with a short-maturity asymptotics for the equivalent log-normal
volatility

(10) $\lim_{T\to 0}\Sigma^{2}_{\rm
LN}(K,T)=\sigma^{2}\frac{\log^{2}(K/S_{0})}{2J_{\rm
BS}(K/S_{0})}=:\Sigma_{0}^{2}(K/S_{0})\,.$

Using the expansion (7) we get

(11) $\Sigma_{0}^{2}(k)=\sigma^{2}\frac{\log^{2}k}{2J_{\rm
BS}(k)}=\frac{1}{3}\sigma^{2}\left(1+\frac{1}{5}\log
k-\frac{1}{84}\log^{2}k-\frac{17}{10,\\!500}\log^{3}k+O(\log^{4}k)\right)$

where $k=K/S_{0}$ and $J_{BS}(k)$ is the rate function appearing in the
statement of Proposition 1. See also Proposition 18 in [16]. A similar
asymptotic result was obtained in [16] in the more general setting of the
local volatility model, and in [20] for a class of jump-diffusion models with
local volatility.

While our study of the equivalent log-normal volatility $\Sigma_{LN}(K,T)$ is
limited to a short maturity expansion, we note an exact prediction which can
be extracted from the results of Ref. [24]. Proposition 1 in this paper
implies the extreme strikes asymptotics
$\lim_{K\to\\{0,\infty\\}}\Sigma_{LN}(K,T)=\sigma$, for any $T>0$.

The equivalent log-normal variance can be expanded in a short maturity
expansion as

(12)
$\Sigma^{2}_{LN}(K,T)=\Sigma_{0}^{2}(e^{x})+2\Sigma_{0}(e^{x})\Sigma_{1}(K,T)T+O(T^{2})\,.$

When including the higher order corrections in $T$, it will be seen to be
convenient to work with the log-moneyness parameter $x=\log(K/A_{\rm fwd})$.
This is expanded as $x=\log k+O(rT)$, such that at leading order in $T$ the
log-moneyness and log-strike are equivalent. When working at higher orders in
$T$ it is important to keep track of the higher order terms in this relation.

In this paper, we compute the $O(T)$ term in the expansion (12). This
correction can be expanded in powers of log-moneyness as

(13)
$\frac{1}{\sigma^{2}}2\Sigma_{0}(e^{x})\Sigma_{1}(K,T)T=(\sigma^{2}T)(s_{0}+s_{1}x+O(x^{2}))+(rT)(r_{0}+r_{1}x+O(x^{2}))\,.$

We give explicit results for the first two terms in this expansion
$r_{0,1},s_{0,1}$. The higher order terms in the $x$-expansion can be
evaluated using the same approach.

We summarize the expansion of the Asian implied volatility including the
$O(T)$ term in the following result. This is the main result of this paper.

######  Proposition 2.

Assume that the asset price follows the Black-Scholes model $dS_{t}=\sigma
S_{t}dW_{t}+rS_{t}dt$. The equivalent log-normal variance of an Asian option
with strike $K$ and maturity $T$ is

$\displaystyle\Sigma_{\rm LN}^{2}(K,T)$ $\displaystyle=$
$\displaystyle\sigma^{2}\Big{\\{}\underbrace{\frac{\log^{2}x}{2J_{\rm
BS}(x)}}\underbrace{-\frac{61}{9,\\!450}(\sigma^{2}T)+\frac{1}{12}(rT)}$
$\displaystyle\hskip 28.45274ptO(1)\hskip 42.67912ptO(T)$ $\displaystyle\hskip
14.22636pt+\underbrace{\Big{[}-\frac{34}{23,\\!625}(\sigma^{2}T)\Big{]}x}+O(Tx^{2})+O(T^{2})\Big{\\}}\,.$
$\displaystyle\hskip 71.13188ptO(Tx)$

where $x=\log(K/A_{\rm fwd})$ is the option log-moneyness.

In particular, for an at-the-money Asian option with strike $K=A_{\rm fwd}$ we
have

(15) $\Sigma^{2}_{\rm LN}(K=A_{\rm
fwd},T)=\sigma^{2}\left(\frac{1}{3}-\frac{61}{9,\\!450}(\sigma^{2}T)+\frac{1}{12}(rT)+O(T^{2})\right)\,.$

###  1.1. Standardization

It is convenient to standardize the pricing problem by reducing it to the
study of the distributional properties of the quantity

(16) $A_{t}^{(\mu)}=\int_{0}^{t}e^{2(B_{s}+\mu s)}ds\,,$

where $B_{t}$ is a standard Brownian motion. Using the rescaling property of
the Brownian motion $B_{\lambda t}=\sqrt{\lambda}B_{t}$, the time-average of
the geometric Brownian motion in (1) can be expressed in terms of
$A_{t}^{(\mu)}$ as [9]

(17)
$A_{T}=\frac{4S_{0}}{\sigma^{2}T}A_{\frac{1}{4}\sigma^{2}T}^{(\frac{2r}{\sigma^{2}}-1)}=S_{0}\frac{1}{\tau}A_{\tau}^{(\mu)}\,,\quad\tau=\frac{1}{4}\sigma^{2}T\,,\quad\mu=\frac{2r}{\sigma^{2}}-1\,.$

The prices of fixed strike Asian options with averaging over the period
$[0,T]$ and strike $K$ can be expressed in terms of the standardized average
$\frac{1}{\tau}A_{\tau}^{(\mu)}$ as [9]

(18) $\displaystyle
C(K,T)=e^{-rT}\mathbb{E}\left[\left(\frac{1}{T}A_{T}-K\right)^{+}\right]=S_{0}e^{-rT}c(k,\tau)$
(19) $\displaystyle
P(K,T)=e^{-rT}\mathbb{E}\left[\left(K-\frac{1}{T}A_{T}\right)^{+}\right]=S_{0}e^{-rT}p(k,\tau)\,.$

with $k=K/S_{0}$ and

(20)
$c(k,\tau):=\mathbb{E}\Big{[}\Big{(}\frac{1}{\tau}A_{\tau}^{(\mu)}-k\Big{)}^{+}\Big{]}\,,\quad
p(k,\tau):=\mathbb{E}\Big{[}\Big{(}K-\frac{1}{\tau}A_{\tau}^{(\mu)}\Big{)}^{+}\Big{]}\,.$

The normalized Asian options $c(k,\tau),p(k,\tau)$ correspond to volatility
$\sigma=2$ and drift $r=\mu+1$, with $S_{0}=1$. The equivalent log-normal
volatility for the normalized options $\sigma_{LN}(k,\tau)$ is defined as

(21) $c(k,\tau)=C_{BS}(k,\tau;a_{\rm fwd}^{(\mu)},\sigma_{\rm
LN}(k,\tau))\,,\quad a_{\rm
fwd}^{(\mu)}:=\mathbb{E}\Big{[}\frac{1}{\tau}A_{\tau}^{(\mu)}\Big{]}=\frac{e^{2(\mu+1)\tau}-1}{2(\mu+1)\tau}\,.$

This is rescaled to the general $S_{0},r,\sigma$ BS model as

(22)
$\Sigma_{LN}^{2}(K,T;S_{0},\sigma)=\frac{1}{4}\sigma^{2}\sigma_{LN}^{2}\Big{(}\frac{K}{S_{0}},\frac{1}{4}\sigma^{2}T\Big{)}\,.$

###  1.2. Outline

The paper is structured as follows. The starting point is the asymptotic
expansion of the density of the time average of the gBM
$\frac{1}{t}A_{t}^{(\mu)}$ given in Proposition 6 in [21], obtained from a
small time expansion of the Hartman-Watson distribution. We collect the
relevant properties of this expansion in Section 2.

The proof of the main result, Proposition 2, is given in Section 3 and is
divided into three parts, organized as separate sections. In Section 3.1 we
compute the subleading correction to the price of Asian options in the BS
model in a small maturity expansion. This is used in Sections 3.2 and 3.3 to
obtain the $O(T)$ correction to the equivalent log-normal implied volatility
of an Asian option by an application of the Gao-Lee transfer result [10]. The
result can be used as the basis of a simple numerical pricing approximation
for Asian options in the Black-Scholes model. In Section 4 we present
numerical tests on benchmark cases in the literature, showing that adding the
subleading correction improves the agreement of the asymptotic expansion with
the benchmark evaluations.

##  2\. The asymptotic distribution of $\frac{1}{t}A_{t}^{(\mu)}$

The starting point for our analysis is a result obtained in [21] for the
leading asymptotics of the density of $\frac{1}{t}A_{t}^{(\mu)}$ as $t\to 0$.
This was obtained by applying Laplace asymptotic methods to a one-dimensional
integral giving this density (due to Yor [23]), combined with an asymptotic
expansion for the Hartman-Watson distribution $\theta(r,t)$ as $t\to 0$.

For completeness, we summarize the main results of [21]. Denote the density of
the normalized average of the geometric Brownian motion (gBM)

(23) $\mathbb{P}\Big{(}\frac{1}{\tau}A_{\tau}^{(\mu)}\in
da\Big{)}=f(a,\tau)\frac{da}{a}\,.$

The density is expressed as [23]

(24)
$f(a,\tau)=e^{-\frac{1}{2}\mu^{2}\tau}a^{\mu-1}\int_{0}^{\infty}\rho^{\mu}e^{-\frac{1+a^{2}\rho^{2}}{2a\tau}}\theta(\rho/\tau,\tau)\frac{d\rho}{\rho}$

The Hartman-Watson function is defined by the integral

(25)
$\theta(r,t)=\frac{r}{\sqrt{2\pi^{3}t}}e^{\frac{\pi^{2}}{2t}}\int_{0}^{\infty}e^{-\frac{\xi^{2}}{2t}}e^{-r\cosh\xi}\sinh\xi\sin\frac{\pi\xi}{t}d\xi$

Proposition 1 of [21] gives an expansion for this function as $t\to 0$ at
fixed $\rho=rt$

(26) $\theta(\rho/t,t)=\frac{1}{2\pi
t}e^{-\frac{1}{t}[F(\rho)-\frac{\pi^{2}}{2}]}G(\rho)(1+\vartheta(\rho,t))$

where the functions $F(\rho),G(\rho)$ are known in closed form, and the error
term is bounded as $|\vartheta(\rho,t)|\leq\frac{1}{70}t$.

The density of the time integral of the gBM $f(a,\tau)$ can be approximated
with the (properly normalized) leading term of this expansion as
$f(a,\tau)=f_{0}(a,\tau)(1+\varepsilon(a,\tau))$ with

(27) $f_{0}(a,\tau):=\frac{1}{n(\tau)}f_{\rm HW}(a,\tau)$

where

(28) $f_{HW}(a,\tau):=\frac{1}{2\pi
t}e^{-\frac{1}{2}\mu^{2}t}a^{\mu-1}\int_{0}^{\infty}\rho^{\mu}G(\rho)e^{-\frac{1}{\tau}H(\rho,a)}\frac{d\rho}{\rho}\,.$

We denoted here

(29) $H(\rho)=\frac{1+a^{2}\rho^{2}}{2a}-\frac{\pi^{2}}{2}+F(\rho)$

and

(30) $n(\tau)=\int_{0}^{\infty}f_{HW}(a,\tau)\frac{da}{a}$

is a normalization factor which ensures that $f_{0}(a,\tau)$ is normalized as
$\int_{0}^{\infty}f_{0}(a,\tau)\frac{da}{a}=1$.

The application of Laplace asymptotic methods to the integral (28) gives the
following result.

######  Proposition 3.

[Proposition 6 in [21]] We have the $\tau\to 0$ asymptotics

(31) $\displaystyle f_{HW}(a,\tau)=\frac{1}{\sqrt{2\pi
t}}g(a,\mu)e^{-\frac{1}{t}J(a)}\frac{da}{a}(1+O(\tau))$

where

(32)
$g(a,\mu):=(a\rho_{*})^{\mu}G(\rho_{*})\frac{1}{\sqrt{H^{\prime\prime}(\rho_{*})}}\frac{1}{\rho_{*}}\,.$

We denote here $J(a)\equiv\inf_{\rho\geq 0}H(\rho)=H(\rho_{*})$ and
$\rho_{*}=\mbox{argmin}H(\rho)$. From (29) it follows that the minimizer
$\rho_{*}$ depends only on $a$ but not on $\mu$.

Substituting the result (31) into (27) gives an approximation $f_{0}(a,\tau)$
for the density of the time integral of the gBM. The error of this
approximation is bounded by the following result.

######  Proposition 4.

The error of the approximation (27) is bounded as

(33) $f_{0}(a,\tau)\frac{-\frac{1}{35}\tau}{1+\frac{1}{70}\tau}\leq
f(a,\tau)-f_{0}(a,\tau)\leq
f_{0}(a,\tau)\frac{\frac{1}{35}\tau}{1-\frac{1}{70}\tau}\,.$

###### Proof.

Using (24) we have

(34) $\displaystyle|f(a,\tau)-f_{HW}(a,\tau)|\leq
e^{-\frac{1}{2}\mu^{2}\tau}a^{\mu-1}\int_{0}^{\infty}\rho^{\mu}e^{-\frac{1+a^{2}\rho^{2}}{2a\tau}}|\theta(\rho/\tau,\tau)-\frac{1}{2\pi\tau}e^{-\frac{1}{t}[F(\rho)-\frac{\pi^{2}}{2}]}G(\rho)|\frac{d\rho}{\rho}$
$\displaystyle\hskip
56.9055pt\leq\frac{1}{2\pi\tau}e^{-\frac{1}{2}\mu^{2}t}a^{\mu-1}\int_{0}^{\infty}\rho^{\mu}G(\rho)e^{-\frac{1}{\tau}H(\rho,a)}|\vartheta(\rho,\tau)|\frac{d\rho}{\rho}\leq\frac{1}{70}\tau
f_{HW}(a,\tau)$

In the last step we used the error bound
$|\vartheta(\rho,\tau)|\leq\frac{1}{70}\tau$.

In a similar way we have

(35)
$\displaystyle|1-n(\tau)|=\left|\int_{0}^{\infty}(f(a,\tau)-f_{HW}(a,\tau))\frac{da}{a}\right|\leq\int_{0}^{\infty}|f(a,\tau)-f_{HW}(a,\tau)|\frac{da}{a}$
$\displaystyle\qquad\leq\frac{1}{70}\tau\int_{0}^{\infty}f_{HW}(a,\tau)\frac{da}{a}=\frac{1}{70}\tau
n(\tau)\,,$

where we used (34) in the last step.

From these two inequalities we get

(36) $\frac{1}{1+\frac{1}{70}\tau}f(a,\tau)\leq
f_{HW}(a,\tau)\leq\frac{1}{1-\frac{1}{70}\tau}f(a,\tau)$

and

(37) $\frac{1}{1+\frac{1}{70}\tau}\leq
n(\tau)\leq\frac{1}{1-\frac{1}{70}\tau}$

Taking their ratio gives

(38) $\frac{1-\frac{1}{70}\tau}{1+\frac{1}{70}\tau}f(a,\tau)\leq
f_{0}(a,\tau)\leq\frac{1+\frac{1}{70}\tau}{1-\frac{1}{70}\tau}f(a,\tau)$

These inequalities can be inverted to give bounds for $f(a,\tau)$ in terms of
$f_{0}(a,\tau)$, which can be expressed as the error bounds (33).

∎

The function $f_{HW}(a,\tau)$ depends on two functions $J(a)$ and $g(a,\mu)$.
The properties of the function $J(a)$ were studied in Sec. 4.1 of [21] where
it was shown that it is simply related to the rate function $J_{BS}(k)$
appearing in the short-maturity asymptotics of Asian options, as

(39) $J(a)=\frac{1}{4}J_{BS}(a)\,.$

The following expansion of $g(a,\mu)$ was obtained in Proposition 10 of [21].
The coefficient $c_{3}$ quoted below is new.

######  Proposition 5.

The function $g(a,\mu)$ has the expansion

$\displaystyle g(a,\mu)$ $\displaystyle=$ $\displaystyle e^{\mu\log
a+(\mu-1)\log\rho_{*}(a)}G(\rho_{*})\frac{1}{\sqrt{H^{\prime\prime}(\rho_{*})}}$
$\displaystyle=$ $\displaystyle\frac{\sqrt{3}}{2}e^{c_{1}\log
a+c_{2}\log^{2}a+c_{3}\log^{3}a+O(\log^{4}a)}\,.$

The first few coefficients $c_{i}$ are

(41) $\displaystyle c_{1}=\frac{3}{4}(\mu+1)-\frac{4}{5}$ (42) $\displaystyle
c_{2}=-\frac{3}{80}(\mu+1)+\frac{57}{1,400}$ (43) $\displaystyle
c_{3}=\frac{1}{350}(\mu+1)-\frac{1}{875}\,.$

##  3\. Subleading corrections to the Asian implied volatility

The equivalent log-normal volatility $\Sigma_{\rm LN}(K,T)$ of an Asian option
in the Black-Scholes model can be expanded in powers of maturity as

(44) $\Sigma_{\rm
LN}(K,T)=\Sigma_{0}(K/S_{0})+T\Sigma_{1}(K,S_{0})+O(T^{2})\,.$

The leading term in this expansion is determined from the short-maturity
asymptotics of the Asian option prices [16]

(45) $\Sigma_{0}^{2}(k)=\sigma^{2}\frac{\log^{2}k}{2J_{\rm
BS}(k)}=\frac{1}{3}\sigma^{2}\left(1+\frac{1}{5}\log
k-\frac{1}{84}\log^{2}k-\frac{17}{10,500}\log^{3}k+O(\log^{4}k)\right)\,.$

See Proposition 18 in [16], where the equivalent log-normal volatility is
denoted $\Sigma_{\rm LN}(K,S_{0})$.

We compute here the subleading term of $O(T)$ to the equivalent log-normal
volatility. The proof proceeds in three steps. In the first step (Sec. 3.1) we
compute the short maturity asymptotics for the reduced Asian option prices. In
the second step (Sec. 3.2) we determine the equivalent log-normal volatility
in the driftless case $r=0$, and in the third step (Sec. 3.3) a non-zero
interest rate is added.

###  3.1. Short maturity asymptotics for Asian option prices

In this section we use the asymptotic distribution of the time average
$\frac{1}{t}A_{t}^{(\mu)}$ in Proposition 3 to compute the price of OTM Asian
options in the Black-Scholes model.

######  Proposition 6.

The leading asymptotics for the OTM Asian options with reduced strike
$k=K/S_{0}$ and maturity $\tau$ is

(46) $\displaystyle
c(k,\tau)=\sqrt{\frac{\tau^{3}}{2\pi}}\frac{g(k,\mu)}{k[J^{\prime}(k)]^{2}}e^{-\frac{1}{\tau}J(k)}(1+O(\tau))\,,\quad
k\geq 1$ (47) $\displaystyle
p(k,\tau)=\sqrt{\frac{\tau^{3}}{2\pi}}\frac{g(k,\mu)}{k[J^{\prime}(k)]^{2}}e^{-\frac{1}{\tau}J(k)}(1+O(\tau))\,,\quad
k\leq 1\,.$

###### Proof.

The reduced Asian option price is expressed as an integral over the exact
distribution of the time-integral of the gBM $f(a,\tau)$, defined in (23)

(48) $c(k,\tau)=\int_{0}^{\infty}(a-k)^{+}f(a,\tau)\frac{da}{a}\,.$

We derive an approximation for $c(k,\tau)$ by performing two successive
approximations:

i) replace $f(a,\tau)$ with its $\tau\to 0$ leading order approximation
$f_{0}(a,\tau)$ defined in (27). Define the corresponding approximation for
the option prices

(49)
$c_{0}(k,\tau):=\int_{0}^{\infty}(a-k)^{+}f_{0}(a,\tau)\frac{da}{a}\,,\quad
p_{0}(k,\tau):=\int_{0}^{\infty}(k-a)^{+}f_{0}(a,\tau)\frac{da}{a}\,.$

The error of this approximation is bounded using the error bound (33) as

(50) $-\frac{1}{35}\tau\frac{1}{1+\frac{1}{70}\tau}c_{0}(k,\tau)\leq
c(k,\tau)-c_{0}(k,\tau)\leq\frac{1}{35}\tau\frac{1}{1-\frac{1}{70}\tau}c_{0}(k,\tau)$

and analogous for $p(k,\tau)$. The approximation error is bounded in absolute
value as

(51)
$|c(k,\tau)-c_{0}(k,\tau)|\leq\frac{1}{35}\tau\frac{1}{1+\frac{1}{70}\tau}c_{0}(k,\tau)\,.$

The approximation $c_{0}(k,\tau)$ is expressed as a double integral with an
integrand known in closed form. This can be easily evaluated numerically, and
offers a simple approximation for pricing Asian options in the Black-Scholes
model, with controlled approximation error. Tests of this approach in Section
5 of [13] demonstrate good agreement with the precise benchmarks of [12].

ii) Next we compute the leading approximation for $c_{0}(k,\tau)$ as $\tau\to
0$ using standard Laplace asymptotic methods for integrals. The result we use
is due to Erdélyi, see Sec. 2.4 in [6], and appears as Theorem 8.1 in Olver
[15]. We give a few details of the application of this result to the Asian
call price asymptotics, using the notations of Theorem 1.2.1 of Nemes [14],
which is reproduced in the Appendix. The theorem applies to our case with the
substitutions:

(52) $\lambda\mapsto 1/\tau\,,\quad f(x)\mapsto J(a)\,,\quad
g(x)\mapsto\frac{1}{\sqrt{2\pi\tau}}(a-k)g(a,\mu)\frac{1}{a}$

and $\alpha=1,\beta=2$.

The technical conditions of the theorem are satisfied: i) The function $J(a)$
is increasing on the integration interval $[k,\infty)$ with $k>1$. ii)
$J(a),g(a,\mu)$ are continuous functions on $a\in[k,\infty)$. iii) The
functions $J(a),g(a,\mu)$ can be expanded around $k>1$ as in (88). For
$J_{BS}(x)$ this follows from the analyticity of this function for $x>0$, see
Sec. 4.1 in [13]. A similar result holds for $g(a,\mu)$ and follows from the
analyticity of $F(\rho),G(\rho)$ for real positive $\rho$ proved in Sec. 4.2
of [13]. The leading coefficients in the expansion (88) are

(53) $a_{0}=J^{\prime}(k)\,,\quad
b_{0}=\frac{1}{\sqrt{2\pi\tau}}g(k,\mu)\frac{1}{k}\,.$

iv) The integrals (49) converge. See Sec. 5 of [13] for numerical evaluations
of $c_{0}(k,\tau)$.

At leading order the Laplace asymptotic expansion (89) gives
$c_{0}(k,\tau)=c_{0}^{L}(k,\tau)(1+O(\tau))$ with

(54)
$c_{0}^{L}(k,\tau):=e^{-\frac{1}{\tau}J(k)}\frac{d_{0}}{(1/\tau)^{2}}=\sqrt{\frac{\tau^{3}}{2\pi}}\frac{g(k,\mu)}{k[J^{\prime}(k)]^{2}}e^{-\frac{1}{\tau}J(k)}$

where we used

(55)
$d_{0}=\frac{b_{0}}{a_{0}^{2}}=\frac{1}{\sqrt{2\pi\tau}}g(k,\mu)\frac{1}{k(J^{\prime}(k))^{2}}\,.$

The correction to the leading order term is of order
$|c_{0}(k,\tau)-c_{0}^{L}(k,\tau)|=c_{0}^{L}(k,\tau)(1+O(\tau))$.

The combined error of the two approximations is

(56)
$\displaystyle|c(k,\tau)-c_{0}^{L}(k,\tau)|\leq|c(k,\tau)-c_{0}(k,\tau)|+|c_{0}(k,\tau)-c_{0}^{L}(k,\tau)|$
$\displaystyle\qquad\leq
c_{0}(k,\tau)(1+O(\tau))+c_{0}^{L}(k,\tau)(1+O(\tau))=c_{0}^{L}(k,\tau)(1+O(\tau))\,.$

This reproduces the quoted result (46). The Asian put option result (47) is
obtained in a similar way. ∎

###  3.2. The driftless case $r=0$

We start with the simpler case $r=0$. Recall that in terms of the normalized
parameters introduced in Sec. 1.1 this corresponds to $\mu=-1$. In the next
step (Sec. 3.3) we include the contribution of a non-zero interest rate.

The starting point is the asymptotic result for option prices of Proposition
6. The leading asymptotic result for an OTM Asian call option has the form
$c(k,\tau)=\tau^{3/2}h(k)e^{-\frac{1}{\tau}J(k)}$ with
$h(k)=\frac{1}{\sqrt{2\pi}}\frac{g(k,\mu)}{k[J^{\prime}(k)]^{2}}$. Recall
$k=K/S_{0}$.

The small-$\tau$ asymptotics of the Asian option can be expressed as an
expansion for the log-price $L=-\log c(k,\tau)$. Using the notations of Gao
and Lee [10], the first terms of this asymptotics are

(57) $L=-\log c(k,\tau)=\frac{1}{\tau}J(k)-\frac{3}{2}\log\tau+\alpha_{0}(k)$

with $\alpha_{0}(k):=-\log h(k)$ which is expanded as

(58)
$\displaystyle\alpha_{0}(k)=-\log\left(\frac{16\sqrt{3}}{18\sqrt{2\pi}}\cdot\frac{k}{\log^{2}k}\right)+2\log\left(1-\frac{3}{10}\log
k+O(\log^{2}k)\right)-\sum_{i=1}^{\infty}c_{i}\log^{i}k$

The second term in this expression is the contribution from $J^{\prime}(k)$
which is expanded by substituting the series expansion (7) for $J_{BS}(k)$ and
differentiating term by term

(59) $J^{\prime}(k)=\frac{1}{4}J^{\prime}_{\rm BS}(k)=\frac{3}{4k}\log
k\left(1-\frac{3}{10}\log k+\frac{109}{1,050}\log^{2}k+O(\log^{3}k)\right)\,.$

The third term in (58) contains the contribution of the exponential factor for
$g(k,\mu)$ in Proposition 5, which is determined by the coefficients $c_{i}$
appearing in the expansion of the exponent around $k=1$.

By Corollary 7.4 in Gao, Lee [10], the asymptotic implied variance is

(60) $\sigma^{2}_{\rm
LN}(k,\tau)=\frac{\log^{2}k}{2J(k)}-\frac{\log^{2}k}{4J^{2}(k)}\left(\log
k+\log\frac{\log^{2}k}{16\pi}+2\alpha_{0}(k)-3\log
J(k)\right)\tau+O(\tau^{2})$

The expression in the brackets in the second term of (60) is expanded around
$k=1$ as

(61) $\displaystyle B(k):=\log
k+\log\frac{\log^{2}k}{16\pi}+2\alpha_{0}(k)-3\log J(k)$
$\displaystyle\quad\,\,=b_{0}+b_{1}\log
k+b_{1L}\log\log^{2}k+b_{2}\log^{2}k+O(\log^{3}k)\,.$

Expanding in $\log k$ gives the coefficients

(62) $\displaystyle b_{0}=0$ (63) $\displaystyle
b_{1}=-(1+\frac{3}{5}+2c_{1})=-\frac{3}{2}(\mu+1)$ (64) $\displaystyle
b_{1L}=0$ (65) $\displaystyle
b_{2}=\frac{293}{2,\\!100}-2c_{2}=\frac{61}{1,\\!050}+\frac{3}{40}(\mu+1)\,.$

We keep the terms proportional to $\mu+1$, although they vanish for $r=0$, in
order to keep track of their contributions for general $r$ in Sec. 3.3.

In the third line we used $c_{1}=-\frac{4}{5}+\frac{3}{4}(\mu+1)$ from (41) to
obtain the null result for $b_{1L}$. In the last line we used
$c_{2}=\frac{57}{1,\\!400}-\frac{3}{80}(\mu+1)$ from (42).

Substituting the expansion (61) into (60) we get the expansion in $\log k$

$\displaystyle\sigma^{2}_{\rm LN}(k,\tau)$ $\displaystyle=$
$\displaystyle\frac{\log^{2}k}{2J(k)}+\left\\{-\frac{16b_{1}}{9\log
k}-\frac{16}{45}\left(2b_{1}+5b_{2}\right)\right.$
$\displaystyle\left.-\frac{8(17b_{1}+420b_{2}+1,\\!050b_{3})}{4,\\!725}\log
k+O(\log^{2}k)\right\\}\tau+O(\tau^{2})\,.$

Note the presence of a singular term $1/\log k$ in the subleading volatility
proportional to $b_{1}$; since $b_{1}=-\frac{3}{2}(\mu+1)$ this divergent term
vanishes for the driftless gBM case $\mu=-1$. After a more careful analysis in
Sec. 3.3 it will be seen to cancel also for the gBM with non-zero drift.

Substituting the expressions for $b_{i}$ from (63), (65) into (60) gives an
explicit result for the ATM implied variance to $O(\tau)$ for the driftless
case $\mu=-1$

(67) $\displaystyle\sigma^{2}_{\rm
LN}(k,\tau)|_{\mu=-1}=\frac{x^{2}}{2J(e^{x})}+\left(-\frac{488}{4,\\!725}+O(x)\right)\tau+O(\tau^{2})\,.$

Rescaling to general $(\sigma,r,T)$ using (22) gives the first terms in the
equivalent log-normal volatility of an Asian option stated in Proposition 2.
Keeping only the ATM expression for the $O(T)$ term this is

(68) $\displaystyle\Sigma_{\rm
LN}^{2}\Big{(}\frac{K}{S_{0}},T\Big{)}|_{r=0}=\Sigma_{0}^{2}\Big{(}\frac{K}{A_{\rm
fwd}}\Big{)}+\sigma^{2}\left(-\frac{61}{9,\\!450}(\sigma^{2}T)+O(x)\right)+O(T^{2})\,.$

###  3.3. Including a non-zero interest rate $r$

In the last step of the proof we include the contribution of the rate $r$ into
the log-strike definition. We will show that this ensures the cancellation of
the divergent term proportional to $\mu+1$ in (3.2), and adds a new finite
term proportional to this factor.

Expanding the log-moneyness of the Asian option to order $O(\tau)$ we have

(69) $x=\log\frac{k}{a_{\rm fwd}^{(\mu)}}=\log k-(\mu+1)\tau+O(\tau^{2})$

where the forward price is the average of $A_{\tau}^{(\mu)}$ obtained using
$r=\mu+1$ in standardized units

(70) $a_{\rm
fwd}^{(\mu)}:=\mathbb{E}\left[\frac{1}{\tau}A_{\tau}^{(\mu)}\right]=\frac{1}{2(1+\mu)\tau}\left(e^{2(1+\mu)\tau}-1\right)\simeq
1+(\mu+1)\tau+O(\tau^{2})\,.$

###### Proof of the general result $r\neq 0$..

The small-$\tau$ expansion of $L=-\log c(k,\tau)$ at fixed $x$ is obtained by
replacing $k\to(1+(\mu+1)\tau)e^{x}$ into (57). To $O(\tau)$ it is sufficient
to replace $\log k\mapsto x+(\mu+1)\tau$. Expanding in $\tau$ we find

$\displaystyle L$ $\displaystyle=$
$\displaystyle\frac{1}{\tau}J(e^{x})+(\mu+1)e^{x}J^{\prime}(e^{x})-\frac{3}{2}\log\tau+\alpha_{0}(e^{x})+O(\tau)$
$\displaystyle=$
$\displaystyle\frac{1}{\tau}J(e^{x})+(\mu+1)\left(\frac{3}{4}x-\frac{9}{40}x^{2}+\frac{109}{1,\\!400}x^{3}+O(x^{4})\right)-\frac{3}{2}\log\tau+\alpha_{0}(e^{x})+O(\tau)$
$\displaystyle:=$
$\displaystyle\frac{1}{\tau}J(e^{x})-\frac{3}{2}\log\tau+\tilde{\alpha}_{0}(e^{x})+O(\tau)\,.$

In the last step we absorbed the second term in the second line into
$\tilde{\alpha}_{0}(e^{x})$. The effect of the new term is to shift the
coefficients $b_{k}$ defined in (61) as

(72) $\displaystyle b_{1}\to\tilde{b}_{1}:=b_{1}+\frac{3}{2}(\mu+1)=0$ (73)
$\displaystyle
b_{2}\to\tilde{b}_{2}:=b_{2}-\frac{9}{20}(\mu+1)=\frac{61}{1,050}-\frac{3}{8}(\mu+1)\,.$

Substituting $b_{i}\to\tilde{b}_{i}$ into (3.2) gives

$\displaystyle\sigma^{2}_{LN}(e^{x},\tau)$ $\displaystyle=$
$\displaystyle\frac{x^{2}}{2J(e^{x})}+\left\\{-\frac{16\tilde{b}_{1}}{9x}-\frac{16}{45}\left(2\tilde{b}_{1}+5\tilde{b}_{2}\right)\right.$
$\displaystyle\left.-\frac{8(17\tilde{b}_{1}+420\tilde{b}_{2}+1,\\!050\tilde{b}_{3})}{4,\\!725}x+O(x^{2})\right\\}\tau+O(\tau^{2})$

Substituting here the explicit results for $\tilde{b}_{i}$ from (72), (73) we
get

(75) $\sigma^{2}_{\rm
LN}(e^{x},\tau)=\frac{x^{2}}{2J(e^{x})}+\left(-\frac{488}{4,\\!725}+\frac{2}{3}(\mu+1)-\frac{544}{23,\\!625}x+O(x^{2})\right)\tau+O(\tau^{2})\,.$

All singular terms as $x\to 0$ cancel out, and the Asian volatility
$\sigma_{\rm LN}(e^{x},\tau)$ is finite and well-defined at the ATM point
$x=0$.

Rescaling to arbitrary volatility $\sigma$ and the actual maturity $T$ using
$\tau\to\frac{1}{4}\sigma^{2}T,\mu+1\to\frac{2r}{\sigma^{2}}$, see (22), gives
the final result

$\displaystyle\Sigma^{2}_{\rm LN}(K,T,S_{0})$ $\displaystyle=$
$\displaystyle\sigma^{2}\Big{\\{}\frac{\log^{2}x}{2J_{\rm
BS}(x)}-\frac{61}{9,\\!450}(\sigma^{2}T)+\frac{1}{12}(rT)$
$\displaystyle\hskip
14.22636pt-\frac{34}{23,\\!625}(\sigma^{2}T)x+O(x^{2}T)+O(T^{2})\Big{\\}}\,.$

This concludes the proof of Proposition 2. ∎

Numerically the ATM implied variance is
$\frac{1}{3}-0.00645(\sigma^{2}T)+0.083(rT)$, such that the $O(rT)$ term
dominates the subleading contribution for most realistic values of the model
parameters. The correction linear in log-moneyness is
$-0.00144(\sigma^{2}T)x$.

###  3.4. Consistency check and an improved estimate

The small maturity limit of the equivalent log-normal volatility of an Asian
option in the Black-Scholes model has been obtained in [17] in a modified
small maturity limit $\sigma^{2}T\to 0$ taken at fixed $\rho=rT$. This limit
takes into account interest rates effects; more precisely it includes
corrections of the order $O((rT)^{n})$ to all orders in $n$.

The result is given in Proposition 19 of [17] and we denote it as $\Sigma_{\rm
LN,\rho}(K,\rho)$. As $\rho\to 0$, this reduces to $\Sigma_{0}(K/S_{0})$ given
in (10).

The Asian volatility $\Sigma_{\rm LN,\rho}(K,\rho)$ includes corrections of
the order $O((rT)^{n})$ to all orders in $n$. At the ATM point this function
simplifies and is given by (see equation (125) in [17])

(77) $\Sigma_{LN,\rho}(K=A_{\rm fwd},\rho)=\sigma\frac{S_{0}}{A_{\rm
fwd}}\sqrt{v(\rho)}=\sigma\frac{\rho}{e^{\rho}-1}\sqrt{v(\rho)}$

with $A_{\rm fwd}=S_{0}\frac{e^{\rho}-1}{\rho}$ and

(78) $v(\rho):=\frac{1}{\rho^{3}}(\rho
e^{2\rho}-\frac{3}{2}e^{2\rho}+2e^{\rho}-\frac{1}{2})=\frac{1}{3}+\frac{5}{12}\rho+\frac{17}{60}\rho^{3}+O(\rho^{4})$

We will use this result to test the coefficient of the $O(rT)$ term in (3.3).

Squaring (77) and expanding in $\rho$ gives

(79) $\Sigma_{LN,\rho}^{2}(K=A_{\rm
fwd},\rho)=\sigma^{2}\frac{\rho^{2}}{(e^{\rho}-1)^{2}}v(\rho)=\sigma^{2}\left(\frac{1}{3}+\frac{1}{12}\rho+\frac{1}{180}\rho^{2}+O(\rho^{3})\right)\,.$

This reproduces the $+\frac{1}{12}(rT)$ correction in Eq. (3.3).

This suggests an improved approximation for $\Sigma_{LN}(K,T)$, obtained by
replacing $\sigma^{2}\frac{x^{2}}{2J_{BS}(e^{x})}$ in (3.3) with
$\Sigma^{2}_{LN,\rho}(K,\rho)$. This approximation includes interest rates
effects through the leading order term, by taking into account corrections of
order $O((rT)^{n})$ to all orders, in addition to the $O(\sigma^{2}T)$
correction computed here. This approximation is somewhat heuristic, as it
neglects e.g. corrections of $O((\sigma^{2}T)^{n})$ with $n>1$. It should be
useful in situations when the $O(rT)$ corrections dominate numerically over
$O(\sigma^{2}T)$.

We denote this improved next-to-leading order (NLO) estimate as
$\Sigma_{LN,NLO}(K,T)$. It is given explicitly by

(80) $\Sigma^{2}_{\rm
LN,NLO}(K,T):=\Sigma_{LN,\rho}^{2}(K,\rho)+\sigma^{2}\Big{(}-\frac{61}{9,\\!450}(\sigma^{2}T)-\frac{34}{23,\\!625}(\sigma^{2}T)\log\frac{K}{A_{\rm
fwd}}\Big{)}\,.$

This approximation can be further improved by adding terms of higher order in
log-moneyness $x=\log(K/A_{\rm fwd})$. Keeping terms up to the linear term in
$x$ should give an accurate approximation in a region of strikes sufficiently
close to the ATM point.

##  4\. Numerical examples

In this section we present a few numerical tests of our results. We can price
Asian options by substituting the equivalent log-normal volatility
$\Sigma_{\rm LN}(K,T)$ of Proposition 2 into the Black-Scholes formula. This
gives the Asian prices

(81) $C(K,T)=e^{-rT}[A_{\rm fwd}\Phi(d_{1})-K\Phi(d_{2})]\,,\quad
P(K,T)=e^{-rT}[K\Phi(-d_{2})-A_{\rm fwd}\Phi(-d_{1})]\,,$

with $A_{\rm fwd}$ given in (9) and

(82) $d_{1,2}=\frac{1}{\Sigma_{\rm LN}(K,T)\sqrt{T}}\left(\log\frac{A_{\rm
fwd}}{K}\pm\frac{1}{2}\Sigma_{\rm LN}^{2}(K,T)T\right)\,.$

Using this approach we evaluate the seven benchmark cases given in Linetsky
[12], and compare them against the precise results obtained in this paper
using a spectral expansion approach. Table 1 shows the results for Asian
option prices obtained from the leading order asymptotic result of [17]
(column $C_{0}(K,T)$) and the improved results obtained keeping also the
$O(T)$ subleading correction. The columns $C_{1}^{ATM}(K,T)$ and $C_{1}^{\rm
lin}(K,T)$ show the results obtained by keeping only the ATM subleading
correction $O(T)$ in (2), and including also the term linear in log-strike
$O(Tx)$, respectively. The last column shows the benchmark results from [12]
obtained using a precise spectral expansion. In brackets we show the relative
error with respect to the benchmark, for each approximation.

The results are shown also in Figure 1. The plots show the equivalent log-
normal volatility $\Sigma_{LN}(K,T)$ including only the ATM subleading
correction (first line of (2)) vs $k=K/S_{0}$ (solid curves), comparing with
the leading order result $\Sigma_{0}(k)$ (dashed curves). The result for
$\Sigma_{\rm LN}(K,T)$ (2) depends only on $(\sigma,r,T)$, so cases 4,5,6 have
a common curve. The dots show the precise benchmark values in the last column
of Table 1 converted to volatility. The vertical lines show the ATM strike
$A_{\rm fwd}(S_{0})/S_{0}$.

The agreement with the benchmark results improves significantly when including
the subleading correction, especially in cases with large $rT$. The error of
$C_{1}^{ATM}(K,T)$ is below 0.02% in all cases, and becomes even smaller for
$C_{1}^{\rm lin}(K,T)$ when including the subleading skew contribution.

The improved approximation (80) which includes corrections of order
$O((rT)^{n})$ to all orders is tested in Table 2 against the same benchmark
cases from [12]. This approximation is expected to perform better for cases
with large $rT$. This is confirmed indeed, as seen for case 3 which has the
largest values of this parameter $rT=0.18$ and agrees with the benchmark to
five digits. The approximation error is below 0.02% in all cases.

Table 1.  Seven benchmark cases for Asian options. $C_{0}(K,T)$ denotes the Asian options obtained using the leading approximation $\Sigma_{0}(k)$ for the equivalent log-normal Asian volatility. $C_{1}^{\rm lin}(K,T)$ denotes the Asian option price obtained by including both terms in the subleading correction (2), and $C_{1}^{ATM}(K,T)$ includes only the ATM subleading correction. The last column shows the precise results of [12] obtained by a spectral expansion. Relative errors relative to the benchmarks are shown in brackets [bps]. Case | $k$ | $r$ | $\sigma$ | $T$ | $C_{0}(K,T)$ | $C_{1}^{ATM}(K,T)$ | $C_{1}^{lin}(K,T)$ | benchmark
---|---|---|---|---|---|---|---|---
1 | 1 | 0.02 | 0.10 | 1 | 0.055923 (-11.2) | 0.055986 (0.0) | 0.055986 (0.0) | 0.055986
2 | 1 | 0.18 | 0.30 | 1 | 0.217054 (-61.0) | 0.218362 (-1.1) | 0.218364 (-1.0) | 0.218387
3 | 1 | 0.0125 | 0.25 | 2 | 0.172163 (-6.2) | 0.172268 (-0.1) | 0.172269 (0.0) | 0.172269
4 | $\frac{2}{1.9}$ | 0.05 | 0.50 | 1 | 0.192895 (-14.4) | 0.193176 (0.1) | 0.193173 (0.0) | 0.193174
5 | 1 | 0.05 | 0.50 | 1 | 0.246125 (-11.8) | 0.246412 (-0.2) | 0.246415 (0.0) | 0.246416
6 | $\frac{2}{2.1}$ | 0.05 | 0.50 | 1 | 0.305927 (-9.6) | 0.306211 (-0.3) | 0.306220 (0.0) | 0.306220
7 | 1 | 0.05 | 0.50 | 2 | 0.349314 (-22.3) | 0.350077 (-0.5) | 0.350087 (-0.2) | 0.350095
Table 2.  The predictions for Asian options prices obtained using the improved approximation for the equivalent log-normal volatility $\Sigma_{LN,NLO}(K,T)$ in (80) including terms of all orders in $O((rT)^{n})$ (NLO), comparing with the benchmarks of Linetsky [12]. The scenarios are the same as in Table 1. Last row shows the relative error in basis points. Case | 1 | 2 | 3 | 4 | 5 | 6 | 7
---|---|---|---|---|---|---|---
NLO | 0.055986 | 0.218385 | 0.172268 | 0.193188 | 0.246409 | 0.306193 | 0.350060
Linetsky | 0.055986 | 0.218387 | 0.172269 | 0.193174 | 0.246416 | 0.306220 | 0.350095
err [bp] | 0 | 0.09 | 0.05 | 0.64 | -0.32 | -1.24 | -1.60

Figure 1.  Asian volatilities $\Sigma_{\rm LN}(K,T)[\%]$ vs $k=K/S_{0}$ for
the seven scenarios in Table 1. The dashed curves show the leading order Asian
volatility $\Sigma_{0}(K/S_{0})$ and the solid curves include the ATM
subleading correction in (2). The dots show the benchmark cases in the last
column of Table 1. The vertical line shows the ATM strike $A_{\rm
fwd}(S_{0})/S_{0}$.

Note added (July 2024). We summarize here the result for the Asian implied
variance, in a form which makes explicit the dependence on the coefficients
$c_{i}$ and is easier to extend to higher orders. This is used to extend
Proposition 2 by including also the convexity term of the subleading Asian
implied variance.

The reduced Asian implied variance to order $O(\tau)$ is

(83) $\displaystyle\sigma_{LN}^{2}(k,\tau)=\frac{\log^{2}k}{2J(k)}$
$\displaystyle+\Big{\\{}\frac{4}{4,\\!725}(1,\\!051+1,\\!680c_{1}+4,\\!200c_{2}-315(\mu+1))$
$\displaystyle\quad+\frac{8}{23,\\!625}(-91+170c_{1}+4,\\!200c_{2}+10,\\!500c_{3})\log
k$
$\displaystyle+\frac{1}{18,\\!191,\\!250}(-250,\\!193-517,\\!440c_{1}+1,\\!047,\\!200c_{2}+25,\\!872,\\!000c_{3}$
$\displaystyle\hskip
71.13188pt+64,\\!680,\\!000c_{4}-39,\\!270(\mu+1))\log^{2}k$
$\displaystyle\quad+O(\log^{3}k)\Big{\\}}\tau+O(\tau^{2})$

The first term is the leading implied variance, to all orders in $\log k$. The
second, third and fourth terms give the ATM level, skew and convexity of the
subleading implied variance, respectively.

Substituting here the coefficients $c_{1-3}$ given in (41) - (43) and

(84) $\displaystyle
c_{4}=\frac{11}{22,\\!400}(\mu+1)+\frac{3,\\!281}{6,\\!160,\\!000}$

gives the following result for the subleading Asian implied variance.

(85)
$\displaystyle\sigma_{LN}^{2}(e^{x},\tau)=\frac{x^{2}}{2J(e^{x})}+\Big{(}-\frac{488}{4,\\!725}+\frac{2}{3}(\mu+1)-\frac{544}{23,\\!625}x$
$\displaystyle\hskip
56.9055pt+\,\Big{(}\frac{12,\\!073}{1,\\!034,\\!500}-\frac{5}{525}(\mu+1)\Big{)}x^{2}+O(x^{3})\Big{)}\tau+O(\tau^{2})$

where $x=\log(K/A_{\rm fwd})$ is the option log-moneyness.

Rescaling to arbitrary Black-Scholes parameters gives the analog of equation
(12), including also the convexity term.

$\displaystyle\Sigma_{\rm LN}^{2}(K,T)$ $\displaystyle=$
$\displaystyle\sigma^{2}\Big{\\{}\underbrace{\frac{\log^{2}x}{2J_{\rm
BS}(x)}}\underbrace{-\frac{61}{9,\\!450}(\sigma^{2}T)+\frac{1}{12}(rT)}$
$\displaystyle\hskip 28.45274ptO(1)\hskip 42.67912ptO(T)$ $\displaystyle\hskip
14.22636pt+\underbrace{\Big{[}-\frac{34}{23,\\!625}(\sigma^{2}T)\Big{]}x}$
$\displaystyle\hskip 71.13188ptO(Tx)$ $\displaystyle\hskip
14.22636pt+\underbrace{\Big{[}\frac{12,\\!073}{16,\\!632,\\!000}(\sigma^{2}T)-\frac{5}{2,\\!016}(rT)\Big{]}x^{2}}+O(Tx^{3})+O(T^{2})\Big{\\}}\,.$
$\displaystyle\hskip 71.13188ptO(Tx^{2})$

##  Appendix A Asymptotic expansion for integrals

The following theorem is due to Erdélyi and is given in Sec. 2.4 of [6] (p.
36). It appears as Theorem 8.1 in Chapter 3.8 of Olver [15]. For convenience
we quote it below in the notations of Theorem 1.2.1 in Nemes [14].

######  Theorem 7.

Consider the integral

(87) $I(\lambda)=\int_{a}^{b}e^{-\lambda f(x)}g(x)dx$

Assume that:

(i) $f(x)>f(a)$ for all $x\in(a,b)$.

(ii) $f^{\prime}(x),g(x)$ are continuous in a neighborhood of $a$.

(iii) the following expansions hold

(88) $\displaystyle f(x)=f(a)+\sum_{k=0}^{\infty}a_{k}(x-a)^{\alpha+k}$
$\displaystyle g(x)=g(a)+\sum_{k=0}^{\infty}b_{k}(x-a)^{\beta+k-1}\,.$

(iv) $I(\lambda)$ converges absolutely for all sufficiently large $\lambda$.

Then

(89) $I(\lambda)=e^{-\lambda
f(a)}\sum_{n=0}^{\infty}\Gamma\Big{(}\frac{n+\beta}{\alpha}\Big{)}\frac{d_{n}}{\lambda^{(n+\beta)/\alpha}}\,.$

The coefficient of the leading order term is

(90) $d_{0}=\frac{b_{0}}{\alpha a_{0}^{\beta/\alpha}}\,.$

## References

  * [1] E. Alòs, E. Nualart and M. Pravosud, On the implied volatility of Asian options under stochastic volatility models, arxiv:2208.01353[q-fin.MF]
  * [2] L.P. Arguin, N.-L. Liu and T.-H. Wang, Most-likely-path in Asian option pricing under local volatility models. International Journal of Theoretical and Applied Finance. 21, 1850029 (2018)
  * [3] D. Dufresne, Laguerre series for Asian and other options, Math. Finance 10, 407-428 (2000)
  * [4] D. Dufresne, The integral of the geometric Brownian motion, Adv. Appl. Probab. 33, 223-241 (2001).
  * [5] D. Dufresne, Bessel processes and a functional of Brownian motion, in M. Michele and H. Ben-Ameur (Eds.), Numerical Methods in Finance 35-57, Springer, 2005.
  * [6] A. Erdélyi, Asymptotic expansions, Dover Publications 1956
  * [7] P. Foschi, Pagliarani, S. and A. Pascucci. Approximations for Asian options in local volatility models. Journal of Computational and Applied Mathematics. 237, 442-459 (2013)
  * [8] G. Fusai and A. Roncoroni, Implementing Models in Quantitative Finance: Methods and Cases, Springer Finance Series, Springer, New York, 2008
  * [9] H. Geman and M. Yor, Bessel processes, Asian options and perpetuities, Mathematical Finance 3, 349-375 (1993)
  * [10] K. Gao and R. Lee, Asymptotics of implied volatility to arbitrary order, Finance and Stochastics 18(2), 349-392 (2014)
  * [11] A. Gulisashvili and E.M. Stein, Asymptotic behavior of distribution densities in models with stochastic volatility.I. Mathematical Finance 20(3), 447-477 (2010)
  * [12] V. Linetsky, Spectral expansions for Asian (Average price) options, Operations Research 52, 856-867 (2004)
  * [13] P. Nándori and D. Pirjol, On the distribution of the time-integral of the geometric Brownian motion, J. Comp. Appl. Math. 402, 113818 (2022)
  * [14] G. Nemes, Asymptotic expansion for integrals, MSc Thesis, L. Ötvös University 2012.
  * [15] F.W.J. Olver, Introduction to asymptotics and special functions, Academic Press 1974.
  * [16] D. Pirjol, L. Zhu, Short maturity Asian options in local volatility models, SIAM J. Finan. Math. 7(1), 947-992 (2016); arXiv:1609.07559[q-fin].
  * [17] D. Pirjol, L. Zhu, Asymptotics for the discrete-time average of the geometric Brownian motion, Adv. Appl. Prob. 49, 1-53 (2017)
  * [18] D. Pirjol, L. Zhu, Sensitivities of Asian options in the Black-Scholes model, International Journal of Theoretical and Applied Finance 21(01), 1850008 (2018)
  * [19] D. Pirjol and L. Zhu. Short maturity Asian options for the CEV model, Probab. Eng. Inform. Sc. 33(2), 258-290, 2019.
  * [20] D. Pirjol and L. Zhu. Asymptotics for short maturity Asian options in jump-diffusion models with local volatility, Quantitative Finance 24(3-4), 433-449, 2024.
  * [21] D. Pirjol, Small-$t$ asymptotic expansion of the Hartman-Watson distribution, Method. Comp. Appl. Prob. 23(4), 1537-1549 (2021); expanded version available as arxiv:2001.09579 [math.PR]
  * [22] L. Rogers and Z. Shi, The value of an Asian option. J. Appl. Prob. 32, 1077-1088 (1995).
  * [23] M. Yor, On some exponential functionals of Brownian motion, Adv. Appl. Probab. 24, 509-531 (1992)
  * [24] L. Zhu, Options with extreme strikes, Risks 3(3), 234-249 (2015)


# Ebola transmission dynamics: will future Ebola outbreaks become cyclic?

David Niyukuri<EMAIL_ADDRESS>Kelly Joëlle Gatore Sinigirira  Jean
De Dieu Kwizera  Salma Omar Abd-Almageid Adam  Department of Mathematics,
University of Burundi, Burundi  Doctoral School, University of Burundi,
Burundi  The South African Department of Science and Technology–National
Research Foundation (DST-NRF) Centre of Excellence in Epidemiological
Modelling and Analysis (SACEMA), Stellenbosch University, Cape Town, South
Africa  Institut Superieur des Cadres Militaires, Bujumbura, Burundi
Department of Mathematics, College of Computer Sciences and Mathematics,
University of Bahri, Karthoum, Sudan

###### Abstract

The Ebola Virus Disease (EVD) can persist in some body fluids after clinical
recovery. In Guinea and the Democratic Republic of the Congo, there are well-
documented cases of EVD re-emergence that were associated with previous
outbreaks. In many cases, male EVD survivors were associated with the re-
introduction of new outbreaks. This has shown that even after controlling an
EVD outbreak, a close biomedical monitoring of survivors and contacts is
critical to avoiding future outbreaks. Thus, in order to explore the main
features of EVD transmission dynamics in the context of re-emergence, we used
a compartmental model by considering vaccination around EVD contacts.
Analytical and numerical analyses of the model were conducted. The model is
mathematically and epidemiologically well-posed. We computed the reproductive
number ($R_{0}$) and the disease equilibrium points (disease-free equilibrium
and endemic equilibrium) for the re-emerging outbreak. The stability analysis
of the model around those equilibrium points was performed. The model
undergoes a backward bifurcation at $R_{0}$ close to 1, regardless $R_{0}<1$,
the disease will not be eradicated. This means that R0 cannot be considered as
an intervention control measure in our model.

###### keywords:

Ebola , transmission dynamic , resurgence

††journal: Epidemics

##  1 Introduction

It is known that the Ebola Virus Disease (EVD) infection can persist in some
body fluids after clinical recovery and clearance of virus from the
bloodstream. In a systematic review, authors highlighted that EVD viral RNA
can still be detected in semen (day 272), aqueous humor (day 63), sweat (day
40), urine (day 30), vaginal secretions (day 33), conjunctival fluid (day 22),
faeces (day 19) and breast milk (day 17) [1].Thus, this has important
implication on EVD transmission dynamics and control strategies. In the past,
sexual transmission from recovered Ebola survivors has been documented [2, 3].

During the EVD outbreak in West Africa in 2014-2015, the World Health
Organization (WHO) declared the end of Ebola virus (EBOV) transmission in the
Republic of Guinea on 29 December 2015, a resurgence of EVD cases was reported
in February 2016 almost after 2 months [4, 5]. Surprisingly, analysis showed
that this resurgence in Guinea was linked to previous outbreak of 500 days
ago. In a longitudinal study in Liberia of a cohort of EVD survivors and their
close contacts to understand EVD sequelae showed that the EVD viral RNA in
semen persisted for a maximum of 40 months [6]. Recent EVD outbreaks in the
Democratic Republic of Congo (DRC) and Guinea in 2021, were close to previous
outbreaks in those respective countries.

In the Democratic Republic of the Congo (DRC), on 7 February 2021, an Ebola
virus disease (EVD) outbreak was declared by the Ministry of Health (MOH) of
the DRC. On 8 October 2021, the Ministry of Health announced an outbreak of
Ebola virus disease (EVD) in Beni Health Zone, North Kivu Province. According
to phylogenetic analysis, these sporadic outbreaks were linked to the
2018-2020 Nord Kivu/Ituri EVD outbreak, hence they were not new spillover
events [7].

In the Republic of Guinea, on 14 February 2021, Guinean health authorities
declared an EVD outbreak, the lineage shows considerably lower divergence than
would be expected during sustained human-to-human transmission, which suggests
a persistent infection with reduced replication or a period of latency [8]. A
phylogenetic tree shows the new virus falls between virus samples from the
2013–16 epidemic [9]. Thus, the virus which caused recent outbreaks barely
differs from the strain seen 5 to 6 years ago, and this becomes a puzzle which
led to different hypothesis. The biology of Ebola virus does not support if
the virus could lay dormant for 5 years in a survivor of the epidemic all that
time. Dan Bausch from the United Kingdom’s Public Health Rapid Support Team
argued that “From the tree, you’d conclude that it is a virus that persisted
in some way in the area, and sure, most likely in a survivor”. He added that
it is also hard to rule out scenarios such as a small, unrecognized chain of
human-to-human transmission. “For example, a 2014 survivor infects his wife a
few years after recovery, who infects another male, who survives and carries
virus for a few years, then infecting another women, who is then seen by a
nurse who dies”—the index case in the new outbreak, he said [10].

Figure 1: Re-introduction of EVD in Guinea and the Democratic Republic of
Congo. In Guinea, the Current EVD outbreak was identified to be close to Ebola
sequence samples in infected individuals in 2014. Same scenario occurred in
February 2016, a re-introduction of EVD from a man survivor was reported in
Guinea. In the Democratic Republic of Congo, the recent EVD outbreak was
associated to the first case which was a woman married to an EVD survivor of
2019.

From those reported cases, we can see that EVD epidemics can exhibit
transmission cycles in which survivors play a pivotal role for outbreak re-
emergence. From above mentioned reported cases, sexual transmission was the
main transmission mode for new EVD outbreak re-introduction. Thus, EVD is
becoming an emerging sexually transmissible infection of great concerns [11].
Thus, the overall aim of this modelling study is to get insights on which
conditions EVD outbreaks will re-emerge from previous outbreaks. Given the
effectiveness of the vaccine for Zaire strains which caused the above
mentioned outbreaks, we will also determine the value of the reproductive
number (R0), the disease-free equilibrium (DFE), and endemic equilibrium (EE)
[12, 13], for the re-emerging outbreak.

##  2 Model description and formulation

To study the transmission cycles of EVD dynamics, we suggested a Susceptible-
Vaccinated-Exposed-Infected-Recovery (SVEIR) type model with 5 compartments
which are described in Figure (2), consisting of coupled nonlinear
differential equations [13, 14].In our model formulation, we assume that the
very first infection of Ebola virus might be a spillover from wild animals to
human subjects. We also assume that recovered individuals can trig new
infections as it has been reported in DRC and Guinea. We perform analytical
analysis to look at the existence of solutions, and stability analysis of the
equilibria of the model.

Figure 2: Schematic representation of mathematical model of EVD transmission.

Explicitly, we have the class of Susceptible (S) humans at time instant t,
this is a class of individuals who are not yet infected but can get infected
by the Ebola virus; Given the availability of EVD vaccination, it is feasible
to consider a specific vaccinated class (V) which is a class of individuals
who are vaccinated. Exposed (E) individuals, they are individuals who were
susceptible and vaccinated and have become exposed after being in contact with
any of Ebola infected person of contaminated body fluids or fomites (infected
inanimate objects and surfaces) [1, 15], these people are placed in this class
which consists of exposed humans at time instant t; Infected (I) humans at
time instant t, exposed who develop the disease; and Recovered (R) class is
for infected individuals who survive and finally recover or acquire immunity
due to vaccination.

The mathematical model of EVD transmission dynamic is given by the
instantaneous change of the flow between compartment, which is described by
the following system of non-linear differential equations 10 which assumes
$\lambda_{s}=\frac{1}{N}(b_{i}\beta_{i}I+b_{r}\beta_{r}R)$ and
$\lambda_{v}=\frac{1}{N}(c_{i}\alpha_{i}I+c_{r}\alpha_{r}R)$.

$\displaystyle\left\\{\begin{array}[]{rcl}\frac{dS}{dt}&=&p-\frac{S}{N}(b_{i}\beta_{i}I+b_{r}\beta_{r}R)-(\mu+\varphi)S\\\
\\\ \frac{dV}{dt}&=&\varphi
S-\frac{V}{N}(c_{i}\alpha_{i}I+c_{r}\alpha_{r}R)-\mu V\\\ \\\
\frac{dE}{dt}&=&\frac{S}{N}(b_{i}\beta_{i}I+b_{r}\beta_{r}R)+\frac{V}{N}(c_{i}\alpha_{i}I+c_{r}\alpha_{r}R)-(\mu+\omega)E\\\
\\\ \frac{dI}{dt}&=&\omega E-(\gamma+\mu+\delta)I\\\ \\\
\frac{dR}{dt}&=&\gamma I-\mu R\end{array}\right.$ (10)

with $N(t)=S(t)+V(t)+E(t)+I(t)+R(t)$, which means that the population change
over time.

Table 1: Parameters for the S-V-E-I-R Ebola Virus Disease dynamic model Parameter |  Definition
---|---
$p$ |  Recruitment rate of individuals into the population
$\varphi$ |  Proportions of susceptible individuals who are vaccinated
$\mu$ |  Natural death rate
$b_{i}$ |  contact rate between infected individuals and susceptible individuals
$\beta_{i}$ |  probability for susceptible individuals of being in contact with infected individuals that
|  occurs an infection.
$b_{r}$ |  contact rate between recovered individuals and susceptible individuals
$\beta_{r}$ |  probability for susceptible individuals of being in contact with recovered individuals that
|  occurs an infection.
$c_{i}$ |  contact rate between infected individuals and vaccinated individuals
$\alpha_{i}$ |  probability for vaccinated individuals of being in contact with infected individuals that
|  occurs an infection.
$c_{r}$ |  contact rate between recovered individuals and vaccinated individuals
$\alpha_{r}$ |  probability for vaccinated individuals of being in contact with recovered individuals that
|  occurs an infection.
$\omega$ |  Exit rate from exposed class.
$\gamma$ |  Exit rate from infected class.
$\delta$ |  disease induced death rate

New EVD infections are results of the interaction between susceptible or
vaccinated and infected individuals respectively at rate $\beta_{i}$ and
$\alpha_{i}$, this gives the rate at which individuals are being exposed to
the infection. We also assume that individuals who are considered recovered
can interact with susceptible and vaccinated individuals respectively at rate
$\beta_{r}$ and $\alpha_{r}$ hence susceptible individuals become exposed by
interacting with EVD survivors, and this can cause re-introduction of EVD.
This is supported by current evidence of resurgence or re-introduction of EVD
through sexual intercourse. The rate of increase or decrease of susceptible
population is given by $p$, which can be by natural deaths, migrations or new
births. Individuals die at rate $\mu$ by natural deaths. EVD death rate is
given by $\delta$. The incubation period also known as latent period which is
the time between infection and sickness is given by $\omega$. All parameters
of the SVEIR model are given in Table 1.

We have seen that new introduction of EVD infections have been driven by men
EVD survivors. Thus, the re-introduction parameter $\alpha_{r}$ and
$\beta_{r}$ should depend also on proportion of sexually active men who
survived the infection. To solve the system of nonlinear differential
equations, we used the numerical method Runge–Kutta of order 4 (RK4). To fit
the model, and estimate some parameters of the model, we used record data of
reported EVD cases in the Republic of Guinea and DRC.

##  3 Model analysis

The analytical analysis of the systems of ODEs (10) shows the dynamical
behaviour of the solutions of the system. With the intrinsic nature of
transmission dynamic through contact, we will estimate the basic reproduction
number R0, analyze different equilibria of the model [12, 13, 16, 17], at the
disease-free equilibrium (DFE), and endemic equilibrium (EE) points. We will
also analyze the attainment of herd immunity.

###  3.1 Bounded and invariant region

The positivity and boundness of the solution are investigated to show that the
model is biologically and mathematically well posed.

Theorem: The feasible region $\Omega$ defined by:

$\displaystyle\Omega$ $\displaystyle=$
$\displaystyle\\{S(t);V(t),E(t),I(t),R(t)\in\mathbb{R}_{+}^{5}:N(t)\leq\dfrac{p}{\mu}\\}$
(11)

with $S(0),V(0),E(0),I(0),R(0)\leq 0$ is positively invariant.

###### Proof.

Adding the equation of the system 10, We get:

$\displaystyle\dfrac{dN}{dt}$ $\displaystyle\leq$ $\displaystyle p-\mu N$ (12)

Solving 12, We have:

$\displaystyle 0\leq N(t)\leq\dfrac{p}{\mu}+N(0)e^{-\mu t}$

$N(0)$ is the initial values of the total population $N(t)$ and

$\displaystyle\lim\limits_{t\rightarrow+\infty}supN(t)\leq\dfrac{p}{\mu}$ (13)

So $\Omega=\\{S(t);V(t),E(t),I(t),R(t)\\}$ is positively invariant, and the
right handside of equations of the model (10) is continuous with continuous
partial derivative in $\Omega$. Hence, there is a unique solution. ∎

###  3.2 Disease-free equilibrium

The disease-free equilibrium (DFE) is the state at which there is no
infectious in the population.
So, $E=I=R=0$ and
$\frac{dS}{dt}=\frac{dV}{dt}=\frac{dE}{dt}=\frac{dI}{dt}=\frac{dR}{dt}=0$
Then, the model admit the disease free equilibrium point
$E_{dfe}=(S_{0},V_{0},0,0,0)$ with

$\displaystyle S_{0}=\frac{p}{\mu+\varphi}\mbox{ and
}V_{0}=\frac{p\varphi}{\mu(\mu+\varphi)}$ (14)

###  3.3 Basic reproduction number

We now determine an important threshold parameter in epidemic modelling known
as the basic reproduction number $(R_{0})$, it gives the average number of
generated infection cases by an infectious individual. Using the New
Generation Number described in (references), we will compute the basic
reproduction number $R_{0}$. Let’s define the transmission matrix $f$ and $v$

$\displaystyle f=\begin{pmatrix}\lambda_{s}S+\lambda_{v}\\\ 0\\\
0\end{pmatrix}\mbox{ and }v=\begin{pmatrix}(\mu+\omega)E\\\ -\omega
E+(\gamma+\mu+\delta)I\\\ -\gamma I+\mu R\end{pmatrix}$ (15)

The corresponding Jacobian matrix $F$ and $V$ of $f$ and $v$ evaluated at the
DFE are respectively:

$\displaystyle
F_{E_{dfe}}=\begin{pmatrix}0&\dfrac{\mu}{\mu+\varphi}(b_{i}\beta_{i}+\dfrac{c_{i}\alpha_{i}\varphi}{\mu})&\dfrac{\mu}{\mu+\varphi}(b_{r}\beta_{r}+\dfrac{c_{r}\alpha_{r}\varphi}{\mu})\\\
0&0&0\\\ 0&0&0\end{pmatrix}\mbox{ and }V=\begin{pmatrix}\mu+\omega&&0\\\
-\omega&(\gamma+\mu+\delta)&0\\\ 0&-\gamma&\mu\end{pmatrix}$ (16)

$FV^{-}{1}$ is the next generation number given by

$\displaystyle\left(\begin{array}[]{c}\dfrac{\omega(\mu b_{i}\beta_{i}+\varphi
c_{i}\alpha_{i})}{(\mu+\varphi)(\mu+\omega)(\mu+\gamma+\delta)}+\dfrac{\omega\gamma(\mu
b_{r}\beta_{r}+\varphi
c_{r}\alpha_{r})}{\mu(\mu+\varphi)(\mu+\omega)(\mu+\gamma+\delta)}\\\ 0\\\
0\end{array}\right.$ $\displaystyle\left.\begin{array}[]{ccc}\dfrac{(\mu
b_{i}\beta_{i}+\phi\alpha_{i}c_{i})}{(\mu+\varphi)(\mu+\gamma+\delta)}+\dfrac{(\mu
b_{r}\beta_{r}+\varphi
c_{r}\alpha_{r})}{\mu(\mu+\varphi)(\mu+\omega)(\mu+\gamma+\delta)}&\dfrac{1}{\mu(\mu+\varphi)}(\mu
b_{r}\beta_{r}+\varphi c_{r}\alpha_{r})\\\ 0&0\\\ 0&0\end{array}\right)$ (21)

It follows that the $R_{0}$ is given by $R_{0}=\rho(FV^{-1})$ Where $\rho$ is
the spectral radius of the new generation infection
So, the basic reproduction number becomes

$\displaystyle R_{0}=\rho FV^{-1}=\dfrac{\omega[\mu(\mu
b_{i}\beta_{i}+c_{i}\alpha_{i}\varphi)+\gamma(\mu
b_{r}\beta_{r}+c_{r}\alpha_{r}\varphi)]}{\mu(\mu+\varphi)(\mu+\omega)(\mu+\gamma+\delta)}$
(22)

###  3.4 Local stability of DFE

Theorem: The disease free equilibrium point, $E_{dfe}$ is locally
asymptotically stable if $R_{0}<1$ and unstable if $R_{0}>1$.

###### Proof.

Using the Jacobian matrix of Equations[Ref], we established the local
stability of the disease-free equilibrium and evaluated it at disease free
equilibrium $E_{dfe}$.This is accomplished by evaluating Equations [Ref]’s
Jacobian matrix at
$E_{dfe}=(\dfrac{p}{\mu+\varphi};\dfrac{p\varphi}{\mu(\mu+\varphi)};0;0,0)$
The eigenvalues of the Jacobian matrix of the model equations at $E_{dfe}$ are
used to determine the model’s local stability. As a result, the Jacobian
matrix of Equations [Ref] is given by

$\displaystyle J_{E_{dfe}}=\begin{pmatrix}-(\mu+\varphi)&0&0&-\dfrac{\mu
b_{i}\beta_{i}}{\mu+\varphi}&-\dfrac{\mu b_{r}\beta_{r}}{\mu+\varphi}\\\ \\\
\varphi&-\mu&0&-\dfrac{\varphi c_{i}\alpha_{i}}{\mu+\varphi}&-\dfrac{\varphi
c_{r}\alpha_{r}}{\mu+\varphi}\\\ \\\ 0&0&-(\mu+\omega)&\dfrac{\mu
b_{i}\beta_{i}+\varphi c_{i}\alpha_{i}}{\mu+\varphi}&\dfrac{\mu
b_{r}\beta_{r}+\varphi c_{r}\alpha_{r}}{\mu+\varphi}\\\ \\\
0&0&\omega&-(\gamma+\mu+\delta)&0\\\ \\\ 0&0&0&-\gamma&-\mu\end{pmatrix}$ (23)

The jacobian matrix have one column with diagonal entries.$-\mu$ are
eigenvalue of the matrix (23). Eliminating the second columns and the
corresponding row, we have:

$\displaystyle\begin{pmatrix}-(\mu+\varphi)&0&-\dfrac{\mu
b_{i}\beta_{i}}{\mu+\varphi}&-\dfrac{\mu b_{r}\beta_{r}}{\mu+\varphi}\\\
0&-(\mu+\omega)&\dfrac{\mu b_{i}\beta_{i}+\varphi
c_{i}\alpha_{i}}{\mu+\varphi}&\dfrac{\mu b_{r}\beta_{r}+\varphi
c_{r}\alpha_{r}}{\mu+\varphi}\\\ 0&\omega&-(\gamma+\mu+\delta)&0\\\
0&0&-\gamma&-\mu\end{pmatrix}$ (24)

Eliminating the first column, we have $-(\mu+\varphi$)as eigenvalue of the
matrix (24). So the Jacobian matrix become:

$\displaystyle\begin{pmatrix}-(\mu+\omega)&\dfrac{\mu b_{i}\beta_{i}+\varphi
c_{i}\alpha_{i}}{\mu+\varphi}&\dfrac{\mu b_{r}\beta_{r}+\varphi
c_{r}\alpha_{r}}{\mu+\varphi}\\\ \omega&-(\gamma+\mu+\delta)&0\\\
0&-\gamma&-\mu\end{pmatrix}$ (25)

Assume that:

$\displaystyle J_{1}$ $\displaystyle=$ $\displaystyle\mu+\omega$
$\displaystyle J_{2}$ $\displaystyle=$ $\displaystyle(\gamma+\mu+\delta)$
$\displaystyle J_{3}$ $\displaystyle=$ $\displaystyle\mu$ $\displaystyle A$
$\displaystyle=$ $\displaystyle\mu b_{i}\beta_{i}+\varphi c_{i}\alpha_{i}$
$\displaystyle B$ $\displaystyle=$ $\displaystyle\mu b_{r}\beta_{r}+\varphi
c_{r}\alpha_{r}$ $\displaystyle C$ $\displaystyle=$ $\displaystyle\mu+\varphi$

The remaining eigenvalues can be calculated by computing all the solutions of
the characteristic equation $|J_{dfe}-\lambda I|=0$
So,

$\displaystyle|J_{dfe}-\lambda
I|=\left|\begin{matrix}-(J_{1}+\lambda)&\dfrac{A}{C}&\dfrac{B}{C}\\\
\omega&-(J_{2}+\lambda)&0\\\ 0&\gamma&-(J_{3}+\lambda)\end{matrix}\right|$
(26)
$\displaystyle(J_{1}+\lambda)(J_{2}+\lambda)(J_{3}+\lambda)-\omega\gamma\dfrac{B}{C}-(J_{3}+\lambda)\omega\dfrac{A}{C}=0$
(27)

Thus, the equation (27) can be written as a polynomial of degree 3:

$\displaystyle\lambda^{3}+A_{2}\lambda^{2}+A_{1}\lambda+A_{0}=0$ (28)

Where

$\displaystyle A_{2}$ $\displaystyle=$ $\displaystyle J_{1}+J_{2}+J_{3}>0$
(29) $\displaystyle A_{1}$ $\displaystyle=$ $\displaystyle
J_{1}J_{2}+J_{1}J_{3}+J_{2}J_{3}-\omega\dfrac{A}{B}>0\mbox{ with
}(\omega\dfrac{A}{B}<J_{1}J_{2}+J_{1}J_{3}+J_{2}J_{3})$ $\displaystyle A_{0}$
$\displaystyle=$ $\displaystyle
J_{1}J_{2}J_{3}-\omega(J_{3}\dfrac{A}{C}+\gamma\dfrac{B}{C})$ $\displaystyle=$
$\displaystyle J_{1}J_{2}J_{3}C(1-R_{0})>0\mbox{ if }R_{0}<1$

According to Descarte’rules of sign [18],it has been observed that all the
coefficient of (28) are strictly positive if $R_{0}<1$. Then, all it’s
eigenvalue have a negative really part.
$E{dfe}$ is locally asymptotically stable by the Poincare Lyapnov theorem [19]
∎

####  3.4.1 Global stability of DFE

Theorem: The disease free equilibrium point, $E_{dfe}$ is globally
asymptotically stable if $R_{0}\leq 1$ and unstable if $R_{0}>1$.

###### Proof.

Let’s define

$L=(E,I,R)=aE+bI+cR$ (30)

as a Lyapunov function for some non negative values of $a$, $b$, and $c$. The
time derivative of $L$ is given by

$\displaystyle\dfrac{dL}{dt}$ $\displaystyle=$ $\displaystyle
a\dfrac{dE}{dt}+b\dfrac{dI}{dt}+c\dfrac{dR}{dt}$ $\displaystyle\leq$
$\displaystyle
a\Big{[}\dfrac{\mu}{\mu+\varphi}(b_{i}\beta_{i}I+b_{r}\beta_{r}R)+\dfrac{\varphi}{\mu+\varphi}(c_{i}\alpha_{i}I+c_{r}\alpha_{r}R)-(\mu+\omega)E\Big{]}$
$\displaystyle+b[\omega E-(\gamma+\mu+\delta)I]+c(\gamma I+\mu R)$
$\displaystyle=$
$\displaystyle[-a(\mu+\omega)+b\omega]E+\big{[}a(\dfrac{\mu}{\mu+\varphi}b_{i}\beta_{i}+\dfrac{\varphi}{\mu+\varphi}c_{i}\alpha_{i})-b(\gamma+\mu+\delta)+c\gamma\big{]}I+$
$\displaystyle\big{[}a(\dfrac{\mu}{\mu+\varphi}b_{r}\beta_{r}+\dfrac{\varphi}{\mu+\varphi}c_{r}\alpha_{r}I)-c\mu\big{]}R$

Computing the Lyapunov function coefficients so that the coefficients of $E$,
$I$, $R$, are all equal to zero yields to

$\displaystyle\left\\{\begin{array}[]{l}a=\dfrac{\omega}{\mu+\omega}\\\ \\\
b=\dfrac{\omega[\mu(\mu b_{i}\beta_{i}+\varphi c_{i}\alpha_{i})+\gamma(\mu
b_{r}\beta_{r}+\varphi
c_{r}\alpha_{r})]}{\mu(\mu+\varphi)(\mu+\omega)(\gamma+\mu+\delta)}\\\ \\\
c=\dfrac{\omega(\mu b_{r}\beta_{r}+\varphi
c_{r}\alpha_{r})}{\mu(\mu+\varphi)(\mu+\omega)}\end{array}\right.$ (36)

The time derivative of the Lyapunov function may be written using these
coefficients as

$\displaystyle\dfrac{dL}{dt}$ $\displaystyle=$
$\displaystyle\dfrac{\omega}{\mu+\omega}\dfrac{dE}{dt}+\dfrac{\omega[\mu(\mu
b_{i}\beta_{i}+\varphi c_{i}\alpha_{i})+\gamma(\mu b_{r}\beta_{r}+\varphi
c_{r}\alpha_{r})]}{\mu(\mu+\varphi)(\mu+\omega)(\gamma+\mu+\delta)}\dfrac{dI}{dt}+\dfrac{\omega(\mu
b_{r}\beta_{r}+\varphi
c_{r}\alpha_{r})}{\mu(\mu+\varphi)(\mu+\omega)}\dfrac{dR}{dt}$
$\displaystyle\leq$
$\displaystyle\dfrac{\omega}{\mu+\omega}\Big{[}\dfrac{\mu}{\mu+\varphi}(b_{i}\beta_{i}I+b_{r}\beta_{r}R)+\dfrac{\varphi}{\mu+\varphi}(c_{i}\alpha_{i}I+c_{r}\alpha_{r}R)-(\mu+\omega)E\Big{]}$
$\displaystyle+R_{0}[\omega E-(\gamma+\mu+\delta)]+\dfrac{\omega(\mu
b_{r}\beta_{r}+\varphi c_{r}\alpha_{r})}{\mu(\mu+\varphi)(\mu+\omega)}(\gamma
I+\mu R)$ $\displaystyle=$ $\displaystyle\omega(R_{0}-1)E$

when $R_{0}<1$ or $R_{0}=1$, it is evident that $\dfrac{dL}{dt}\leq 0$.
Furthermore $\dfrac{dL}{dt}=0$ iff $E=I=R=0$ and the largest compact invariant
set in $\omega$ such taht $\dfrac{dL}{dt}\leq 0$ is the $E_{dfe}$.
As a result, according to the Lyapunov-Lasalle Invariance Principle [20], the
disease free equilibrium is globally stable in $\Omega$ if $R_{0}\leq 0$

∎

###  3.5 Existence of endemic equilibrium

The endemic equilibrium (EE) is defined as the point at which the disease
remains present in the population. To compute it, we consider
$\dfrac{dS}{dt}=\dfrac{dV}{dt}=\dfrac{dE}{dt}=\dfrac{dI}{dt}=\dfrac{dR}{dt}=0$
and

$\displaystyle\left\\{\begin{array}[]{rcl}p-\lambda_{s}S^{*}-(\mu+\varphi)S^{*}&=&0\\\
\\\ \varphi S^{*}-\lambda_{v}V^{*}-\mu V^{*}&=&0\\\ \\\
\lambda_{s}S^{*}+\lambda_{v}V^{*}-(\mu+\omega)E^{*}&=&0\\\ \\\ \omega
E^{*}-(\gamma+\mu+\delta)I^{*}&=&0\\\ \\\ \gamma I^{*}-\mu
R^{*}&=&0\end{array}\right.$ (46)

We obtain

$\displaystyle\left\\{\begin{array}[]{rcl}S^{*}&=&\dfrac{p}{\lambda_{s}+(\mu+\varphi)}\\\
\\\ V^{*}&=&\dfrac{p\varphi}{[\lambda_{s}+(\mu+\varphi)](\lambda_{v}+\mu)}\\\
\\\ E^{*}&=&\dfrac{\lambda_{s}p(\lambda_{v}+\mu)+\lambda_{v}\varphi
p}{(\mu+\omega)[\lambda_{s}+(\mu+\varphi)](\lambda_{v}+\mu)}\\\ \\\
I^{*}&=&\dfrac{\omega[\lambda_{s}p(\lambda_{v}+\mu)+\lambda_{v}\varphi]}{(\mu+\omega)(\gamma+\mu+\delta)[\lambda_{s}+(\mu+\varphi)](\lambda_{v}+\mu)}\\\
\\\
R^{*}&=&\dfrac{\gamma\omega[\lambda_{s}p(\lambda_{v}+\mu)+\lambda_{v}\varphi
p]}{\mu(\mu+\omega)(\gamma+\mu+\delta)[\lambda_{s}+(\mu+\varphi)](\lambda_{v}+\mu)}\\\
\\\ \end{array}\right.$ (57)

###  3.6 Local stability of EE

We can calculate the stability of the endemic point using the eigenvalue
approach and evaluate it at $E^{*}$, but our system will require sophisticated
computations, which could lead to many errors. Then, at the DFE [21], we apply
the center Manifold theory created by Chavez and Song for bifurcation
analysis.
From system (10), let us rewritten the variables $S,V,E,I,R,$ in the
dimensionless states as $S=x_{1},V=x_{2},E=x_{3},I=x_{4},R=x_{5}$ and let
${X=(x_{1},x_{2},x_{3},x_{5}x_{5})^{T}}$, the expressions of $X$ in vector
notation.
The model system (10) can be rewritten in terms of
$\dfrac{dX}{dt}=(f_{1},f_{2},f_{3},\cdots f_{5})^{T}$
where

$\displaystyle\left\\{\begin{array}[]{rcl}\dfrac{dx_{1}}{dt}&=&p-\dfrac{x_{1}}{N}(b_{i}\beta_{i}x_{4}+b_{r}\beta_{r}x_{5})-(\mu+\varphi)x_{1}\\\
\\\ \dfrac{dx_{2}}{dt}&=&\varphi
x_{1}-\dfrac{x_{2}}{N}(c_{i}\alpha_{i}x_{4}+c_{r}\alpha_{r}x_{5})-\mu x_{2}\\\
\\\
\dfrac{dx_{3}}{dt}&=&\dfrac{x_{1}}{N}(b_{i}\beta_{i}x_{4}+b_{r}\beta_{r}x_{5})+\dfrac{x_{2}}{N}(c_{i}\alpha_{i}x_{4}+c_{r}\alpha_{r}x_{5})-(\mu+\omega)x_{3}\\\
\\\ \dfrac{dx_{4}}{dt}&=&\omega x_{3}-(\gamma+\mu+\delta)x_{4}\\\ \\\
\dfrac{dx_{5}}{dt}&=&\gamma x_{4}-\mu x_{5}\end{array}\right.$ (67)

and Where $N=x_{1}+x_{2}+x_{3}+x_{4}+x_{5}$

This can be used to show that when $R_{0}=1$, the bifurcation parameter is
$\beta_{i}=\beta_{i}^{*}$ , it follows that

$\displaystyle 1=\dfrac{\omega[\mu(\mu
b_{i}\beta_{i}+c_{i}\alpha_{i}\varphi)+\gamma(\mu
b_{r}\beta_{r}+c_{r}\alpha_{r}\varphi)]}{\mu(\mu+\varphi)(\mu+\omega)(\mu+\gamma+\delta)}$
(68)

so

$\displaystyle\beta_{i}^{*}=\dfrac{\mu(\mu+\varphi)(\mu+\omega)(\mu+\gamma+\delta)-\omega\Big{[}\gamma(\mu
b_{r}\beta_{r}+c_{r}\alpha_{r}\varphi)+c_{i}\alpha_{i}\varphi\Big{]}}{\omega\mu^{2}b_{i}}$
(69)

The center manifold theory will aid in the description of the system’s
dynamics at $\beta_{i}^{*}$, and the jacobian matrix of system (67) at the DFE
$E_{dfe}$ is provided by:

$\displaystyle J=\begin{pmatrix}-(\mu+\varphi)&0&0&-\dfrac{\mu
b_{i}\beta_{i}}{\mu+\varphi}&-\dfrac{\mu b_{r}\beta_{r}}{\mu+\varphi}\\\ \\\
\varphi&-\mu&0&-\dfrac{\varphi c_{i}\alpha_{i}}{\mu+\varphi}&-\dfrac{\varphi
c_{r}\alpha_{r}}{\mu+\varphi}\\\ \\\ 0&0&-(\mu+\omega)&\dfrac{\mu
b_{i}\beta_{i}+\varphi c_{i}\alpha_{i}}{\mu+\varphi}&\dfrac{\mu
b_{r}\beta_{r}+\varphi c_{r}\alpha_{r}}{\mu+\varphi}\\\ \\\
0&0&\omega&-(\gamma+\mu+\delta)&0\\\ \\\ 0&0&0&\gamma&-\mu\end{pmatrix}$ (70)

Assume that:
$J_{1}=\mu+\varphi$,
$J_{2}=\mu$,$J_{3}=\mu+\omega$,$J_{4}=\gamma+\mu+\delta$,$J_{5}=\mu$
The right eigenvector corresponding to the zero eigenvalue
$W=(w_{1},w_{2},w_{3},w_{4},w_{5})^{T}$, is computed using $J.W=0$ which gives
:

$\displaystyle J=\begin{pmatrix}-(\mu+\varphi)&0&0&-\dfrac{\mu
b_{i}\beta_{i}}{\mu+\varphi}&-\dfrac{\mu b_{r}\beta_{r}}{\mu+\varphi}\\\ \\\
\varphi&-\mu&0&-\dfrac{\varphi c_{i}\alpha_{i}}{\mu+\varphi}&-\dfrac{\varphi
c_{r}\alpha_{r}}{\mu+\varphi}\\\ \\\ 0&0&-(\mu+\omega)&\dfrac{\mu
b_{i}\beta_{i}+\varphi c_{i}\alpha_{i}}{\mu+\varphi}&\dfrac{\mu
b_{r}\beta_{r}+\varphi c_{r}\alpha_{r}}{\mu+\varphi}\\\ \\\
0&0&\omega&-(\gamma+\mu+\delta)&0\\\ \\\
0&0&0&\gamma&-\mu\end{pmatrix}\begin{pmatrix}w_{1}\\\ \\\ w_{2}\\\ \\\
w_{3}\\\ \\\ w_{4}\\\ \\\ w_{5}\\\ \\\ \end{pmatrix}$ (71)

we deduce

$\displaystyle\left\\{\begin{array}[]{l}-J_{1}w_{1}-\dfrac{\mu
b_{i}\beta_{i}}{\mu+\varphi}w_{4}-\dfrac{\mu
b_{r}\beta_{r}}{\mu+\varphi}w_{5}=0\\\ \varphi w_{1}-J_{2}w_{2}+\dfrac{\varphi
c_{i}\alpha}{\mu+\varphi}w_{4}-\dfrac{\varphi
c_{i}\alpha}{\mu+\varphi}w_{5}=0\\\ -J_{3}w_{3}+\dfrac{\mu
b_{i}\beta_{i}+\varphi c_{i}\alpha}{\mu+\varphi}w_{4}+\dfrac{\mu
b_{r}\beta_{r}+\varphi c_{i}\alpha}{\mu+\varphi}w_{5}=0\\\ \omega
w_{3}-J_{4}w_{4}=0\\\ \gamma w_{4}-J_{5}w_{5}=0\end{array}\right.$ (77)

By setting $W_{5}=1$ and solving (77) we get:

$\displaystyle\left\\{\begin{array}[]{rcl}w_{1}&=&-\dfrac{(\gamma
b_{i}\beta_{i}+\mu b_{r}\beta_{r})]}{(\mu+\varphi)^{2}}\\\
w_{2}&=&-\dfrac{\varphi\Big{[}\mu c_{i}\alpha_{i}+\gamma
c_{r}\alpha_{r}-\gamma(\mu+\varphi)^{]}}{\gamma\mu(\mu+\varphi)}\\\
w_{3}&=&\dfrac{\mu(\gamma+\mu+\delta)}{\gamma\omega}\\\
w_{4}&=&\dfrac{\mu}{\gamma}\\\ w_{5}&=&1\end{array}\right.$ (83)

$J$ has also a left eigenvector $V=(v_{1},v_{2},v_{3},v_{4},v_{5})$ associated
with the zero eigenvalue and satisfying $V.J=0$. To find this, we use the
jacobian matrix (70).
So,

$\displaystyle
J=\begin{pmatrix}v_{1}&v_{2}&w_{3}&v_{4}&v_{5}\end{pmatrix}\begin{pmatrix}-(\mu+\varphi)&0&0&-\dfrac{\mu
b_{i}\beta_{i}}{\mu+\varphi}&-\dfrac{\mu b_{r}\beta_{r}}{\mu+\varphi}\\\ \\\
\varphi&-\mu&0&-\dfrac{\varphi c_{i}\alpha_{i}}{\mu+\varphi}&-\dfrac{\varphi
c_{r}\alpha_{r}}{\mu+\varphi}\\\ \\\ 0&0&-(\mu+\omega)&\dfrac{\mu
b_{i}\beta_{i}+\varphi c_{i}\alpha_{i}}{\mu+\varphi}&\dfrac{\mu
b_{r}\beta_{r}+\varphi c_{r}\alpha_{r}}{\mu+\varphi}\\\ \\\
0&0&\omega&-(\gamma+\mu+\delta)&0\\\ \\\ 0&0&0&\gamma&-\mu\end{pmatrix}$ (84)

Which yields to

$\displaystyle\left\\{\begin{array}[]{l}-J_{1}v_{1}+\varphi v_{2}=0\\\
-J_{2}v_{2}=0\\\ -J_{3}v_{3}+\omega v_{4}=0\\\ -\dfrac{\mu
b_{i}\beta_{i}}{\mu+\varphi}v_{1}-\dfrac{\varphi
c_{i}\alpha_{i}}{\mu+\varphi}v_{2}+\dfrac{\mu b_{i}\beta_{i}+\varphi
c_{i}\alpha_{i}}{\mu+\varphi}v_{3}-J_{4}v_{4}-\gamma v_{5}=0\\\ -\dfrac{\mu
b_{r}\beta_{r}}{\mu+\varphi}v_{1}-\dfrac{\varphi
c_{r}\alpha_{r}}{\mu+\varphi}+\dfrac{\mu b_{r}\beta_{r}+\varphi
c_{r}\alpha_{r}}{\mu+\varphi}-J_{5}v_{5}=0\end{array}\right.$ (90)

The left eigenvector is computed and $v_{2}=0,v_{1}=0$,and setting $v_{3}=1$
we get:

$\displaystyle\left\\{\begin{array}[]{l}v_{1}=0\\\ v_{2}=0\\\ v_{3}=1\\\
v_{4}=\dfrac{\mu+\omega}{\omega}\\\ v_{5}=\dfrac{\mu b_{r}\beta_{r}+\varphi
c_{r}\alpha_{r}}{\mu(\mu+\varphi)}\end{array}\right.$ (96)

The theorem of Chavez and Song (2004) is repeated here for convenience, and it
may be used to demonstrate the local stability of the endemic equilibrium
point around $R_{0}=1$.

Theorem 2: Consider the following general system of ordinary differential
equations with a parameter $\beta^{*}_{i}$

$\displaystyle\dfrac{dx}{dt}=f(x,\beta^{*}_{i}),f:\mathbb{R}^{n}\times\mathbb{R}\rightarrow\mathbb{R}\mbox{
and }f\in\mathbb{C}^{2}(\mathbb{R}^{n}\times\mathbb{R})$

where $0$ is an equilibrium point of the system, (that is
$f(0.\beta^{*}_{i})=0$ for all $\beta^{*}_{i}$ and assume:
$A_{1}:J=D_{x}f(0,0)=\dfrac{df_{i}}{dx_{i}}(0,0)$ is the linearization matrix
of System (70) around the equilibrium point. Zero is a simple eigenvalue of
$A$ and all other eigenvalues of $J$ have negative real parts;
$A2$ : Matrix $J$ has a nonnegative right eigenvector $w$ and a left
eigenvector $v$ corresponding to the zero eigenvalue.

Let $f_{k}$ be the $k^{th}$ component of $f$ and

$\displaystyle
a=\sum^{n}_{k,i,j=1}v_{k}w_{i}w_{j}\dfrac{\partial^{2}f_{k}}{\partial
x_{i}\partial x_{j}}(0,0)$ (97) $\displaystyle
b=\sum^{n}_{k,i=1}v_{k}w_{i}\dfrac{\partial^{2}f_{k}}{\partial
x_{i}\partial\beta^{*}_{i}}(0,0)$ (98)

Then the local dynamics of the system around the equilibrium piont $0$ are
totally determined by the sign of $a$ and $b$ .

  1. 1.

$a>0$, $b>0$. When $\varphi<0$ with $|\beta_{i}^{*}|<<1$, 0 is locally
asymptotically stable, and there exists a positive unstable equilibrium; when
$0<\beta^{*}_{i}<<1$, $0$ is unstable and there exists a negative and locally
asymptotically stable equilibrium;

  2. 2.

$a<0$, $b<0$. When $\beta_{i}^{*}<0$ with $|\beta_{i}^{*}|<<1$, $0$ is
unstable; when $0<\beta_{i}^{*}<<1$, $0$ is locally asymptotically stable, and
there exists a positive unstable equilibrium;

  3. 3.

$a>0$, $b<0$. When $\beta_{i}^{*}<0$ with $|\beta_{i}^{*}|<<1$, $0$ is
unstable, and there exists a locally asymptotically stable negative
equilibrium; when $0<\beta_{i}^{*}<<1$, $0$ is stable, and a positive unstable
equilibrium appears;

  4. 4.

$a<0$, $b>0$. When $\beta_{i}^{*}$ changes from negative to positive, $0$
changes its stability from stable to unstable. Correspondingly a negative
unstable equilibrium becomes positive and locally asymptotically stable.
Particularly, if $a>0$ and $b>0$, then, a backward bifurcation occurs at
$\beta_{i}^{*}=0.$

Computation of $a$ and $b$

Since $v_{1}=v_{2}=0$ for $k=1,2$, we consider only $k=3,4,5$, the values of
$a$ and $b$ are obtained from:

$\displaystyle a$ $\displaystyle=$
$\displaystyle\sum_{i,j=1}^{5}v_{3}w_{i}w_{j}\dfrac{\partial^{2}f_{3}}{\partial
x_{i}\partial
x_{j}}+\sum_{i,j=1}^{5}v_{4}w_{i}w_{j}\dfrac{\partial^{2}f_{4}}{\partial
x_{i}\partial
x_{j}}+\sum_{i,j=1}^{5}v_{5}w_{i}w_{j}\dfrac{\partial^{2}f_{5}}{\partial
x_{i}\partial x_{j}}$ (99) $\displaystyle b$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{5}v_{3}w_{i}\dfrac{\partial^{2}f_{3}}{\partial
x_{i}\partial\beta_{i}^{*}}+\sum_{i=1}^{5}v_{4}w_{i}\dfrac{\partial^{2}f_{4}}{\partial
x_{i}\partial\beta_{i}^{*}}+\sum_{i=1}^{5}v_{5}w_{i}\dfrac{\partial^{2}f_{5}}{\partial
x_{i}\partial\beta_{i}^{*}}$ (100)

The corresponding non-zero second order partial derivatives at disease-free
equilibrium for system (67) are provided by:

$\displaystyle\dfrac{\partial^{2}f_{3}}{\partial x_{1}\partial x_{4}}$
$\displaystyle=$
$\displaystyle\dfrac{\varphi\mu}{p(\mu+\varphi)}(b_{i}\beta_{i}-c_{i}\alpha_{i})$
(101) $\displaystyle\dfrac{\partial^{2}f_{3}}{\partial x_{1}\partial x_{5}}$
$\displaystyle=$
$\displaystyle\dfrac{\varphi\mu}{p(\mu+\varphi)}(b_{r}\beta_{r}-c_{r}\alpha_{r})$
(102) $\displaystyle\dfrac{\partial^{2}f_{3}}{\partial x_{2}\partial x_{4}}$
$\displaystyle=$
$\displaystyle-\dfrac{\mu^{2}}{p(\mu+\varphi)}(b_{i}\beta_{i}-c_{i}\alpha_{i})$
(103) $\displaystyle\dfrac{\partial^{2}f_{3}}{\partial x_{2}\partial x_{5}}$
$\displaystyle=$
$\displaystyle-\dfrac{\mu^{2}}{p(\mu+\varphi)}(b_{r}\beta_{r}-c_{r}\alpha_{r})$
(104) $\displaystyle\dfrac{\partial^{2}f_{3}}{\partial x_{4}\partial x_{5}}$
$\displaystyle=$ $\displaystyle-\dfrac{\mu}{p(\mu+\varphi)}(\mu
b_{i}\beta_{i}-\varphi c_{i}\alpha_{i})-(\mu b_{r}\beta_{r}-\varphi
c_{r}\alpha_{r})$ (105) $\displaystyle\dfrac{\partial^{2}f_{3}}{\partial
x_{4}\partial\beta^{*}_{i}}$ $\displaystyle=$
$\displaystyle\dfrac{\mu}{\mu+\varphi}b_{i}$ (106)

By replacing (101), (102), (103), (104), (105) in $a$ (99)

$\displaystyle a$ $\displaystyle=$ $\displaystyle
2w_{1}w_{4}\dfrac{\partial^{2}f_{3}}{\partial x_{1}\partial
x_{4}}+2w_{1}w_{5}\dfrac{\partial^{2}f_{3}}{\partial x_{1}\partial
x_{5}}+2w_{2}w_{4}\dfrac{\partial^{2}f_{3}}{\partial x_{2}\partial
x_{4}}+2w_{2}w_{5}\dfrac{\partial^{2}f_{3}}{\partial x_{2}\partial
x_{5}}+2w_{4}w_{5}\dfrac{\partial^{2}f_{3}}{\partial x_{4}\partial x_{5}}$
$\displaystyle=$ $\displaystyle-\dfrac{2\mu^{2}\varphi(\gamma
b_{i}\beta_{i}+\mu
b_{r}\beta_{r})}{p(\mu+\varphi)^{3}}(b_{i}\beta_{i}-c_{i}\alpha_{i})-\dfrac{2\varphi\mu(\gamma
b_{i}\beta_{i}+\mu
b_{r}\beta_{r})}{p(\mu+\varphi)^{3}}(b_{r}\beta_{r}-c_{r}\alpha_{r})$
$\displaystyle+\dfrac{2\varphi\mu^{2}[\mu c_{i}\alpha_{i}+\gamma
c_{r}\alpha_{r}-\gamma(\mu+\varphi)^{2}]}{\gamma^{2}p(\mu+\varphi)^{2}}(b_{r}\beta_{r}-c_{r}\alpha_{r})$
$\displaystyle-\dfrac{2\mu^{2}}{\gamma p(\mu+\varphi)}[(\mu
b_{i}\beta_{i}+\varphi c_{i}\alpha_{i})+\mu(b_{r}\beta_{r}+\varphi
c_{r}\alpha_{r})]$ $\displaystyle=$ $\displaystyle-2(\gamma b_{i}\beta_{i}+\mu
b_{r}\beta_{r})\dfrac{\varphi\mu}{p(\mu+\varphi)^{3}}\Big{[}\dfrac{\mu}{\gamma}(b_{i}\beta_{i}-c_{i}\alpha_{i})+(b_{r}\beta_{r}-c_{r}\alpha_{r})\Big{]}$
$\displaystyle-2\mu\varphi\Big{[}\dfrac{\mu c_{i}\alpha_{i}+\gamma
c_{r}\alpha_{r}-\gamma(\mu+\varphi)^{2}}{\gamma
p(\mu+\varphi)^{2}}\Big{]}\Big{[}\dfrac{\mu}{\gamma}(b_{i}\beta_{i}-c_{i}\alpha_{i})+(b_{r}\beta_{r}-c_{r}\alpha_{r})\Big{]}$
$\displaystyle-\dfrac{2\mu^{2}}{\gamma p(\mu+\varphi)}[(\mu
b_{i}\beta_{i}+\varphi c_{i}\alpha_{i})+(\mu b_{r}\beta_{r}+\varphi
c_{r}\alpha_{r})]$ $\displaystyle=$
$\displaystyle\dfrac{2\varphi\mu^{2}}{\gamma
p(\mu+\varphi)^{3}}c_{i}\alpha_{i}(\gamma b_{i}\beta_{i}+\mu
b_{r}\beta_{r})+\dfrac{2\varphi\mu}{p(\mu+\varphi)^{3}}c_{r}\alpha_{r}(\gamma
b_{i}\beta_{i}+\mu b_{r}\beta_{r})$
$\displaystyle+\dfrac{2\varphi\mu^{2}}{\gamma^{2}p(\mu+\varphi)^{2}}b_{i}\beta_{i}(\mu
c_{i}\alpha_{i}+\gamma c_{r}\alpha_{r})+\dfrac{2\varphi\mu}{\gamma
p(\mu+\varphi)^{2}}b_{r}\beta_{r}(\mu c_{i}\alpha_{i}+\gamma
c_{r}\alpha_{r})+\dfrac{2\varphi\mu^{2}}{\gamma
p}c_{i}\alpha_{i}+\dfrac{2\varphi\mu}{p}c_{r}\alpha_{r}$
$\displaystyle-\dfrac{2\varphi\mu^{2}}{\gamma
p(\mu+\varphi)^{3}}b_{i}\beta_{i}(\gamma b_{i}\beta_{i}+\mu
b_{r}\beta_{r})-\dfrac{2\varphi\mu}{p(\mu+\varphi)^{3}}b_{r}\beta_{r}(\gamma
b_{i}\beta_{i}+\mu b_{r}\beta_{r})$
$\displaystyle-\dfrac{2\varphi\mu^{2}}{\gamma^{2}p(\mu+\varphi)^{2}}c_{i}\alpha_{i}(\mu
c_{i}\alpha_{i}+\gamma c_{r}\alpha_{r})-\dfrac{2\varphi\mu}{\gamma
p(\mu+\varphi)^{2}}c_{r}\alpha_{r}(\mu c_{i}\alpha_{i}+\gamma
c_{r}\alpha_{r})-\dfrac{2\varphi\mu^{2}}{\gamma p}b_{i}\beta_{i}$
$\displaystyle-\dfrac{2\varphi\mu}{p}b_{r}\beta_{r}-\dfrac{2\mu^{2}}{\gamma
p(\mu+\varphi)}[(\mu b_{i}\beta_{i}+\varphi c_{i}\alpha_{i})+(\mu
b_{r}\beta_{r}+\varphi c_{r}\alpha_{r})]$ $\displaystyle=$
$\displaystyle\dfrac{2\varphi\mu}{\gamma p(\mu+\varphi)^{3}}(\gamma
b_{i}\beta_{i}+\mu b_{r}\beta_{r})(\mu c_{i}\alpha_{i}+\gamma
c_{r}\alpha_{r})+\dfrac{2\varphi\mu}{\gamma^{2}p(\mu+\varphi)^{2}}(\mu
c_{i}\alpha_{i}+\gamma c_{r}\alpha_{r})(\mu b_{i}\beta_{i}+\gamma
b_{r}\beta_{r})$ $\displaystyle+\dfrac{2\varphi\mu}{\gamma p}(\mu c_{i}\alpha
c_{i}+\gamma c_{r}\alpha_{r})-\dfrac{2\varphi\mu}{\gamma
p(\mu+\varphi)^{3}}(\gamma b_{i}\beta_{i}+\mu b_{r}\beta_{r})(\mu
b_{i}\beta_{i}+\gamma b_{r}\beta_{r})$
$\displaystyle-\dfrac{2\varphi\mu}{\gamma^{2}p(\mu+\varphi)^{2}}(\mu
c_{i}\alpha_{i}+\gamma c_{r}\alpha_{r})^{2}-\dfrac{2\varphi\mu}{\gamma p}(\mu
b_{i}\beta_{i}+\gamma b_{r}\beta_{r})$ $\displaystyle-\dfrac{2\mu^{2}}{\gamma
p(\mu+\varphi)}[(\mu b_{i}\beta_{i}+\varphi c_{i}\alpha_{i})+(\mu
b_{r}\beta_{r}+\varphi c_{r}\alpha_{r})]$

So, $a=A-B$ with

$\displaystyle A$ $\displaystyle=$ $\displaystyle\dfrac{2\varphi\mu}{\gamma
p(\mu+\varphi)^{3}}(\gamma b_{i}\beta_{i}+\mu b_{r}\beta_{r})(\mu
c_{i}\alpha_{i}+\gamma
c_{r}\alpha_{r})+\dfrac{2\varphi\mu}{\gamma^{2}p(\mu+\varphi)^{2}}(\mu
c_{i}\alpha_{i}+\gamma c_{r}\alpha_{r})(\mu b_{i}\beta_{i}+\gamma
b_{r}\beta_{r})$ $\displaystyle+\dfrac{2\varphi\mu}{\gamma p}(\mu c_{i}\alpha
c_{i}+\gamma c_{r}\alpha_{r})$ $\displaystyle B$ $\displaystyle=$
$\displaystyle\dfrac{2\varphi\mu}{\gamma p(\mu+\varphi)^{3}}(\gamma
b_{i}\beta_{i}+\mu b_{r}\beta_{r})(\mu b_{i}\beta_{i}+\gamma
b_{r}\beta_{r})-\dfrac{2\varphi\mu}{\gamma^{2}p(\mu+\varphi)^{2}}(\mu
c_{i}\alpha_{i}+\gamma c_{r}\alpha_{r})^{2}$
$\displaystyle-\dfrac{2\varphi\mu}{\gamma p}(\mu b_{i}\beta_{i}+\gamma
b_{r}\beta_{r})-\dfrac{2\mu^{2}}{\gamma p(\mu+\varphi)}[(\mu
b_{i}\beta_{i}+\varphi c_{i}\alpha_{i})+(\mu b_{r}\beta_{r}+\varphi
c_{r}\alpha_{r})]$

Possibilities of the sign of a

  * 1.

Case 1: if $A>B$, then $a>0$

  * 2.

Case 2: if $A<B$, then $a<0$

By replacing (106) in $b$ (100)

$\displaystyle b$ $\displaystyle=$ $\displaystyle
2v_{3}w_{4}\dfrac{\partial^{2}f_{3}}{\partial x_{4}\partial\beta_{i}^{*}}$
$\displaystyle=$ $\displaystyle\dfrac{2\mu^{2}}{\gamma(\mu+\varphi)}b_{i}>0$

Lemma If $a>0$ then the model (10) undergoes a backward bifurcation at $R_{0}$
close to $1$, if not then the endemic equilibrium point is locally
asymptotically stable for $R_{0}>1$.

In addition, when $R_{0}<1$, the existence of a backward bifurcation means
that there is a potential coexistence of an endemic equilibrium and a disease-
free equilibrium. In this case, the strategy of reducing the basic
reproduction number to a value less than unity would not be sufficient for the
eradication of Ebola disease even with high coverage of vaccination. Such
situation will have an impact on herd immunity as we will see below.

##  4 Discussion

In this paper an ordinardinary differential equation (ODE) model for the
transmission of EVD, with 5 variables(S,V,E,I,R) was analyzed. The model takes
into account a changing total human population over the time, it includes new
recruits in the susceptible class. We proved that the model is mathematically
and epidemiologically well posed, we also proved the existence of the disease
free equilibrium and endemic equilibrium points. The basic reproduction
number, $R_{0}$, was used to assess the stability of the disease-free steady
state and the presence of endemic equilibria.

We showed that when $R_{0}<1$, then the DFE is locally asymptotically stable
and unstable when $R_{0}>1$. Our model undergoes a backward bifurcation at
$R_{0}$ close to $1$ which shows a possible cohabitation between the EE and
the DFE, In general this shows us that it can not be possible to prove that
the EE is unique and stable for $R_{0}>1$. We noticed that $R_{0}$ cannot be
considered as a modified control intervention measure in our model regardless
of $R_{0}<1$, the disease will not be eradicated.

## References

  * [1] AA Chughtai, M Barnes, and CR Macintyre.  Persistence of ebola virus in various body fluids during convalescence: evidence and implications for disease transmission and control.  Epidemiology & Infection, 144(8):1652–1660, 2016.
  * [2] Anna Thorson, Pierre Formenty, Clare Lofthouse, and Nathalie Broutet.  Systematic review of the literature on viral persistence and sexual transmission from recovered ebola survivors: evidence and recommendations.  BMJ open, 6(1):e008859, 2016.
  * [3] Jessica L Abbate, Carmen Lia Murall, Heinz Richner, and Christian L Althaus.  Potential impact of sexual transmission on ebola virus epidemiology: Sierra leone as a case study.  PLoS neglected tropical diseases, 10(5):e0004676, 2016.
  * [4] Boubacar Diallo, Daouda Sissoko, Nicholas J Loman, Hadja Aïssatou Bah, Hawa Bah, Mary Claire Worrell, Lya Saidou Conde, Ramata Sacko, Samuel Mesfin, Angelo Loua, et al.  Resurgence of ebola virus disease in guinea linked to a survivor with virus persistence in seminal fluid for more than 500 days.  Clinical infectious diseases, 63(10):1353–1356, 2016.
  * [5] Lorenzo Subissi, Mory Keita, Samuel Mesfin, Giovanni Rezza, Boubacar Diallo, Steven Van Gucht, Emmanuel Onuche Musa, Zabulon Yoti, Sakoba Keita, Mamoudou Harouna Djingarey, et al.  Ebola virus transmission caused by persistently infected survivors of the 2014–2016 outbreak in west africa.  The Journal of infectious diseases, 218(suppl_5):S287–S291, 2018\.
  * [6] Prevail III Study Group.  A longitudinal study of ebola sequelae in liberia.  New England Journal of Medicine, 380(10):924–934, 2019.
  * [7] The new case of Ebola Virus Disease in Beni is linked to the 2018-2020 Nord Kivu/Ituri EVD outbreak and does not represent a new spillover event.  https://virological.org/t/oct-2021-evd-case-in-drc-linked-to-2018-2020-nord-kivu-evd-outbreak/762, 2021\.  Accessed: 2021-10-26.
  * [8] Alpha Kabinet Keita, Fara R Koundouno, Martin Faye, Ariane Düx, Julia Hinzmann, Haby Diallo, Ahidjo Ayouba, Frederic Le Marcis, Barré Soropogui, Kékoura Ifono, et al.  Resurgence of ebola virus in 2021 in guinea suggests a new paradigm for outbreaks.  Nature, 597(7877):539–543, 2021.
  * [9] Irina Tiper, Moussa Kourout, Carolyn Fisher, Krishnamurthy Konduru, Anjan Purkayastha, Gerardo Kaplan, and Robert Duncan.  Tracking ebolavirus genomic drift with a resequencing microarray.  Plos one, 17(2):e0263732, 2022.
  * [10] Kai Kupferschmidt.  New ebola outbreak likely sparked by a person infected 5 years ago.  Science, 372:1342–1344, 2021.
  * [11] Ashma Surani, Yogesh S Marfatia, Ashutosh Pal, and Ruchi Shah.  Ebola virus: An emerging sexually transmissible infection pathogen.  Indian Journal of Sexually Transmitted Diseases and AIDS, 39(1):65, 2018.
  * [12] N Anggriani, AK Supriatna, and E Soewono.  The existence and stability analysis of the equilibria in dengue disease infection model.  In Journal of Physics: Conference Series, volume 622, page 012039\. IOP Publishing, 2015.
  * [13] Muhammad Rafiq, Waheed Ahmad, Mujahid Abbas, and Dumitru Baleanu.  A reliable and competitive mathematical analysis of ebola epidemic model.  Advances in Difference Equations, 2020(1):1–24, 2020.
  * [14] Harout Boujakjian.  Modeling the spread of ebola with seir and optimal control.  SIAM Undergraduate Research Online, 9:299–310, 2016.
  * [15] Daniel G Bausch, Jonathan S Towner, Scott F Dowell, Felix Kaducu, Matthew Lukwiya, Anthony Sanchez, Stuart T Nichol, Thomas G Ksiazek, and Pierre E Rollin.  Assessment of the risk of ebola virus transmission from bodily fluids and fomites.  The Journal of infectious diseases, 196(Supplement_2):S142–S147, 2007.
  * [16] Chinwendu E Madubueze, Anande R Kimbir, and Terhemen Aboiyar.  Global stability of ebola virus disease model with contact tracing and quarantine.  Applications and Applied Mathematics: An International Journal (AAM), 13(1):25, 2018.
  * [17] Muhammad Tahir, Syed Inayat Ali Shah, Gul Zaman, and Sher Muhammad.  Ebola virus epidemic disease its modeling and stability analysis required abstain strategies.  Cogent Biology, 4(1):1488511, 2018.
  * [18] Bruce Anderson, Jeffrey Jackson, and Meera Sitharam.  Descartes’ rule of signs revisited.  The American Mathematical Monthly, 105(5):447–451, 1998.
  * [19] Chris P Tsokos.  On the classical stability theorem of poincaré-lyapunov with a random parameter.  Proceedings of the Japan Academy, 45(9):781–785, 1969.
  * [20] Gang Huang, Yasuhiro Takeuchi, Wanbiao Ma, and Daijun Wei.  Global stability for delay sir and seir epidemic models with nonlinear incidence rate.  Bulletin of mathematical biology, 72:1192–1207, 2010.
  * [21] F Berezovskaya, Georgy Karev, Baojun Song, and Carlos Castillo-Chavez.  A simple epidemic model with surprising dynamics.  Mathematical Biosciences & Engineering, 2(1):133–152, 2004.