Continuously monitored barrier options under Markov processes
Abstract.
In this paper we present an algorithm for pricing barrier options in onedimensional Markov models. The approach rests on the construction of an approximating continuoustime Markov chain that closely follows the dynamics of the given Markov model. We illustrate the method by implementing it for a range of models, including a local Lévy process and a local volatility jumpdiffusion. We also provide a convergence proof and error estimates for this algorithm.
1. Introduction
1.1. Background and motivation
Barrier options are among the most popular exotic derivatives. Such contracts form effective risk management tools, and are liquidly traded in the Foreign Exchange markets. The most liquid barrier options in FX markets are continuously monitored single or doublenotouch options and knockin or knockout calls and puts (see e.g. Hakala and Wystup [28], Lipton [45], [46], Wystup [56]). The main challenge in the risk management of large portfolios of barrier options faced by trading desks that make a market in these securities is to be able to price and hedge the barrier products in models that are flexible enough to describe the observed option prices (i.e. calibrate to the vanilla market price quotes).
It is by now well established that the classical BlackScholes model lacks the flexibility to fit accurately to observed option price data (see e.g. Gatheral [23] and the references therein). A variety of models have been proposed to provide an improved description of the dynamics of the price of the underlying that can more accurately describe the option surface. Parametric diffusion models like the CEV process [14] have additional flexibility to fit the vanilla skew at a single maturity for as many options as there are free parameters in the model. The seminal idea (developed by Dupire [18] and Gyöngy [26]) that allows one to construct a model that can describe the entire implied volatility surface (across all strikes and maturities) is that of local volatility models, where a nonparametric form of the local volatility function is constructed from the option price data. It has been shown that in practice such models imply unrealistic dynamics of the option prices (see the formula for the implied volatility in a local volatility model given in [27]). The ramification is an unrealistic amount of vega risk, which is expensive to hedge. Therefore, even though in a local volatility model barrier options can be priced using a PDE solver, this modelling framework alone is not suitable for the risk management of a large portfolio of barrier options.
At the other end of the spectrum are the jump processes with stationary and independent increments, which can fit very well the volatility smile at a single maturity (see e.g. [12] and the references therein). A variety of models in the exponential Lévy class have been proposed in the literature: CGMY [8], KoBoL [6], generalised hyperbolic [19], NIG [4] and Kou [37]. Exponential Lévy processes are simple examples of Markov processes whose law is uniquely determined by the distribution of the process at a single time. Since the set of call option prices at a fixed maturity for all strikes uniquely determines the marginal riskneutral distribution at that maturity, calibration to option prices at multiple maturities in principle fixes the corresponding marginals. It has been reported (see e.g. [12]) that Lévy processes lack the flexibility of calibrating simultaneously across a range of strikes and maturities. Several generalisations within the onedimensional Markov framework have been proposed.
If the stationarity assumption is relaxed while the property of independence of increments is retained, one arrives at the class of exponential additive processes, which have recently been shown to calibrate well to several maturities in equity markets. The Sato process introduced in Carr et al. [10] is an example of such an additive model used in financial modelling.
The independent increments property of a process implies that its transition probabilities are translation invariant in space, so that they only depend on the difference between the end and starting value of the process. It is well known that the distribution of a logasset price depends in a nonlinear way on the starting point (e.g. in equity markets it has been observed that if the current price is high, then the volatility is low and vice versa). To capture this effect one is led to consider Markov jumpprocesses whose increments are not independent. As a generalisation of local volatility models, the class of local Lévy processes introduced by Carr et al. [9] allows the modeller to modulate the intensity of the jumps as well as their distribution depending on where the underlying asset is trading. A local volatility jumpdiffusion with similar structural properties was calibrated to the implied volatility surface in Andersen and Andreasen [2] and He et al. [29]. Due to the presence of jumps and the absence of stationarity and independence of increments, the problem of obtaining the firstpassage probabilities for such a general class of processes is computationally less tractable.
There exists currently a good deal of literature on numerical methods for the pricing of barriertype options. It is well known that in this case a straightforward Monte Carlo simulation algorithm will be timeconsuming and yield unstable results for the prices and especially the sensitivities. The knockin/out features in the barrier option payoffs lead to slower convergence of the Monte Carlo algorithm. To address this problem the following (semi)analytical approaches have been developed for specific models:
The method (a) employs the explicit spectral decompositions for this class of diffusion models, whereas the approach (b) exploits the independence and stationarity of the increments of the Lévy process, and the socalled WienerHopf factorisation. Since both of these approaches hinge on special structural properties of the underlying processes, it is not clear if and how they can be extended to more general Markovian models.
A different approach, pioneered by Kushner (see e.g. [41]), is the discrete time Markov chain approximation method. Originally developed for the numerical solution of stochastic optimal control problems in continuous time, this method consists of approximating the system of interest by a discrete time chain that closely follows its dynamics, and solving the problem of interest for this chain. An application to the pricing of American type options is given in Kushner [40]. Bally et al. [3] develop a quantization method to efficiently value American options on baskets of assets under a local volatility model. Using a discrete time Markov chain, Duan et al. [16] price a discretely monitored barrier option in the BlackScholes and NGARCH models. Rogers and Stapleton [51] investigate an efficient binomial tree method for barrier option pricing (see also references therein for related methods). Related are PDE and PIDE finite difference discretization methods that have been investigated by various authors; Zvan et al. [57] consider barrier and related options in the BlackScholes model; Tavella and Randall [55] present an overview of PDE finite difference methods for the pricing of financial instruments; Wang et al. [32] develop a robust scheme for American options under a CGMY model, and Cont and Voltchkova [13] follow a viscosity approach for the PIDEs connected to European and barrier option under Lévy models. Markov chains have also been employed to directly model the evolution of price processes; Albanese and Mijatović [1] model the stochasticity of risk reversals and carried out a calibration study in FX markets under a certain continuoustime Markov chain constructed to model the FX spot process.
1.2. Contribution of the current paper
In this paper we consider the problem of pricing barrier options in the setting of onedimensional Markov processes, which in particular includes the case of local Lévy models as well as local volatility jumpdiffusions and additive processes. The presented approach is probabilistic in nature and is based on the following two elementary observations: (i) given a Markov asset price process it is straightforward to construct a continuoustime Markov chain model whose law is close to that of , by approximating the generator of the process with an intensity matrix; (ii) the corresponding firstpassage problem for a continuoustime Markov chain can be solved explicitly via a closedform formula that only involves the generator matrix of the chain . More precisely, for a given Markov asset price model on the statespace with corresponding generator the algorithm for the pricing of any barrier product (including rebate options, which depend on the position at the moment of firstpassage) consists of the following two steps:

Construct a finite statespace and a generator matrix for the chain that approximates the operator on .
The form of the generators of Markov processes that commonly arise in pricing theory (including the local Lévy class) is well known from general theory (and are reviewed in Section 2). The statespace in step (i) is taken to be nonuniform with a higher density of points in the relevant areas such as spot and barrier levels. Subsequently, the generator matrix is defined by matching the instantaneous moments of the Markov processes and the chain on the statespace . This criterion implies in particular that the chain has the same average drift as the asset price process . Step (ii) of the algorithm consists of the evaluation of the closedform formulas for the firstpassage probabilities that can be derived employing continuoustime Markov chain theory (see Theorem 1 below). The evaluation of this formula consists of exponentiation of either the matrix or , which can be performed using the Padé approximation algorithm that is implemented in standard packages such as Matlab (see also [30]). The outputs of the algorithm yield arbitragefree prices in a certain continuous time Markov chain model for the risky asset price process. This feature is a consequence of the fact that the algorithm is based on the construction of an approximate stochastic model. This property is not shared by many other numerical methods used in practice that are based on purely analytical considerations.
We implemented this algorithm for a number of models that include the features of local volatility as well as jumps. We obtained an accurate match with the numerical results under diffusion and Lévy models that were considered elsewhere in the literature (see Section 6). We prove that by refining the grid the prices generated by this approach converge to those of the limiting model. We also establish, under additional regularity assumptions, an error bound that is linear in the spatial mesh size and the truncation error (see Section 5). We showed that an additional logarithmic factor may arise in this error bound when the Lévy density has a pole of order two at the origin. Numerical experiments (reported in Section 6) appear to suggest that, for a number of models, the error actually decays quadratically in the spatial mesh size. An extension to the case of timedependent characteristics is presented in an unabridged version of this paper [49], which also includes extra numerical examples, as well as some sample code.
There is good deal of literature devoted to the study of the (weak) convergence of Markov chains to limiting processes. However the (sharp) rates of convergence of prices generated by the Markov chain approximation to those of the limiting model are rarely available, especially for barrier options. For discrete time Markov chains some explicit rates have been established (see e.g. Broadie et al. [7] and Gobet and Menozzi [25]). Establishing the sharp rates of convergence for specific models remains an open question, left for future research.
The remainder of the paper is organized as follows. In Section 2 we define the class of models and barrier option contracts that is considered, and state some preliminary results about Markov processes. Section 3 presents the formulas for the firstpassage quantities of the continuoustime Markov chains. In Section 4 we describe the discretization algorithm to construct the intensity matrix of the chain . Section 5 states the convergence results, which are proved in Appendix A. Numerical results are presented in Section 6 and Section 7 concludes the paper.
2. Problem setting: Barrier options for Markov processes
The problem under consideration is that of the valuation of general barrier options, which can be formulated as follows. Given a random process modelling the price evolution of a risky asset, nonnegative payoff and rebate functions and , and a set specifying the range of values for which the contract ‘knocks out’, it is of interest to evaluate the expected discounted value of the random cash flow associated to a general barrier option contract
(2.1) 
where denotes the indicator of a set and
is the first time that enters the set . Furthermore, it is relevant to quantify the sensitivities of this value with respect to different parameters such as the spot value . The cash flow in (2.1) consists of a payment in the case the contract has not knocked out by the time , and a rebate if it has. Examples of commonly traded options included in this setting are the downandout, upandout and double knockout options. In particular, by taking we retrieve the case of a standard European claim with payoff at maturity .
We will consider this valuation problem in a Markovian setting, assuming that the underlying is a Markov process with statespace defined on some filtered probability space , where denotes the standard filtration generated by . Thus, takes values in and satisfies the Markov property:
(2.2) 
for all and bounded Borel functions , where denotes the expectation under the probability measure and is given by
(2.3) 
By taking expectations in (2.2) we see that the family forms a semigroup:
Informally, these conditions state that the expected value of the random cash flow occurring at time conditional on the available information up to time depends on the past via the value only. Setting the rate of discounting equal to a nonnegative constant , for any pair of nonnegative Borel functions and the expected discounted value of the barrier cash flow (2.1) at the epoch , the earlier of maturity and the first entrance time , is given by
(2.4) 
If represents the price of a tradeable asset, is the riskfree rate, is the dividend yield and the process is a martingale, standard arbitrage arguments imply that no arbitrage is introduced if expression (2.4) is used as the current price of the option with payoff (2.1).
Before proceeding we review some key concepts of the standard Markovian setup that will be needed in the sequel. For background on the (general) theory of Markov processes we refer to the classical works Chung and Walsh [11], Ethier and Kurtz [20], Itô and McKean [33] and Rogers and Williams [52] (the latter two in particular treat the case of Markov processes with continuous sample paths). In what follows we will restrict to be in a subclass of Markov processes for which, if the function is continuous and tends to zero at infinity, the expected payoff has the following properties: it depends continuously on the spot and on expiry and also decays to zero when tends to infinity. More precisely, denoting by the set of continuous functions on that tend to zero at infinity, we make the following assumption:
Assumption 1.
is a Feller process on , that is, for any , the family , with defined in (2.3), satisfies:

for any ;

for any .
The Feller property guarantees that there exists a version of the process with càdlàg paths satisfying the strong Markov property. In particular, a Feller process is a Hunt process.
Throughout the paper we will take the knockout set to be of the form
(2.5) 
which includes the cases of double and single barrier options—the latter by taking or To rule out degeneracies we will make the following assumption on the behaviour of at the boundary points and :
Assumption 2.
For every we have where .
This assumption states that the first entrance times into and its interior coincide almost surely, if the spot is equal to . If , a sufficient condition for Assumption 2 to be satisfied is for ; that is, when started at or , the process immediately enters .
The class of Feller processes satisfying Assumption 2 includes many of the models employed in quantitative finance such as (Feller)diffusions, jumpdiffusions with nongenerate diffusion coefficient and Lévy processes whose Lévy measure admits a density.
The family is determined by its infinitesimal generator that is defined as
(2.6) 
for any function
for which the righthand side of (2.6)
converges in the strong sense.
We next give a few examples of Feller processes with their generators.
Example 1.
A diffusion asset price model evolves under a riskneutral measure according to the stochastic differential equation (SDE)
(2.7) 
where
is the initial price, and
is a given measurable function. To guarantee the absence of arbitrage
we assume that is chosen such that
the discounted process is a martingale,
If, in addition, infinity is not entrance
(2.8) 
where denotes the derivative of with respect to (see [20, Sec. 8.1]).
Example 2.
The price process in an exponential Lévy model given by
(2.9) 
where and are constants representing the interest rate and dividend yield and is a Lévy process, such that for all . By construction, is a martingale. Further, if and only if the Lévy measure integrates at infinity, that is,
(2.10) 
The law of is determined by its characteristic exponent , which is related to the characteristic function of by and which, under condition (2.10), has the LévyKhintchine representation
where is the characteristic triplet, with and the Lévy measure, which satisfies the integrability condition . The process is a Feller process with an infinitesimal generator acting on as (cf. Sato [53, Thm. 31.5])
where
Example 3.
More generally, one may specify the price process by directly prescribing its generator to act on sufficiently regular functions as
(2.11) 
where is given in (2.8) and
where for every , is a (Lévy) measure with support in such that
The discounted process is a local martingale. Sufficient conditions on and to guarantee the existence of a Feller process corresponding to this generator were established in Kolokoltsov [36, Thm. 1.1].
A key step in solving the barrier valuation problems is the observation that the expected values of knockout options and general barrier options can be expressed in terms of the marginal distributions of two Markov processes associated to . Given a Markov process the processes that have the same dynamics as before entering , but are stopped or jump to the graveyard state upon entering the set , respectively, are themselves Markov processes. Let denote the killed process and the process killed at rate that is stopped upon entering the set . Then the value of the barrier options can be expressed as
where we assume that and the function is defined as . To calculate the valuefunctions of barrier options written on the underlying price process we thus need to identify and . This can be achieved by employing the infinitesimal generators of the semigroups and associated to the Markov processes and which are explicitly expressed in terms of the generator as follows:
Lemma 1.
(i) For any , where is the domain of the generator , we have
(2.12) 
where the convergence is pointwise. If is a Feller process, and
(2.13) 
then , where is the infinitesimal generator of the semigroup .
(ii) Let and assume that satisfies
(2.14) 
Then
(2.15) 
where the convergence is pointwise. If is a Feller process and
then , where is the infinitesimal generator of the semigroup .
Lemma 1 is a straightforward consequence of the definition of the infinitesimal generator and the HilleYosida theorem, see e.g. [53, Lemma 31.7] (see [49] for the complete proof of Lemma 1). If and are themselves Feller processes, the relations between and , and between and can formally be expressed as follows:
(2.16) 
Equation (2.16) can be given a precise meaning if, for example, and can be defined as a selfadjoint operator on a separable Hilbert space (see e.g. Ch. XII in Dunford and Schwarz [17], or Hille and Philips [31]). By determining the spectral decompositions of and one can construct spectral expansions of and , which in the case of a discrete spectrum reduces to a series expansion. See Linetsky [42, 43, 44] for a development of this spectral expansion approach for onedimensional diffusion models in finance, and an overview of related literature.
When (asymmetric) jumps are present, the operator is nonlocal and not selfadjoint, and the spectral theory has been less well developed, with fewer explicit results. Here we will follow a different approach: we will approximate by a finite state Markov chain, and show that for the approximating chain a matrix analog of the identities (2.12)–(2.16) holds true, where the infinitesimal generators and can be easily obtained from . We give a selfcontained development of this approach in Section 3.
3. Exit probabilities for continuoustime Markov chains
Given a Markov price process of interest, the idea is to construct a continuoustime finite state Markov chain whose dynamics are “close” to those of , and to calculate the relevant expectations for this approximating chain. In this section we will focus on the latter; we will return to the question of how to construct such a chain in Section 4. Assume therefore we are given a finite state continuoustime Markov chain . From Markov chain theory it is well known that the chain is completely specified by its statespace (or grid) and its generator matrix , which is an square matrix with zero row sums and nonpositive diagonal elements, if has elements. Given the generator matrix the family of transition matrices of , defined by for , is given by
In particular, the expexted discounted payoff at maturity is then given by
(3.1) 
for and any function . Here and throughout the paper we will identify any square matrix and any vector in with functions
where the vectors denote the corresponding standard basis vectors of and stands for transposition.
The generator can be retrieved from the family of matrices defined above by differentiation at , that is,
In order for itself to form a pricing model, defined under a martingale measure, we will suppose in addition that is nonnegative and that the discounted process is a martingale; more precisely we assume that is a martingale where
(3.2) 
with the hitting time of the “boundary” which consists of the smallest and largest elements of the state space , and . The Markov property of implies that the process given in (3.2) is a martingale precisely if
(3.3) 
where the function is given by . Below we will show how to express the exit probabilities of the chain using matrix exponentiation in a way that is identical in form to the expected value (3.1) of a European payoff. To that end, we partition into a ‘continuation’ set and a ‘knockout’ set , where
(3.4) 
and define the first exit time of from by
(3.5) 
where we use the convention and where we will take the set as in (2.5).
The value of a general barrier knockout option with a rebate depends on the joint distribution of the exit time from and the positions of the underlying at maturity and at the moment of exit. The corresponding quantities for the chain can be expressed in terms of two transformations of , namely the chain that is killed upon exiting and the chain that is absorbed at that instance, respectively. Correspondingly, we associate to the generator matrix , two matrices: the matrix , where , and the matrix , defined by
(3.6)  
(3.7) 
We can now state the key result of this section:
Theorem 1.
For any , and and any function it holds that
(3.8) 
In particular, for and with for we have that
(3.9)  
(3.10)  
(3.11) 
where , the restriction of to .
Formulas (3.9)–(3.11) give a simple way of computing barrier option prices by a single matrix exponentiation. The expectation in (3.9) can be obtained by computing the spectral decomposition of the matrix and applying the formula . The powerful Padé approximation method for matrix exponentiation, described in [30], can also be used to compute efficiently the matrix exponentials in Theorem 1. Since the statespace is finite, Theorem 1 is a corollary of Lemma 1. We present next a direct probabilistic derivation.
Proof. To prove equation (3.8), we will verify that the expected value of an ArrowDebreu barrier security that pays 1 precisely if is in the state at the earlier of the maturity and the knockout time is given by
(3.12) 
For a given time grid with denote by the expected value of the corresponding discretely monitored ArrowDebreu security and let
be the corresponding time at which the barrier is crossed. Since the paths of the chain are piecewise constant, it follows that and as tends to infinity. Hence the expected values of the discretely monitored ArrowDebreu securities converge to the expected value of the continuously monitored one,
Clearly, since is the knockout set, it holds for all that
where is a square matrix of size with if and zero else. Further, for an application of the Markov property of shows that
Combining the two cases, iterating the argument and using the differentiability of at shows that
since . When tends to zero, this expression converges to , which completes the proof of (3.12) and hence implies (3.8). Equation (3.11) follows then directly by applying (3.8) to . Noting that for any and any we have
we get that
which yields (3.9). Finally, (3.10) is a direct consequence of (3.1), (3.9) and the fact that a European option with payoff is equal to the sum of a knockout barrier option and a knockin barrier option with the same payoff and same knockout/knockin levels.
4. Construction of the Markov chain
The Markov chain approximation algorithm for homogeneous Feller processes can now be described as follows:

Construct an approximating Markov chain:

specify a (nonuniform) grid ;

define a generator matrix of a Markov chain with statespace .

A suitable choice of the grid in step (1a) is essential for the effectiveness of the above pricing algorithm. The construction of an optimal grid (according to some criterion) is a topic of separate study, which will not be pursued further in this paper. One of the features of a good grid is that it has sufficient resolution in regions of interest, such as the current spot value and the barrier levels, which is a necessary condition for constructing a Markov chain market model that approximates well the dynamics of the given price process. Another desirable feature is that the grid “covers” a sufficiently large part of the state space, which is needed to control the truncation error that arises when approximating an infinite state space by a finite state space. To employ a uniform grid that satisfies these conditions would be computationally expensive. The use of adaptive meshes for option pricing was proposed in Figlewski and Gao [21]. Here we employ the following procedure for generating a suitable nonuniform grid , based on an algorithm from [55]:

Pick and the density parameters , , and the smallest and largest values of the grid , where .

Define GenerateSubGrid, for , where .

Set ,
where the subgrid is generated by the following procedure:
GenerateSubGrid

Compute , .

Define the lower part of the grid by the formula for . Note that .

Define the upper part of the grid using the formula for . Note that .
Return
The nonuniform statespace for the Markov chain is constructed by concatenating the three subgrids that are generated by specifying lower, middle and upper points and density parameters and . A smaller density corresponds to a grid that is more concentrated around the middle point . Observe that the algorithm above places the current spot and the barrier levels on the grid and that the resolution of the grid around and can be controlled by the density parameters ,. A Matlab implementation of this grid generator can be downloaded from [48]. The remainder of this section will be devoted to step (1b) in the algorithm described above.
4.1. Jump processes with statedependent characteristics
The construction of the generator matrix of the approximating Markov chain is now carried out in two steps: we first define the jump matrix , which corresponds to the discretization of the jump measure , and then characterize a tridiagonal generator matrix by stipulating that the Markov chain with the generator has the same instantaneous moments as the process .
We start by building the statespace with elements using the algorithm described in the beginnig of this section. Define the sets
where the “boundary” consist of the smallest (i.e. ) and largest (i.e. ) elements in and the “interior” is the complement of the boundary. For any given we associate to the set defined by
The set consists of the relative jump sizes of jumps starting from and arriving at any other point in . The th row of the jump part of the generator of is obtained by discretizing the jump measure on the set . In particular let and define a function such that
(4.1) 
where . A possible natural choice for would be the midpoint of the interval
We can now define the jump part of the generator as
(4.2) 
where and
(4.3) 
Note that the function generates a partititon of . The jump intensities of the chain, defined in formula (4.2), are obtaiend by integrating the Lévy measure over the corresponding part of the partition. For we set for all . It is clear that the matrix constructed in this way is a generator matrix.
In the second step we match the first and second instantaneous moments of the asset price process . In other words the chain must satisfy conditions
(4.4) 
for all starting states . Note that condition (4.4) implicitly assumes that the second instantaneous moment of exists. This is the case if the jump measure satisfies the following condition
(4.5) 
which we now assume to hold.
The task now is to find a tridiagonal generator matrix such that the chain generated by the sum satisfes (4.4). The tridiagonal matrix therefore has to satisfy the following conditions
(4.6)  
(4.7)  
(4.8) 