In the current environment of economic uncertainty related to the spread of COVID-19, it could be of interest to develop financial instruments that allow market participants to express their views on the economic impact of an epidemic. The simplest such instruments are infection options, which are options on the fraction of infected individuals in a population affected by an epidemic 1 . Infection options would allow dynamic hedging of risk exposure to pandemic linked catastophe bonds. Investors purchasing catastrophe bonds provide insurance to issuers and lose invested princpal if certain conditions are met 2 .
In this note we propose a mathematical model for pricing and risk management of such options. One of the challenges is to formulate a model which, at the very least, exhibits typical stylized facts about the spread of highly infectious 1 According to https://www.worldometers.info/, as of March 17, the fractions of populations currently infected by COVID-19 is 0.0018% for the USA, 0.0313% for Switzerland, and 0.0431% for Italy. 2 In June 2017 the World Bank issued pandemic linked catastrophe bonds and derivatives worth $425 million. diseases such as COVID-19, while remaining numerically tractable, intuitive, and interpretable. The basic category of models of mathematical epidemiology [1] are compartmental models, which divide populations into homogeneous compartments. The model we selected is a stochastic differential equation extension of a member of that family of models, namely the Susceptible-Infected-Recovered (SIR) model.
The SIR model of infectious diseases is among the simplest models, and was originally published in 1927 by Kermack and McKendrick [6] . The classic SIR model is formulated in terms of three state variables: The model dynamics is given as the following dynamical system:
with the initial condition:
The constant parameters β > 0 and γ > 0 are called the infection and recovery rates, respectively. Notice that this dynamics obeys the conservation law
consistent with the assumption that the variables x, y, and z represent population fractions. This means that the variable z is, in a way, redundant, as its current value does not affect the dynamics of x and y, and it can be computed in a straightforward manner from (3).
Stochastic extensions of the SIR model go back to Bartlett [3] , who formulated a stochastic jump process describing the evolution of an epidemic. Various other stochastic extensions have been proposed since, see e.g. [2] . Here, we will focus on two stochastic extensions of the SIR model, which utilize stochastic differential equations, see e.g. [5] and references therein. These models are particularly adept at option modeling. The first of the models, a one factor stochastic SIR model is obtained in the following manner. Let W t denote the standard Brownian motion, and letẆ t denote the (generalized) white noise process. We assume that the infection rate β, rather than being constant, is subject to random shocks, namely β t = β + σẆ t , while the recovery rate γ remains constant. Here σ > 0 is a constant volatility parameter. This leads to the following stochastic dynamics:
with the initial condition analogous to (2) . Notice that the conservation law
continues to hold in the stochastic model, and so we can focus attention on the dynamics described by X and Y only. Consider now a time horizon T > 0, and ξ, η > 0. The backward Kolmogorov equation for the Green's function G = G T,ξ,η (t, x, y) corresponding to the stochastic system (4) reads:
We assume that G T,ξ,η (t, x, y) satisfies the following terminal condition:
Similarly, a two factor stochastic SIR model can be formulated in the following manner. In addition to allowing the infection rate being stochastic, we assume that the recovery rate is subject to random shocks delivered by a second Brownian motion B t , γ t = γ + ζḂ t . Again, ζ > 0 is a volatility parameter. The two Brownian motions are not required to be independent, and we let ρ denote the correlation coefficient between them,
This leads to the following stochastic dynamics:
with the initial condition analogous to (2) . The two factor stochastic model continues to obey the conservation law (5) . Notice that the two factor model has two additional parameters compare to the one factor model, namely the volatility of the recovery rate ζ and the correlation coefficient ρ in (8).
The backward Kolmogorov equation for the stochastic system (9) reads as follows. Set
ThenĠ − βxy∇ x G + (βx − γ)y∇ y G
where the Green's function G = G T,ξ,η (t, x, y) satisfies terminal condition (7).
The payoff of the call option on the infection fraction struck at K is given by (Y T − K) + multiplied by the notional amount on the contract. Here, as usual, 0) . Likewise, the payoff of the put option on the infection fraction struck at K is (K − Y T ) + multiplied by the notional amount on the contract Infection options are written on an underlying which is not a financial asset. For this reason, they are much closer in spirit to weather options (or real options), than to traditional equity or FX options. Consequently, there is no natural risk neutral approach to their valuation.
For this reason, we adopt an approach that is based on the process (4) (or, alternatively, the more complicated process (9)), which is written under the physical measure. Today's value of an infection call option is thus given by
where r is the discount rate, and E denotes expected value with respect to the probability measure generated by the process (4) (or (9)). In terms of the Green's function G, this can be expressed as a two dimensional integral over the terminal probability distribution:
Note that while the current value of the option depends on the current reading of the susceptible fraction, its value at expiration is not set. Analogous formulas hold for put options. Closed form solutions to equations (6) and (11) are not available and the evaluation of the expected values requires numerical computations. Two methods are of practical importance.
Substituting t → T − t we convert the (two-dimensional) backward Kolmogorov equation to a diffusion equation, which can be solved through a number of standard methods, such as the ADI method, see e.g. [7] .
Alternatively, one can calculate the expected values using Monte Carlo simulations [4] . Specifically, a convenient discretization to (4) is formulated as follows. In order to guarantee the positivity of X t and Y t throughout the simulation, we set X t = exp(−l t ) and Y t = exp(−m t ). Denoting by δ the magnitude of the time step, we are lead to the following Euler scheme:
for n = 0, 1, . . .. Here, z n are independent variates drawn from the standard normal distribution. At each iteration step one has to floor l n+1 and m n+1 at 0, l n+1 = max(l n+1 , 0), m n+1 = max(m n+1 , 0), so that X t , Y t ≤ 1. Similarly, the two factor system (9) can be discretized as follows:
l n+1 = l n + (βδ + σ √ δ z n )e −mn + 1 2 σ 2 δe −2mn , m n+1 = m n + γδ − (βδ + σ √ δ z n )e −ln + ζ √ δ b n + 1 2 η(e −ln ) 2 δ,
with η(x) defined in (10), for n = 0, 1, . . .. The two dimensional random normal vectors (z n , b n ) T are independently drawn from the Gaussian distribution with covariance