Continuously monitored barrier options under Markov processes

# Continuously monitored barrier options under Markov processes

## Abstract.

In this paper we present an algorithm for pricing barrier options in one-dimensional Markov models. The approach rests on the construction of an approximating continuous-time 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 jump-diffusion. We also provide a convergence proof and error estimates for this algorithm.

Acknowledgements: We would like to thank the anonymous referees, Bjorn Eriksson, Mike Giles, Vassili Kolokoltsov, Steven Kou, Sergei Levendorskii, Dilip Madan, Jan Obloj and Johannes Stolte as well as the participants of the 2009 Leicester workshop on Spectral and Cubature Methods in Finance and Econometrics for useful suggestions and constructive comments, which led to improvements of the paper. Research supported by EPSRC grant EP/D039053. This research was carried out while AM was based at Imperial College London.

## 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 double-no-touch options and knock-in or knock-out 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 Black-Scholes 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 non-parametric 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 risk-neutral 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 one-dimensional 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 log-asset price depends in a non-linear 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 jump-processes 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 jump-diffusion 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 first-passage 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 barrier-type options. It is well known that in this case a straightforward Monte Carlo simulation algorithm will be time-consuming and yield unstable results for the prices and especially the sensitivities. The knock-in/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:

1. spectral expansions for several parametric diffusion models (Davydov and Linetsky [15], Lipton [45]),

2. transform based approaches for exponential Lévy models (Boyarchenko and Levendorskii [5], Geman and Yor [24], Jeannin and Pistorius [35], Kou and Wang [38], Sepp [54]),

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 so-called Wiener-Hopf 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 Black-Scholes 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 Black-Scholes 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 continuous-time 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 one-dimensional Markov processes, which in particular includes the case of local Lévy models as well as local volatility jump-diffusions 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 continuous-time Markov chain model whose law is close to that of , by approximating the generator of the process with an intensity matrix; (ii) the corresponding first-passage problem for a continuous-time Markov chain can be solved explicitly via a closed-form formula that only involves the generator matrix of the chain . More precisely, for a given Markov asset price model on the state-space with corresponding generator the algorithm for the pricing of any barrier product (including rebate options, which depend on the position at the moment of first-passage) consists of the following two steps:

1. Construct a finite state-space and a generator matrix for the chain that approximates the operator on .

2. To value knock-out and rebate options, obtain the matrices and by the procedure in Figure 1 and apply closed-form formulas in terms of these matrices (given in equations (3.9) and (3.11) below).

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 state-space in step (i) is taken to be non-uniform 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 state-space . 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 closed-form formulas for the first-passage probabilities that can be derived employing continuous-time 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 arbitrage-free 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 time-dependent 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 first-passage quantities of the continuous-time 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, non-negative 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) g(ST)I{τA>T}+h(SτA)I{τA≤T},

where denotes the indicator of a set and

 τA=inf{t≥0:St∈A}

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 down-and-out, up-and-out and double knock-out 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 state-space defined on some filtered probability space , where denotes the standard filtration generated by . Thus, takes values in and satisfies the Markov property:

 (2.2) E[f(St+s)|Ft]=Psf(St),

for all and bounded Borel functions , where denotes the expectation under the probability measure and is given by

 (2.3) Psf(x):=Ex[f(Ss)]:=E[f(Ss)|S0=x].

By taking expectations in (2.2) we see that the family forms a semigroup:

 Pt+sf=Pt(Psf),for all s,t≥0, and P0f=f.

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 non-negative constant , for any pair of non-negative 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) Ex[e−rTg(ST)I{τA>T}]+Ex[e−rτAh(SτA)I{τA≤T}].

If represents the price of a tradeable asset, is the risk-free 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 knock-out set to be of the form

 (2.5) A=[0,ℓ]∪[u,∞),0≤ℓ

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, jump-diffusions with non-generate 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) Lf(x):=limt↓01t(Ptf−f)(x)

for any function for which the right-hand side of (2.6) converges in the strong sense.1 The set of such functions is called the domain of the operator and is dense in . These fundamental facts about semigroups and their generators can be found in [20, Ch. 1].

We next give a few examples of Feller processes with their generators.

###### Example 1.

A diffusion asset price model evolves under a risk-neutral measure according to the stochastic differential equation (SDE)

 (2.7) \rm dStSt=γ\rm dt+σ(St)\rm dWt,

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 entrance2 for and is a continuous function, then is a Feller process, and its infinitesimal generator acts on 3 as

 (2.8) LDf(x)=σ(x)2x22f′′(x)+γxf′(x),

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) St:=S0e(r−d)teLtE[eLt]

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) ∫(1,∞)eyν(\rm dy)<∞.

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évy-Khintchine representation

 Ψ(s)=ics−σ2s22+∫R(eisy−1−isy)ν(\rm dy),

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])

 Lf(x)=σ2x22f′′(x)+γxf′(x)+∫R[f(xey)−f(x)−xf′(x)(ey−1)]ν(\rm dy),

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) Lf(x)=LDf(x)+LJf(x),

where is given in (2.8) and

 LJf(x)=∫(−1,∞)[f(x(1+y))−f(x)−f′(x)xy]ν(x,\rm dy),

where for every , is a (Lévy) measure with support in such that

 ∫(−1,∞)min{y2,|y|}ν(x,\rm dy)<∞.

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 knock-out 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

 Ex[g(ST)I{τA>T}]=Ex[g(ˆSAT)] =: ˆPATg(x), Ex[e−rTg(ST)I{τA>T}]+Ex[e−rτAh(SτA)I{τA≤T}] = Ex[e−r(T∧τA)f(ST∧τA)] = Ex[f(˜SAT)]=:˜PATf(x),

where we assume that and the function is defined as . To calculate the value-functions 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) limt↓0t−1(˜PAtf(x)−f(x)) = ˜k(x):={0,x∈A,(L−r)f(x),x∉A,

where the convergence is pointwise. If is a Feller process, and

 (2.13) limx→∂A{Lf(x)−rf(x)}=0,where % ∂A is the boundary of A,

then , where is the infinitesimal generator of the semigroup .

(ii) Let and assume that satisfies

 (2.14) Px(τA≤t)=o(t) as t↘0org(SτA)=0Px-a.s.

Then

 (2.15) limt↓0t−1(ˆPAtg(x)−g(x)) = Lg(x),

where the convergence is pointwise. If is a Feller process and

 g|A=0andlimx→∂ALg(x)=0,

then , where is the infinitesimal generator of the semigroup .

Lemma 1 is a straightforward consequence of the definition of the infinitesimal generator and the Hille-Yosida 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) ˜PAt=exp(t˜LA),ˆPAt=exp(tˆLA).

Equation (2.16) can be given a precise meaning if, for example, and can be defined as a self-adjoint 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 one-dimensional diffusion models in finance, and an overview of related literature.

When (asymmetric) jumps are present, the operator is non-local and not self-adjoint, 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 self-contained development of this approach in Section 3.

## 3. Exit probabilities for continuous-time Markov chains

Given a Markov price process of interest, the idea is to construct a continuous-time 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 continuous-time Markov chain . From Markov chain theory it is well known that the chain is completely specified by its state-space (or grid) and its generator matrix , which is an square matrix with zero row sums and non-positive diagonal elements, if has elements. Given the generator matrix the family of transition matrices of , defined by for , is given by

 Pt=exp(tΛ).

In particular, the expexted discounted pay-off at maturity is then given by

 (3.1) Ex[e−rTϕ(XT)] = e−rT⋅(exp(TΛ)ϕ)(x)

for and any function . Here and throughout the paper we will identify any square matrix and any vector in with functions

 A:G×G→R, A(x,y):=e′xAey,x,y∈G,and ϕ:G→R, ϕ(x):=e′xϕ,x∈G,

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,

 Pt=I+tΛ+o(t) as t↘0.

In order for itself to form a pricing model, defined under a martingale measure, we will suppose in addition that is non-negative and that the discounted process is a martingale; more precisely we assume that is a martingale where

 (3.2) Mt={e−γtXt,t<ζ,e−γζXζ,t≥ζ,

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) (Λη)(x)=γη(x)for x∈G∖∂G,

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 pay-off. To that end, we partition into a ‘continuation’ set and a ‘knock-out’ set , where

 (3.4) ˆG := {x∈G:x∈Ac},

and define the first exit time of from by

 (3.5) τ:=inf{t∈R+:Xt∉ˆG},

where we use the convention and where we will take the set as in (2.5).

The value of a general barrier knock-out 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) ˜Λr(x,y) := ⎧⎪ ⎪⎨⎪ ⎪⎩Λ(x,y)−rif x∈ˆG, % x=y,Λ(x,y)if x∈ˆG, y∈G, x≠y,% 0if x∈ˆGc, y∈G, (3.7) ˆΛ(x,y) := ˆΛ0(x,y):=Λ(x,y)if x∈ˆG, y∈ˆG.

We can now state the key result of this section:

###### Theorem 1.

For any , and and any function it holds that

 (3.8) Ex[e−r(T∧τ)ϕ(XT∧τ)]=(exp(T˜Λr)ϕ)(x).

In particular, for and with for we have that

 (3.9) Ex[ψ(XT)I{τ>T}] = (exp(TˆΛ)ψ)(x)for anyx∈ˆG, (3.10) Ex[ϕ(XT)I{τ≤T}] = (exp(ΛT)ϕ)(x)−(exp(TˆΛ)ˆϕ)(x)IˆG(x)  for anyx∈G, (3.11) Ex[e−rτξ(Xτ)I{τ≤T}] = (exp(T˜Λr)ξ)(x)for anyx∈G,

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 state-space 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 Arrow-Debreu barrier security that pays 1 precisely if is in the state at the earlier of the maturity and the knock-out time is given by

 (3.12) Ex[e−r(T∧τ)I{XT∧τ=y}]=(exp(T˜Λr))(x,y)for all x,y∈G.

For a given time grid with denote by the expected value of the corresponding discretely monitored Arrow-Debreu security and let

 τn=inf{s∈Tn:Xs∉ˆG}

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 Arrow-Debreu securities converge to the expected value of the continuously monitored one,

 ˜P(n)T(x,y)=Ex[e−r(T∧τn)I{XT∧τn=y}]⟶Ex[e−r(T∧τ)I{XT∧τ=y}].

Clearly, since is the knock-out set, it holds for all that

 ˜P(n)t(x,y) = {1if x∈ˆGc, x=y0if x∈ˆGc, x≠y = (I−¯¯¯I)(x,y)for all x∈ˆGc, y∈G,

where is a square matrix of size with if and zero else. Further, for an application of the Markov property of shows that

 ˜P(n)T(x,y) = e−rΔt∑z∈GPΔt(x,z)˜P(n)T−Δt(z,y) = (¯¯¯I(e−rΔtPΔt)˜P(n)T−Δt)(x,y).

Combining the two cases, iterating the argument and using the differentiability of at shows that

 ˜P(n)T(x,y) = ((I−¯¯¯I+¯¯¯Ie−rΔtPΔt)˜P(n)T−Δt)(x,y) = ((I−¯¯¯I+¯¯¯Ie−rΔtPΔt)T/Δt)(x,y) = ((I+¯¯¯Ie−rΔt(PΔt−I)+¯¯¯I(e−rΔt−1))T/Δt)(x,y) = ((I+Δt(˜Λ0−r¯¯¯I)+o(Δt))T/Δt)(x,y),

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

 (ˆΛψ)(x)=(˜Λ0ψ0)(x),where ψ0:G→R is given by ψ0(y):=ψ(y)IˆG(y),

we get that

 (exp(TˆΛ)ψ)(x)=(exp(T˜Λ0)ψ0)(x),

which yields (3.9). Finally, (3.10) is a direct consequence of (3.1), (3.9) and the fact that a European option with pay-off is equal to the sum of a knock-out barrier option and a knock-in barrier option with the same pay-off and same knock-out/knock-in levels.

## 4. Construction of the Markov chain

The Markov chain approximation algorithm for homogeneous Feller processes can now be described as follows:

1. Construct an approximating Markov chain:

1. specify a (non-uniform) grid ;

2. define a generator matrix of a Markov chain with state-space .

2. Compute the barrier option prices using formulae (3.8)–(3.11).

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 non-uniform grid , based on an algorithm from [55]:

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

2. Define GenerateSubGrid, for , where .

3. Set ,

where the subgrid is generated by the following procedure:

GenerateSubGrid

1. Compute , .

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

3. Define the upper part of the grid using the formula for . Note that .

Return

The non-uniform state-space 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 state-dependent 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 tri-diagonal generator matrix by stipulating that the Markov chain with the generator has the same instantaneous moments as the process .

We start by building the state-space with elements using the algorithm described in the beginnig of this section. Define the sets

 ∂G:={x1,xN}andGo:=G∖∂G,

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

 Gx := {zx−1:z∈G}.

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) αx(yi)∈(yi,yi+1)fori∈{1,…,N−1}

where . A possible natural choice for would be the mid-point of the interval

 αx(yi)=yi+yi+12.

We can now define the jump part of the generator as

 (4.2) ΛJ(x,x(1+yi)) := ∫αx(yi)αx(yi−1)ν(x,% \rm dy)

where and

 (4.3) ΛJ(x,x) := −∑z∈G∖{x}ΛJ(x,z).

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) Ex[(SΔt−S0)j] = Ex[(XΔt−X0)j]+o(Δt),forx∈Go,j∈{1,2}

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) ∫(−1,∞)y2ν(x,\rm dy)<∞for % allx∈E

which we now assume to hold.

The task now is to find a tri-diagonal generator matrix such that the chain generated by the sum satisfes (4.4). The tri-diagonal matrix therefore has to satisfy the following conditions

 (4.6) ∑z∈GΛD(x,z) = 0andΛD(x,z)+ΛJ(x,z)≥0∀z∈G∖{x}, (4.7) ∑z∈GΛD(x,z)(z−x) = (r−d)x−∑z′∈GΛJ(x,z′)(z′−x), (4.8) ∑z∈GΛD(x,z)(z−x)2 = x2[σ(x)2+∫∞−1y2ν(x,\rm dy)]−∑z′∈GΛJ