# Valuation of barrier options using sequential Monte Carlo

## Pavel V Shevchenko and Pierre Del Moral

#### Need to know

• This paper presents Sequential Monte Carlo (SMC) method for pricing barrier options.
• The authors demonstrate that SMC is more efficient then standard Monte Carlo for pricing barrier options.
• General formulas and references are provided for pricing other exotic options using SMC.

#### Abstract

ABSTRACT

Sequential Monte Carlo (SMC) methods have been used successfully in many applications in engineering, statistics and physics. However, they are seldom used in financial option pricing literature and its practice. We present an SMC method for pricing barrier options with continuous and discrete monitoring of the barrier condition. With our method, simulated asset values rejected due to the barrier condition are resampled from asset samples that do not breach the barrier condition, improving the efficiency of the option price estimator, while with the standard Monte Carlo method many simulated asset paths can be rejected by the barrier condition, making it harder to estimate the option price accurately. We compare the SMC with the standard Monte Carlo method and demonstrate that there is little extra effort required to implement the former compared with the latter, while the improvement in price estimation can be significant. Both methods result in unbiased estimators for the price converging to the true value as 1 / √M , where M is the number of simulations (asset paths). However, the variance of the SMC estimator is smaller and does not grow with the number of time steps, unlike standard Monte Carlo. In this paper, we demonstrate that the SMC can successfully be used for pricing barrier options. SMC methods can also be used for pricing other exotic options and for cases with many underlying assets and additional stochastic factors, such as stochastic volatility; we provide general formulas and references.

#### Introduction

Sequential Monte Carlo (SMC) methods (also referred to as particle methods) have been used successfully in many applications in engineering, statistics and physics for many years, especially in signal processing, state-space modeling and estimating rare event probability. The SMC method coincides with the quantum Monte Carlo method introduced as an heuristic-type scheme in physics by Enrico Fermi in 1948 while studying neutron diffusions. From a mathematical point of view, SMC methods can be seen as mean-field particle interpretations of Feynman–Kac models. (For a detailed analysis of these stochastic models and applications, we refer the reader to Del Moral (2004, 2013) and the references therein.) The application of these particle methods in mathematical finance is relatively recent. For instance, using the rare event interpretation of SMC, Carmona et al (2009) and Del Moral and Patras (2011) proposed an SMC algorithm for computing the probabilities of simultaneous defaults in large credit portfolios, and Targino et al (2015) developed the SMC for capital allocation problems. There are many articles that use the SMC for estimation of stochastic volatility, jump diffusion and state-space price models (see, for example, Johannes et al 2009; Peters et al 2013). The applications of SMC methods in option pricing have been studied recently by the second author in a series of articles (Carmona et al 2012; Del Moral et al 2011, 2012a, 2012b; Jasra and Del Moral 2011). However, these methods are not widely known by option pricing practitioners or in the option pricing literature.

The aim of this paper is to provide a simple explanation of the SMC method and illustrate its efficiency. For simplicity, we consider barrier options with a simple geometric Brownian motion for the underlying asset. The SMC can also be used for pricing other exotic options and different underlying stochastic processes; we provide general formulas and references.

Barrier options are widely used in trading. The option is extinguished (knocked out) or activated (knocked in) when an underlying asset reaches a specified level (barrier). Many more complex related instruments such as bivariate barrier, ladder, step-up or step-down barrier options have become very popular in over-the-counter markets. In general, these options can be considered as having a payoff dependent upon the path extremums of the underlying assets. A variety of closed-form solutions for such instruments on a single underlying asset have been obtained in the classical Black–Scholes settings of constant volatility, interest rate and barrier level (see, for example, Heynen and Kat 1994b; Kunitomo and Ikeda 1992; Rubinstein and Reiner 1991). If the barrier option is based on two assets, then a practical analytical solution can be obtained for some special cases considered in Heynen and Kat (1994a) and He et al (1998).

In practice, however, numerical methods are used to price barrier options for a number of reasons, eg, if the assumptions of constant volatility and drift are relaxed or the payoff is too complicated. Numerical schemes such as binomial and trinomial lattices (Hull and White 1993; Kat and Verdonk 1995) or finite-difference schemes (Dewynne and Wilmott 1994) can be applied to the problem. However, the implementation of these methods can be difficult. Also, if more than two underlying assets are involved in the pricing equation, then these methods are not practical.

The Monte Carlo simulation method is a good general pricing tool for such instruments; for a review of advanced Monte Carlo methods for barrier options, see Gobet (2009). Many studies address finding the extremums of the continuously monitored assets by sampling assets at discrete dates. The standard discrete-time Monte Carlo approach is computationally expensive, as a large number of sampling dates and simulations are required. Loss of information about all parts of the continuous-time path between sampling dates introduces a substantial bias for the option price. The bias decreases very slowly as $1/\smash{\sqrt{N}}$ for $N\gg 1$, where $N$ is the number of equally spaced sampling dates (see Broadie et al (1997), which also shows how to calculate the discretely monitored barrier option approximately via a continuous barrier case with some shift applied to the barrier). Also, extrapolation of Monte Carlo estimates to the continuous limit is usually difficult due to finite sampling errors. For the case of a single underlying asset, it was shown by Andersen and Brotherton-Racliffe (2006) and Beaglehole et al (1997) that the bias can be eliminated by a simple conditioning technique, the so-called Brownian bridge simulation. This method is based on the simulation of a one-dimensional Brownian bridge extremum between the sampled dates according to a simple analytical formula for the distribution of the extremum (or just multiplying the simulated option payoff by the conditional probability of the path not crossing the barrier between the sampled dates); see also Glasserman (2004), pp. 368–370. The technique is very efficient in the case of an underlying asset following the standard lognormal process, because only one time step is required to simulate the asset path and its extremum if the barrier, drift and volatility are constant over time. A closely related method of sampling the underlying asset conditional on not crossing a barrier is studied in Glasserman and Staum (2001). The method of Brownian bridge simulation can also be applied in the case of multiple underlying assets, as studied by Shevchenko (2003). Importance sampling and the control variates method can be applied to reduce the variance of the barrier option price Monte Carlo estimator; for a textbook treatment, see Glasserman (2004). To improve time discretization scheme convergence in the case of more general underlying stochastic processes, Giles (2008a, 2008b) introduced a multilevel Monte Carlo path simulation method for the pricing of financial options, including barrier options, that improves the computational efficiency of Monte Carlo path simulation by combining results using different numbers of time steps. Gobet and Menozzi (2010) developed a procedure for multidimensional stopped diffusion processes, accounting for boundary correction by shifting the boundary, which can be used to improve the barrier option Monte Carlo estimates in the case of multi-asset and multi-barrier options with more general underlying processes.

However, the coefficient of variation of the Monte Carlo estimator grows when the number of asset paths rejected by the barrier condition increases (ie, the probability of the asset path reaching maturity without breaching the barrier decreases, for example, when barriers get closer to the asset spot). This can be improved by the SMC method, which resamples asset values rejected by the barrier condition from the asset samples that do not breach the barrier condition at each barrier monitoring date. Both SMC and Monte Carlo estimators are unbiased and converge to the true value as $1/\smash{\sqrt{M}}$, where $M$ is the number of simulations (asset paths), but the SMC has smaller variance.

In this paper we present an SMC algorithm and compare SMC and Monte Carlo estimators. We focus on the case of one underlying asset for easy illustration, but the algorithm can easily be adapted for the case with many underlying assets and with additional stochastic factors, such as stochastic volatility. Note that we do not address the error due to time discretization, but instead improve the accuracy of the option price sampling estimator for a given time discretization.

The organization of the paper is as follows. Section 2 describes the model and notation. In Section 3 we provide the basic formulas for Feynman–Kac representation underlying the SMC method. Section 4 presents SMC and Monte Carlo algorithms and corresponding option price estimators. The use of importance sampling to improve SMC estimators is discussed in Section 5. Numerical examples are presented in Section 6. Concluding remarks are given Section 7.

## 2 Model

Assume that underlying asset $S_{t}$ follows the risk-neutral process

 $\mathrm{d}S_{t}=S_{t}\mu\,\mathrm{d}t+S_{t}\sigma\,\mathrm{d}W_{t},$ (2.1)

where $\mu=r-q$ is the drift, $r$ is the risk-free interest rate, $q$ is the continuous dividend rate (corresponding to the foreign interest rate if $S_{t}$ is the exchange rate, or to continuous dividends if $S_{t}$ is the stock), $\sigma$ is the volatility and $W_{t}$ is the standard Brownian motion. The interest rate can be a function of time, and drift and volatility can be functions of time and the underlying asset. In this paper, we do not consider time discretization errors; for simplicity, we hereafter assume that model parameters are piecewise constant functions of time.

### 2.1 Pricing the barrier option

Today’s fair price of the continuously monitored knock-out barrier option with a lower barrier $L_{t}$ and an upper barrier $U_{t}$ can be calculated as the expectation with respect to the risk-neutral process (2.1): given information today at $t_{0}=0$ (ie, conditional on $S_{0}=s_{0}$),

 $Q_{\mathrm{C}}=B_{0,T}\mathbb{E}(h(S_{T})1_{\mathcal{A}_{t}}(S_{t})_{t\in[0,T]% }),\quad B_{0,T}=\exp\bigg({-}\int_{0}^{T}r(\tau)\,\mathrm{d}\tau\bigg),$ (2.2)

where $\smash{B_{0,T}}$ is the discounting factor from maturity $T$ to $t_{0}=0$; $\smash{1_{\mathcal{A}}(x)}$ is the indicator function, which equals 1 if $x\in\mathcal{A}$ and 0 otherwise; $h(x)$ is the payoff function, ie, $h(x)=\max(x-K,0)$ for a call option and $h(x)=\max(K-x,0)$ for a put option, where $K$ is the strike price; and $\mathcal{A}_{t}=(L_{t},U_{t})$. All standard barrier structures (lower barrier only, upper barrier only or several window barriers) can be obtained by setting $L_{t}=0$ or $U_{t}=\infty$ for the corresponding time periods.

Assume that drift, volatility and barriers are piecewise constant functions of time for the time discretization $0=t_{0}. Denote the corresponding asset values by $S_{0},S_{1},\dots,S_{N}$, the lower and upper barriers by $L_{1},\dots,L_{N}$ and $U_{1},\dots,U_{N}$, respectively, and the drift and volatility by $\mu_{1},\dots,\mu_{N}$ and $\sigma_{1},\dots,\sigma_{N}$, respectively. That is, $L_{1}$ is the lower barrier for time period $[t_{0},t_{1}]$; $L_{2}$ is the lower barrier for $[t_{1},t_{2}]$, and so on, and similarly for the upper barrier, drift and volatility. If there is no lower or upper barrier during $[t_{n-1},t_{n}]$, then we set $L_{n}=0$ and $U_{n}=\infty$, respectively.

We denote the transition density from $S_{n}$ to $S_{n+1}$ by $f(S_{n+1}\mid S_{n})$, which is just a lognormal density in the case of process (2.1), with solution

 $S_{n}=S_{n-1}\exp((\mu_{n}-\tfrac{1}{2}\sigma^{2}_{n})\delta t_{n}+\sigma_{n}% \sqrt{\delta t_{n}}Z_{n}),\quad n=1,\dots,N,$ (2.3)

where $\delta t_{n}=t_{n}-t_{n-1}$ and $Z_{1},\dots,Z_{N}$ are independent and identically distributed random variables from the standard normal distribution.

In the case of the barrier monitored at $t_{0},t_{1},\dots,t_{N}$ (discretely monitored barrier), the option price (2.2) simplifies to

 $Q_{\mathrm{D}}=B_{0,T}\mathbb{E}\bigg(h(S_{N})\prod_{n=1}^{N}1_{(L_{n},U_{n})}% (S_{n})\bigg).$ (2.4)

This is a biased estimate of the continuously monitored barrier option $Q_{\mathrm{C}}$ such that $Q_{\mathrm{D}}\rightarrow Q_{\mathrm{C}}$ for $\delta t_{n}\rightarrow 0$; see Broadie et al (1997), which also shows how to calculate the discretely monitored barrier option approximately via the continuous barrier option price, with some shift applied to the barrier in the case of one-dimensional Brownian motion (for the high-dimensional case and more general processes, see Gobet and Menozzi 2010).

In the case of a continuously monitored barrier, the barrier option price expectation (2.2) can be written as

 $Q_{\mathrm{C}}=B_{0,T}\int_{L_{1}}^{U_{1}}\mathrm{d}s_{1}\,f(s_{1}\mid s_{0})g% (s_{0},s_{1})\cdots\int_{L_{N}}^{U_{N}}\mathrm{d}s_{N}\,f(s_{N}\mid s_{N-1})g(% s_{N-1},s_{N})h(s_{N}),$ (2.5)

where $g(S_{n-1},S_{n})$ is the probability of no barrier hit within $[t_{n-1},t_{n}]$, conditional on $S_{n}\in(L_{n},U_{n})$ and $S_{n-1}\in(L_{n-1},U_{n-1})$. For a single barrier level $B_{n}$ (either lower $B_{n}=L_{n}$ or upper $B_{n}=U_{n}$) within $[t_{n-1},t_{n}]$,

 $g(S_{n-1},S_{n})=1-\exp\bigg({-}2\frac{\ln(S_{n}/B_{n})\ln(S_{n-1}/B_{n})}{% \sigma^{2}_{n}\delta t_{n}}\bigg);$ (2.6)

and there is a closed-form solution for the case of a double barrier within $[t_{n-1},t_{n}]$:

 $\displaystyle g(S_{n-1},S_{n})$ $\displaystyle=1-\sum_{m=1}^{\infty}[R_{n}(\alpha_{n}m-\gamma_{n},x_{n})+R_{n}(% -\alpha_{n}m+\beta_{n},x_{n})]$ $\displaystyle\qquad+\sum_{m=1}^{\infty}[R_{n}(\alpha_{n}m,x_{n})+R_{n}(-\alpha% _{n}m,x_{n})],$ (2.7)

where

 $\displaystyle x_{n}=\ln\frac{S_{n}}{S_{n-1}},\quad\alpha_{n}=2\ln\frac{U_{n}}{% L_{n}},\quad\beta_{n}=2\ln\frac{U_{n}}{S_{n-1}},$ $\displaystyle\gamma_{n}=2\ln\frac{S_{n-1}}{L_{n}},\quad R_{n}(z,x)=\exp\bigg({% -}\frac{z(z-2x)}{2\sigma^{2}_{n}\delta t_{n}}\bigg).$

Typically, a small number of terms in the above summations are enough to obtain high accuracy (in the actual implementation, the number of terms can be adapted to achieve the required accuracy; the smaller the time step $\delta t_{n}$, the fewer terms are needed). Formulas (2.6) and (2.1) can easily be obtained from the well-known distribution of the maximum and minimum of a Brownian motion (see, for example, Borodin and Salminen 1996; Karatzas and Shreve 1991); these formulas can also be found in Shevchenko (2011).

The integral (2.5) can be rewritten as

 $\displaystyle Q_{\mathrm{C}}$ $\displaystyle=B_{0,T}\int_{0}^{\infty}\mathrm{d}s_{1}\,f(s_{1}\mid s_{0})g(s_{% 0},s_{1})1_{(L_{1},U_{1})}(s_{1})\cdots\int_{0}^{\infty}\mathrm{d}s_{N}\,f(s_{% N}\mid s_{N-1})g(s_{N-1},s_{N})h(s_{N})1_{(L_{U},U_{N})}(s_{N})$ $\displaystyle=B_{0,T}\mathbb{E}\bigg(h(S_{N})\prod_{n=1}^{N}(1_{(L_{n},U_{n})}% (S_{n})g(S_{n-1},S_{n}))\bigg).$ (2.8)

An alternative expression for the barrier option that might provide a more efficient numerical estimate is presented by (2.11) in the next section. It is not analyzed in this paper and is the subject of a future study.

### 2.2 An alternative solution for the barrier option

The integral for barrier option price (2.1) can be rewritten in terms of the Markov chain $\smash{\hat{S}_{n}}$, starting at $\smash{\hat{S}_{0}=S_{0}}$, with elementary transitions

 $\Pr(\hat{S}_{n}\in\mathrm{d}s_{n}\mid\hat{S}_{n-1}=s_{n-1}):=\frac{\Pr(S_{n}% \in\mathrm{d}s_{n}\mid S_{n-1}=s_{n-1})1_{(L_{n},U_{n})}(s_{n})}{\Pr(S_{n}\in(% L_{n},U_{n})\mid S_{n-1}=s_{n-1})}.$ (2.9)

 $\hat{S}_{n}=\hat{S}_{n-1}\exp(a_{n}+b_{n}\hat{Z}_{n}),$ (2.10)

with

 $a_{n}:=(\mu_{n}-\tfrac{1}{2}\sigma^{2}_{n})\delta t_{n}\quad\text{and}\quad b_% {n}:=\sigma_{n}\sqrt{\delta t_{n}}.$

In addition, given the state variable $\smash{\hat{S}_{n-1}}$, $\smash{\hat{Z}_{n}}$ stands for a standard Gaussian random variable restricted to the set $\smash{(A_{n}(\hat{S}_{n-1}),B_{n}(\hat{S}_{n-1}))}$, with

 $A_{n}(\hat{S}_{n-1}):=\bigg[\ln{\bigg(\frac{L_{n}}{\hat{S}_{n-1}}\bigg)}-a_{n}% \bigg]\frac{1}{b_{n}}\quad\text{and}\quad B_{n}(\hat{S}_{n-1}):=\bigg[\ln{% \bigg(\frac{U_{n}}{\hat{S}_{n-1}}\bigg)}-a_{n}\bigg]\frac{1}{b_{n}}.$

Let

 $\varPhi(x):=\int_{-\infty}^{x}\frac{1}{\sqrt{2\pi}}\mathrm{e}^{-y^{2}/2}\,% \mathrm{d}y$

be the standard normal (Gaussian) distribution function, and let $\varPhi^{-1}(\cdot)$ be the inverse function. In this notation, we have that

 $\displaystyle\Pr(S_{n}\in(L_{n},U_{n})\mid S_{n-1}=s_{n-1})$ $\displaystyle=\Pr(Z_{n}\in(A_{n}(s_{n-1}),B_{n}(s_{n-1}))\mid S_{n-1}=s_{n-1})$ $\displaystyle=\varPhi(B_{n}(s_{n-1}))-\varPhi(A_{n}(s_{n-1})).$

We can also simulate the transition $\smash{\hat{S}_{n-1}\leadsto\hat{S}_{n}}$ by sampling a uniform random variable $\mathcal{U}_{n}$ by taking

 $\hat{Z}_{n}:=\varPhi^{-1}[\varPhi(A_{n}(\hat{S}_{n-1}))+\mathcal{U}_{n}~{}(% \varPhi(B_{n}(\hat{S}_{n-1}))-\varPhi(A_{n}(\hat{S}_{n-1})))]$

in (2.10). If we set

 $\varphi_{k-1}(s_{k-1}):=\Pr(S_{k}\in(L_{k},U_{k})\mid S_{k-1}=s_{k-1})=\varPhi% (B_{k}(s_{k-1}))-\varPhi(A_{k}(s_{k-1})),$

then we have that

 $\displaystyle\bigg\{\prod_{k=1}^{n}1_{(L_{k},U_{k})}(s_{k})\bigg\}\bigg\{\prod% _{k=1}^{n}\Pr(S_{k}\in\mathrm{d}s_{k}\mid S_{k-1}=s_{k-1})\bigg\}$ $\displaystyle=\bigg\{\prod_{k=1}^{n}\Pr(S_{k}\in(L_{k},U_{k})\mid S_{k-1}=s_{k% -1})\bigg\}\bigg\{\prod_{k=1}^{n}\Pr(\hat{S}_{k}\in\mathrm{d}s_{k}\mid\hat{S}_% {k-1}=s_{k-1})\bigg\}$ $\displaystyle=\bigg\{\prod_{k=1}^{n}{\varphi}_{k-1}(s_{k-1})\bigg\}\bigg\{% \prod_{k=1}^{n}\Pr(\hat{S}_{k}\in\mathrm{d}s_{k}\mid\hat{S}_{k-1}=s_{k-1})% \bigg\},$

from which we conclude that

 $Q_{\mathrm{C}}=B_{0,T}\mathbb{E}\bigg(h(\hat{S}_{N})\prod_{n=1}^{N}~{}\hat{G}_% {n-1}(\hat{S}_{n-1},\hat{S}_{n})\bigg),$ (2.11)

with the $[0,1]$-valued potential functions

 $\hat{G}_{n-1}(\hat{S}_{n-1},\hat{S}_{n}):=\varphi_{n-1}(\hat{S}_{n-1})g(\hat{S% }_{n-1},\hat{S}_{n}).$ (2.12)

Explicitly, the option price integral becomes

 $\displaystyle Q_{\mathrm{C}}$ $\displaystyle=B_{0,T}\int_{0}^{1}\mathrm{d}w_{1}(\varPhi(\tilde{U}_{1})-% \varPhi(\tilde{L}_{1}))g(s_{0},s_{1})\cdots\int_{0}^{1}\mathrm{d}w_{N}(\varPhi% (\tilde{U}_{N})-\varPhi(\tilde{L}_{N}))g(s_{N-1},s_{N})h(s_{N})$ $\displaystyle=B_{0,T}\int_{0}^{1}\cdots\int_{0}^{1}\mathrm{d}w_{1}\cdots% \mathrm{d}w_{N}h(s_{N})\prod_{n=1}^{N}(\varPhi(\tilde{U}_{n})-\varPhi(\tilde{L% }_{n}))g(s_{n-1},s_{n}),$ (2.13)

where

 $\displaystyle\tilde{U}_{n}$ $\displaystyle=\bigg(\ln\bigg(\frac{U_{n}}{s_{n-1}}\bigg)-(\mu_{n}-\tfrac{1}{2}% \sigma^{2}_{n})\delta t_{n}\bigg)\frac{1}{\sigma_{n}\sqrt{\delta t_{n}}},$ $\displaystyle\tilde{L}_{n}$ $\displaystyle=\bigg(\ln\bigg(\frac{L_{n}}{s_{n-1}}\bigg)-(\mu_{n}-\tfrac{1}{2}% \sigma^{2}_{n})\delta t_{n}\bigg)\frac{1}{\sigma_{n}\sqrt{\delta t_{n}}},$ $\displaystyle z_{n}$ $\displaystyle=\varPhi^{-1}[\varPhi(\tilde{L}_{n})+w_{n}(\varPhi(\tilde{U}_{n})% -\varPhi(\tilde{L}_{n}))],$ $\displaystyle s_{n}$ $\displaystyle=s_{n-1}\exp((\mu_{n}-\tfrac{1}{2}\sigma^{2}_{n})\delta t_{n}+% \sigma_{n}\sqrt{\delta t_{n}}z_{n})$

are calculated from $w_{1},\dots,w_{N}$ recursively for $n=1,2,\dots,N$ for given $s_{0}$.

This alternative solution for the barrier option might provide a more efficient numerical estimate but it is not analyzed in this paper.

## 3 Feynman–Kac representations

In this section, we provide the basic option price formulas under the Feynman–Kac representation underlying the SMC method; for a detailed introduction to this topic, see Carmona et al (2012).

### 3.1 Description of the models

Given that the transition-valued sequence

 $X_{n}=(S_{n},S_{n+1}),\quad n=0,\dots,N-1,$

forms a Markov chain, the option price expectation in the case of a continuously monitored barrier (2.1) can be written as

 $Q_{\mathrm{C}}=B_{0,T}\mathbb{E}\bigg(H(X_{N})\prod_{n=0}^{N-1}G_{n}(X_{n})% \bigg),$ (3.1)

with the extended payoff functions

 $H(X_{N})=H(S_{N},S_{N+1}):=h(S_{N})$

and the potential functions

 $G_{n}(X_{n})=g(S_{n},S_{n+1})1_{(L_{n+1}U_{n+1})}(S_{n+1}),\quad n=0,1,\dots,N% -1.$

These potential functions measure the chance to stay within the barriers during the interval $[t_{p},t_{p+1}]$. Equation (3.1) is the Feynman–Kac formula for discrete-time models (see Carmona et al 2012), which is used to develop the SMC option price estimator.

In this notation, the discretely monitored barrier option expectation (2.4) takes the following form:

 $Q_{\mathrm{D}}=B_{0,T}\mathbb{E}\bigg(H(X_{N})\prod_{n=0}^{N-1}\tilde{G}_{n}(X% _{n})\bigg),$ (3.2)

with the indicator potential functions

 $\tilde{G}_{n}(X_{n})=1_{(L_{n+1}U_{n+1})}(S_{n+1}),\quad n=0,1,\dots,N-1.$

We end this section with a Feynman–Kac representation of the alternative formulas for barrier option expectation presented by (2.11). In this case, if we consider the transition-valued Markov chain sequence

 $\hat{X}_{n}=(\hat{S}_{n},\hat{S}_{n+1}),\quad n=0,\dots,N-1,$

based on the modified underlying asset process $\hat{S}_{n}$ given by (2.10), then we can rewrite (2.11) as follows:

 $Q_{\mathrm{C}}=B_{0,T}\mathbb{E}\bigg(H(\hat{X}_{N})\prod_{n=0}^{N-1}\hat{G}_{% n}(\hat{X}_{n})\bigg),$ (3.3)

with the potential function $\hat{G}_{n}$ defined in (2.12). We observe that the above expression has exactly the same form as (3.1) by replacing $(X_{n},G_{n})$ with $(\hat{X}_{n},\hat{G}_{n})$.

Once the option price expectation is written as a Feynman–Kac representation then it is straightforward to develop SMC estimators as described in the following sections.

### 3.2 Some preliminary results

In this section, we review some key formulas related to unnormalized Feynman–Kac models. We provide a brief description of the evolution semigroup of Feynman–Kac measures. This section also presents some key multiplicative formulas describing the normalizing constants in terms of normalized Feynman–Kac measures. These mathematical objects are essential to define and analyze particle approximation models. For instance, the particle approximations of normalizing constants are defined by mimicking the multiplicative formula discussed above, replacing the normalized probability distributions by the empirical measures of the particle algorithm. We also emphasize that the bias and the variance analysis of these particle approximations are described in terms of the Feynman–Kac semigroups. A more thorough discussion on these stochastic models is provided in Del Moral (2004), Section 2.7.1 and Del Moral (2013), Section 3.2.2.

First, we observe that (3.1) can be written in the following form:

 $Q_{\mathrm{C}}=B_{0,T}\gamma_{N}(H)=B_{0,T}\gamma_{N}(1)\eta_{N}(H),$ (3.4)

with the Feynman–Kac unnormalized $\gamma_{N}$ and normalized $\eta_{N}$ measures given for any function $\varphi$ by the formulas

 $\gamma_{N}(\varphi)=\mathbb{E}\bigg(\varphi(X_{N})\prod_{n=0}^{N-1}G_{n}(X_{n}% )\bigg)\quad\text{and}\quad\eta_{N}(\varphi)=\frac{\gamma_{N}(\varphi)}{\gamma% _{N}(1)}.$ (3.5)

Note that the sequence of nonnegative measures $\smash{(\gamma_{n})_{n\geq 0}}$ satisfies, for any bounded measurable function $\varphi$, the recursive linear equation

 $\gamma_{n}(\varphi)=\gamma_{n-1}(\mathcal{Q}_{n}(\varphi)),$ (3.6)

with the integral operator

 $\mathcal{Q}_{n}(\varphi)(x)=G_{n-1}(x)K_{n}(\varphi)(x),$ (3.7)

where

 $K_{n}(\varphi)(x)=\mathbb{E}(\varphi(X_{n})\mid X_{n-1}=x)=\int K_{n}(x,% \mathrm{d}y)\varphi(y)$ (3.8)

and $K_{n}(X_{n-1},\mathrm{d}x):=\Pr(X_{n}\in\mathrm{d}x\mid X_{n-1})$ is the Markov transition in the chain $X_{n}$.

We prove this claim using the fact that

 $\displaystyle\gamma_{n}(\varphi)$ $\displaystyle=\mathbb{E}\bigg(\mathbb{E}\bigg(\varphi(X_{n})\prod_{p=0}^{n-1}G% _{p}(X_{p})\biggm|(X_{0},\dots,X_{n-1})\bigg)\bigg)$ $\displaystyle=\mathbb{E}\bigg(\mathbb{E}\bigg(\varphi(X_{n})\mid(X_{0},\dots,X% _{n-1})\bigg)\prod_{p=0}^{n-1}G_{p}(X_{p})\bigg)$ $\displaystyle=\mathbb{E}\bigg(\mathbb{E}(\varphi(X_{n})\mid X_{n-1})\prod_{p=0% }^{n-1}G_{p}(X_{p})\bigg)$ $\displaystyle=\mathbb{E}\bigg(G_{n-1}(X_{n-1})K_{n}(\varphi)(X_{n-1})\prod_{p=% 0}^{n-2}G_{p}(X_{p})\bigg)$ $\displaystyle=\gamma_{n-1}(G_{n-1}K_{n}(\varphi)).$ (3.9)

By construction, we also have that

 $\displaystyle\gamma_{N}(1)=\mathbb{E}\bigg(\prod_{n=0}^{N-1}G_{n}(X_{n})\bigg)% =\mathbb{E}\bigg(G_{N-1}(X_{N-1})\prod_{n=0}^{N-2}G_{n}(X_{n})\bigg)=\gamma_{N% -1}(G_{N-1}).$ (3.10)

This yields

 $\gamma_{N}(1)=\gamma_{N-1}(1)\frac{\gamma_{N-1}(G_{N-1})}{\gamma_{N-1}(1)}=% \gamma_{N-1}(1)\eta_{N-1}(G_{N-1}),$ (3.11)

from which we conclude that

 $\gamma_{N}(1)=\prod_{0\leq n (3.12)

and therefore

 $Q_{\mathrm{C}}=B_{0,T}\bigg[\prod_{0\leq n (3.13)

which is used for SMC estimators by replacing $\eta_{n}$ with its empirical approximation, as described in the following sections.

## 4 Monte Carlo estimators

In this section, we present Monte Carlo and SMC estimators and their corresponding algorithms to calculate the option price in the case of continuously and discretely monitored barrier conditions.

### 4.1 Standard Monte Carlo

Using process (2.3), simulate independent asset path realizations

 $\bm{S}^{(m)}=(S_{1}^{(m)},\dots,S_{N}^{(m)}),\quad m=1,\dots,M.$

Then, the unbiased estimator for the continuously monitored barrier option price integral (2.1) is a standard average of option price payoff realizations over the simulated paths

 $\displaystyle\hat{Q}^{\text{MC}}_{\mathrm{C}}$ $\displaystyle=B_{0,T}\frac{1}{M}\sum_{m=1}^{M}\bigg(h(S_{N}^{(m)})\prod_{n=1}^% {N}\{g(S^{(m)}_{n-1},S^{(m)}_{n})1_{[L_{n},U_{n}]}(S^{(m)}_{n})\}\bigg)$ $\displaystyle=B_{0,T}\frac{1}{M}\sum_{m=1}^{M}\bigg(H(X_{N}^{(m)})\prod_{n=0}^% {N-1}G_{n}(X_{n}^{(m)})\bigg),$ (4.1)

with $\smash{X_{n}^{(m)}=(S^{(m)}_{n},S^{(m)}_{n+1})}$, and the unbiased estimator for discretely monitored barrier option (2.4) is

 $\hat{Q}^{\text{MC}}_{\mathrm{D}}=B_{0,T}\frac{1}{M}\sum_{m=1}^{M}\bigg(H(X_{N}% ^{(m)})\prod_{n=0}^{N-1}\tilde{G}_{n}(X_{n}^{(m)})\bigg),$ (4.2)

with $X_{n}^{(m)}=(S^{(m)}_{n},S^{(m)}_{n+1})$.

### 4.2 Sequential Monte Carlo

Another unbiased estimator for option price integral (2.1) can be obtained using (3.13) via the SMC method with the following algorithm.

###### Algorithm 4.1.
1. (1)

Initial step

1. (a)

(Proposition.) For the initial time step, $I_{0}=[t_{0},t_{1}]$, simulate $M$ independent realizations

 $X_{0}^{(m)}:=(S_{0}^{(m)},S_{1}^{(m)}),\quad m=1,\dots,M,$

using process (2.3); these are referred to as $M$ (transition-type) particles. Set

 $G_{0}(X_{0}^{(m)})=1_{(L_{1},U_{1})}(S^{(m)}_{1})g(S_{0}^{(m)},S^{(m)}_{1})$

for each $1\leq m\leq M$.

2. (b)

(Acceptance–rejection.) Sample $M$ random and $[0,1]$-valued uniform variables $\smash{U^{(m)}_{0}}$. The rejected transition-type particles $\smash{X_{0}^{(m)}}$ are those for which $\smash{G_{0}(X_{0}^{(m)})}<\smash{U^{(m)}_{0}}$. The particles $\smash{X_{0}^{(m)}}$ for which $\smash{U^{(m)}_{0}}\leq\smash{G_{0}(X_{0}^{(m)})}$ are accepted. Note that a transition-type particle $\smash{X_{0}^{(m)}}$ such that $\smash{S_{1}^{(m)}}\not\in(L_{1},U_{1})$ is instantly rejected (since its weight $\smash{G_{0}(X_{0}^{(m)})}=0$ is null), and a transition-type particle $\smash{X_{0}^{(m)}}$ such that $\smash{S_{1}^{(m)}}\in(L_{1},U_{1})$ is rejected with a probability $1-\smash{G_{0}(X_{0}^{(m)})}$.

3. (c)

(Recycling–selection.) Resample each rejected transition-type particle $\smash{X_{0}^{(m)}}$ by resampling its $\smash{S_{1}^{(m)}}$ component from the discrete distribution with density function

 $f(s_{1})=\sum_{m=1}^{M}\frac{G_{0}(X_{0}^{(m)})}{\sum_{k=1}^{M}G_{0}(X_{0}^{(k% )})}\delta(s_{1}-S_{1}^{(m)}),$ (4.3)

where $\delta(y-y_{0})$ is a point mass function centered at $y_{0}$ (ie, the Dirac $\delta$-function, which is zero everywhere except from $y=y_{0}$, and its integral over any interval containing $y_{0}$ is equal to 1). In other words, when a transition-type particle, say $\smash{X_{0}^{(r)}}$, is rejected for some index $r$, we replace it by one of the particles $X_{0}^{(m)}$ randomly chosen with respect to its weight

 $\frac{G_{0}(X_{0}^{(m)})}{\sum_{k=1}^{M}G_{0}(X_{0}^{(k)})}.$

Efficient and simple sampling of the rejected particle from the discrete density (4.3) can be accomplished by Algorithm 4.4 in Section 4.3.

At the end of the acceptance–rejection-recycling scheme, we have $M$ (transition-type) particles that we denote by

 $\tilde{X}_{0}^{(m)}=(S_{0}^{(m)},\tilde{S}_{1}^{(m)}),\quad m=1,\dots,M.$
###### Remark 4.2.

By definition of (4.3) we note that transition-type particles $\smash{X_{0}^{(m)}}$ such that $\smash{S_{1}^{(m)}\not\in(L_{1},U_{1})}$ have a null weight. Therefore, they cannot be selected to replace the rejected ones. Moreover, the transition-type particles $\smash{X_{0}^{(m)}}$ such that $\smash{S_{1}^{(m)}\in(L_{1},U_{1})}$ with a large probability $\smash{G_{0}(X_{0}^{(m)})}$ of not hitting the barrier within $[t_{0},t_{1}]$ are more likely to be selected (in replacement of the rejected ones).

2. (2)

Step $0\leadsto 1$

1. (a)

(Proposition.) For the second time step, $I_{1}=[t_{1},t_{2}]$, simulate $M$ independent realizations

 $X_{1}^{(m)}:=(\tilde{S}_{1}^{(m)},S_{2}^{(m)}),\quad m=1,\dots,M,$

starting from the end points $\smash{\tilde{S}_{1}^{(m)}}$ of the selected transitions $\smash{\tilde{X}_{0}^{(m)}}$ at the previous step, using the process evolution (2.3); these are referred to as $M$ (transition-type) particles $X_{1}^{(m)}$ at time $1$. Set

 $G_{1}(X_{1}^{(m)})=1_{(L_{2},U_{2})}(S^{(m)}_{2})~{}\times~{}g(\tilde{S}_{1}^{% (m)},S_{2}^{(m)})$

for each $1\leq m\leq M$.

2. (b)

(Acceptance–rejection.) Sample $M$ random and $[0,1]$-valued uniform variables $\smash{U^{(m)}_{1}}$. The rejected transition-type particles $\smash{X_{1}^{(m)}}$ are those for which $\smash{G_{1}(X_{1}^{(m)})}<\smash{U^{(m)}_{1}}$ and the particles $\smash{X_{1}^{(m)}}$ for which

 $U^{(m)}_{1}\leq G_{1}(X_{1}^{(m)})$

are accepted. That is, a transition-type particle $\smash{X_{1}^{(m)}}$ such that $\smash{S_{2}^{(m)}}\not\in\smash{(L_{2},U_{2})}$ is instantly rejected and $\smash{X_{1}^{(m)}}$ such that $\smash{S_{2}^{(m)}}\in(L_{2},U_{2})$ is rejected with a probability $1-\smash{G_{1}(X_{1}^{(m)})}$.

3. (c)

(Recycling–selection.) Resample each rejected transition-type particle $X_{1}^{(m)}$ by resampling its $\smash{S_{2}^{(m)}}$ component from the discrete distribution, with density

 $f(s_{2})=\sum_{m=1}^{M}\frac{G_{1}(X_{1}^{(m)})}{\sum_{k=1}^{M}G_{1}(X_{1}^{(k% )})}\delta(s_{2}-S_{2}^{(m)}),$ (4.4)

using, for example, the efficient and simple Algorithm 4.4 in Section 4.3.

At the end of the acceptance–rejection-recycling scheme, we have $M$ (transition-type) particles denoted by $\smash{\tilde{X}_{1}^{(m)}=(\tilde{S}_{1}^{(m)},\tilde{S}_{2}^{(m)})}$, $1\leq m\leq M$. A remark similar to Remark 4.2 is also applied here: transition-type particles $\smash{X_{1}^{(m)}}$ such that $\smash{S_{2}^{(m)}\not\in(L_{2},U_{2})}$ have a null weight and therefore they cannot be selected in replacement of the rejected ones. Moreover, the transition-type particles $\smash{X_{1}^{(m)}}$ such that $\smash{S_{2}^{(m)}\in(L_{2},U_{2})}$ with a high probability $\smash{g(\tilde{S}_{1}^{(m)},S_{2}^{(m)})}$ of not hitting the barrier within $[t_{0},t_{1}]$ are more likely to be selected (to replace the rejected ones).

3. (3)

Repeat steps (a)–(c) in Step 2 for time steps $[t_{2},t_{3}],\dots,[t_{N-1},t_{N}]$.

Calculate the final unbiased option price estimator as

 $\displaystyle\hat{Q}^{\text{SMC}}_{\mathrm{C}}=B_{0,T}\bigg[\prod_{n=0}^{N-1}% \frac{1}{M}\sum_{m=1}^{M}G_{n}(X_{n}^{(m)})\bigg]\frac{1}{M}\sum_{m=1}^{M}H(X_% {N}^{(m)}).$ (4.5)

That is, $\eta_{N}(H)$ is replaced by its empirical approximation

 $\frac{1}{M}\sum_{m=1}^{M}H(X_{N}^{(m)})$

and $\eta_{n}(G_{n})$ is replaced by its empirical approximation

 $\frac{1}{M}\sum_{m=1}^{M}G_{n}(X_{n}^{(m)})$

in (3.13).

Note that $\smash{H(X_{N}^{(m)})=h(\tilde{S}_{N}^{(m)})}$, ie, payoff at maturity is calculated using particles $\smash{\tilde{S}_{N}^{(m)}}$ after rejection-recycling at maturity $t_{N}=T$. The proof of the unbiasedness properties of these estimators is provided in Section 4.4.

In much the same way, an unbiased estimator of $Q_{\mathrm{D}}$ defined in (3.2) is given by

 $\hat{Q}^{\text{SMC}}_{\mathrm{D}}=B_{0,T}\bigg[\prod_{n=0}^{N-1}\frac{1}{M}% \sum_{m=1}^{M}\tilde{G}_{n}(X_{n}^{(m)})\bigg]\frac{1}{M}\sum_{m=1}^{M}H(X_{N}% ^{(m)}),$ (4.6)

where $\smash{(X^{(m)}_{n})_{0\leq n\leq N}}$, $1\leq m\leq M$, is obtained by the above algorithm, with potential functions $\smash{(G_{n})_{0\leq n\leq N}}$ replaced by the indicator potential functions $\smash{(\tilde{G}_{n})_{0\leq n\leq N}}$.

In both cases, all the particles may exit the barrier after some proposition stage. In this case, we use the convention that the above estimates are null. One way to solve this problem is to consider the Feynman–Kac description (3.3) for the alternative option price expression (2.11) presented in Section 2.2. In this context, an unbiased estimator of $Q_{\mathrm{C}}$ is given by

 $\hat{\hat{Q}}^{\text{SMC}}_{\mathrm{C}}=B_{0,T}\bigg[\prod_{n=0}^{N-1}\frac{1}% {M}\sum_{m=1}^{M}\hat{G}_{n}(\hat{X}_{n}^{(m)})\bigg]\frac{1}{M}\sum_{m=1}^{M}% H(\hat{X}_{N}^{(m)}),\vspace*{2pt}$ (4.7)

where $\smash{(\hat{X}^{(m)}_{n})_{0\leq n\leq N}}$, $1\leq m\leq M$, is obtained by the above algorithm for $\smash{(X^{(m)}_{n})_{0\leq n\leq N}}$, with potential functions $\smash{(G_{n})_{0\leq n\leq N}}$ replaced by the potential functions $\smash{(\hat{G}_{n})_{0\leq n\leq N}}$, while the process for $S_{n}$ is replaced by the process $\hat{S}_{n}$ described in Section 2.2.

###### Remark 4.3.

As we mentioned in the introduction of Section 3.2, the particle estimate in (4.5) is defined as in (3.13) by replacing the normalized Feynman–Kac measures $\eta_{n}$ with particle empirical approximations. Formula (4.6) (respectively, (4.7)) follows the same line of argument as that based on the Feynman–Kac model (3.2) (respectively, (3.3)).

Figure 1 presents an illustration of the algorithm with $M=6$ particles. In this particular case, we simulate six particles at time $t_{1}$ (starting from $S_{0}$). Then particle $\smash{S_{1}^{(4)}}$ is rejected and resampled (moved to position $\smash{S_{1}^{(1)}}$), particle $\smash{S_{1}^{(6)}}$ is rejected and moved to position $\smash{S_{1}^{(3)}}$. Then two particles located at $\smash{S_{1}^{(3)}}$ will generate two particles at $t_{2}$, two particles located at $\smash{S_{1}^{(1)}}$ will generate two particles at $t_{2}$, etc. For each time slice, including the last $t_{N}$, after resampling, we have six particles above the barrier. Note that it is possible that $\smash{S_{1}^{(1)}}$, $\smash{S_{1}^{(2)}}$, $\smash{S_{1}^{(3)}}$ and $\smash{S_{1}^{(5)}}$ are also rejected in the case of a continuously monitored barrier.

### 4.3 Sampling from discrete distribution

For the benefit of the reader, in this section we present a simple, efficient algorithm for sampling the rejected particles from the discrete density required during the recycle–selection step of the SMC algorithm described in Section 4.2, ie, sampling from discrete densities (4.3) and (4.4).

In general, sampling of $R$ independent random variables $\smash{(Y^{(r)})_{1\leq r\leq R}}$ from a weighted discrete probability density function

 $f(x)=\sum_{m=1}^{M}p_{m}\delta(x-x_{m})$ (4.8)

can be done in the usual way by the inverse distribution method. That is,

 $F(x)=\frac{1}{M}\sum_{m=1}^{M}1_{[x_{m},\infty)}(x)$

is a distribution corresponding to discrete density (4.8), and $X=F^{-1}(U)$ is a sample from $F(x)$ if $U$ is from the uniform $(0,1)$ distribution. It is important to use a computationally efficient method for the sampling of $R$ variables. If the order of samples is not important (as in the case of the recycling–selection steps of the SMC algorithm in Section 4.2), then, for example, we can sample $(R+1)$ independent exponential random variables $\smash{(\mathcal{E}_{r})_{1\leq r\leq(R+1)}}$, with unit parameter and set

 $\mathcal{T}_{r}=\sum_{1\leq s\leq r}\mathcal{E}_{s}\quad\text{and}\quad% \mathcal{V}_{r}=\frac{\mathcal{T}_{r}}{\mathcal{T}_{R+1}},\quad r=1,2,\dots,R+1.$ (4.9)

The random variables $(\mathcal{V}_{1},\dots,\mathcal{V}_{R})$ calculated in such a way are the order statistics of $R$ independent random variables uniformly distributed on $(0,1)$, which is a well-known property of the Poisson process (see, for example, Bartoli and Del Moral (2001), Example 3.6.9 and Section 2.6.2 or Daley and Vere-Jones (2003), Exercise 2.1.2). Then, the sampling of $\smash{(Y^{(r)})_{1\leq r\leq R}}$, by calculating $Y^{(r)}=F^{-1}(\mathcal{V}_{r})$, can be accomplished using the following synthetic pseudocode.

###### Algorithm 4.4.
1. 1.

$k=1$ and $r=1$

2. 2.

While $r\leq R$

1. $\bullet$

While $\mathcal{V}_{r}

• $Y^{(r)}=x_{k}$

• $r=r+1$

2. $\bullet$

End while

3. $\bullet$

$k=k+1$

3. 3.

End while

The computational cost of this sampling scheme is linear with respect to $R$. In particular, to simulate from the probability density (4.3), set

 $p_{m}=\frac{G_{0}(X_{0}^{(m)})}{\sum_{k=1}^{M}G_{0}(X_{0}^{(k)})}\quad\text{% and}\quad x_{m}=S_{1}^{(m)},$

and to simulate from the discrete distribution (4.4), set

 $p_{m}=\frac{G_{1}(X_{1}^{(m)})}{\sum_{k=1}^{M}G_{1}(X_{1}^{(k)})}\quad\text{% and}\quad x_{m}=S_{2}^{(m)}$

in (4.8).

### 4.4 Unbiasedness properties

The SMC estimator for the option price (4.5) can be written as

 $\displaystyle\hat{Q}^{\text{SMC}}=B_{0,T}\gamma^{M}_{N}(1)\eta^{M}_{N}(H),$ (4.10)

with the empirical measures $\smash{\eta^{M}_{N}}$ given by

 $\eta^{M}_{N}(H)=\frac{1}{M}\sum_{m=1}^{M}H(X_{N}^{(m)})\$ (4.11)

and normalizing constants

 $\gamma^{M}_{N}(1)=\prod_{p=0}^{N-1}\frac{1}{M}\sum_{m=1}^{M}G_{p}(X_{p}^{(m)})% =\prod_{p=0}^{N-1}\eta^{M}_{p}(G_{p}).$ (4.12)

In this notation, the $M$-particle approximations of the Feynman–Kac measures $\gamma_{N}$ for any function $\varphi$ are given by

 $\gamma^{M}_{N}(\varphi):=\gamma^{M}_{N}(1)\times\eta^{M}_{N}(\varphi)\implies% \hat{Q}^{\text{SMC}}=B_{0,T}\times\gamma^{M}_{N}(H).$ (4.13)

Here, $\smash{\eta_{N}^{M}}$ and $\smash{\gamma_{N}^{M}}$ are particle empirical approximations of Feynman–Kac measures $\eta_{N}$ and $\gamma_{N}$ in the option price formula (3.13).

The objective of this section is to show that the $M$-particle estimates $\smash{\hat{Q}^{\text{SMC}}}$ for continuous and discrete cases (4.5) and (4.6) are unbiased. The unbiased property is not so obvious, mainly because it is based on biased $M$-empirical measures $\smash{\eta^{M}_{N}}$. It is clearly beyond the scope of this study to present a quantitative analysis of these biased measures; we refer the reader to the monographs by Del Moral (2004, 2013), and references therein. For instance, we can prove that

 $\sup_{\|\varphi\|\leq 1}\|\mathbb{E}(\eta^{M}_{N}(\varphi))-\eta_{N}(\varphi)% \|\leq\frac{c(N)}{M}$ (4.14)

for some finite positive constant $c(N)$ whose values only depend on the time horizon $N$. That is, $\smash{\eta^{M}_{N}(\varphi)}$ converges to $\smash{\eta_{N}(\varphi)}$ as $M$ increases. The unnormalized particle measures $\smash{\gamma_{N}^{M}}$ in (4.5)–(4.7) are unbiased. On the other hand, the empirical measures $\smash{\eta_{N}^{M}(\varphi)}$ can be expressed in terms of the ratio of two unnormalized quantities, $\smash{\gamma_{N}^{M}(\varphi)}$ and $\smash{\gamma_{N}^{M}(1)}$. Taking into consideration the fluctuation of these unnormalized particle models, the estimate of the bias (4.14) is obtained using an elementary Taylor-type expansion of the first order of this ratio.

To prove that $\smash{\gamma_{N}^{M}(H)}$ is unbiased, ie, $\smash{\hat{Q}^{\text{SMC}}}$ is unbiased, recall that the particles evolve sequentially using a selection and a mutation transition. Thus, we have the conditional expectation formula

 $\displaystyle\mathbb{E}(\eta^{M}_{N}(H)\mid(X_{0}^{(m)},\dots,X_{N-1}^{(m)})_{% 1\leq m\leq M})$ $\displaystyle=\mathbb{E}(H(X_{N}^{(1)})\mid(X_{0}^{(m)},\dots,X_{N-1}^{(m)})_{% 1\leq m\leq M})$ $\displaystyle=\sum_{1\leq m\leq M}\frac{G_{N-1}(X_{N-1}^{(m)})}{\sum_{1\leq k% \leq M}G_{N-1}(X_{N-1}^{(k)})}K_{N}(H)(X^{(m)}_{N-1}),$ (4.15)

where $K_{N}$ is the Markov transition integral operator of the chain $X_{n}^{(m)}$, $n=1,\dots,N-1$, defined in (3.8). The weighted mixture of Markov transitions expresses the fact that the particles are selected using the potential functions before exploring the solution space using the mutation transitions. This implies that

 $\displaystyle\mathbb{E}(\gamma^{M}_{N}(H)\mid(X_{0}^{(m)},\dots,X_{N-1}^{(m)})% _{1\leq m\leq M})$ $\displaystyle=\bigg[\prod_{p=0}^{N-1}\eta^{M}_{p}(G_{p})\bigg]\frac{1}{N}\sum_% {1\leq m\leq M}\frac{G_{N-1}(X_{N-1}^{(m)})}{(1/N)\sum_{1\leq k\leq M}G_{N-1}(% X_{N-1}^{(k)})}K_{N}(H)(X^{(m)}_{N-1})$ $\displaystyle=\bigg[\prod_{p=0}^{N-2}\eta^{M}_{p}(G_{p})\bigg]\eta^{M}_{N-1}(% \mathcal{Q}_{N}(H)),$ (4.16)

with the one-step Feynman–Kac semigroup $\smash{\mathcal{Q}_{N}}$ introduced in (3.7). That is,

 $\mathbb{E}(\gamma^{M}_{N}(H)\mid(X_{0}^{(m)},\dots,X_{N-1}^{(m)})_{1\leq m\leq M% })=\gamma^{M}_{N-1}(\mathcal{Q}_{N}(H)),$ (4.17)

and therefore

 $\mathbb{E}(\gamma^{M}_{N}(H))=\mathbb{E}(\gamma^{M}_{N-1}(\mathcal{Q}_{N}(H))).$ (4.18)

For $N=0$, we use the convention $\prod_{\emptyset}=1$ so that

 $\gamma^{M}_{0}=\eta^{M}_{0}\implies\mathbb{E}(\gamma^{M}_{0}(\varphi))=\mathbb% {E}(\eta^{M}_{0}(\varphi))=\eta_{0}(\varphi)=\gamma_{0}(\varphi)$

for any function $\varphi$. Iterating (4.17) backward in time, we obtain the evolution equation of the unnormalized Feynman–Kac distributions defined in (3.6). Next, for the convenience of the reader, we provide a more detailed proof of the unbiased property, and we further assume that

 $\mathbb{E}(\gamma^{M}_{n}(\varphi))=\gamma_{n}(\varphi)$ (4.19)

at some rank $n$, for any $M\geq 1$ and any $\varphi$. In this case, arguing as above, we have

 $\mathbb{E}(\gamma^{M}_{n+1}(\varphi))=\mathbb{E}(\gamma^{M}_{n}(\mathcal{Q}_{n% +1}(\varphi))).$ (4.20)

Under the induction hypothesis, this implies that

 $\mathbb{E}(\gamma^{M}_{n+1}(\varphi))=\gamma_{n}(\mathcal{Q}_{n+1}(\varphi))=% \gamma_{n+1}(\varphi).$ (4.21)

This ends the proof of the unbiasedness property of $\smash{\hat{Q}^{\text{SMC}}}$. The results about standard errors of these SMC unbiased estimators can be found, for example, in Cérou et al (2011). While it is beyond the scope of this paper to go into the details of the theoretical results on the variance of empirical approximations of normalized Feynman–Kac measures, it is important to note that the standard error of the SMC estimator is proportional to $1/\smash{\sqrt{M}}$, which is the same as for the standard Monte Carlo estimator. However, while for the Monte Carlo estimator the proportionality coefficient is easily estimated as the standard deviation of the simulated asset path payoffs, for the SMC estimator there is no simple expression, and we have to run independent calculations of the SMC estimator to find its standard error; numerical experiments are presented in Section 6.

## 5 Importance-sampling models

The Feynman–Kac representation formulas (3.1) and their particle interpretations discussed in Section 4.2 are far from being unique. For instance, using (2.1), for any nonnegative probability density functions $\bar{f}(s_{n}\mid s_{n-1})$, we also have that

 $Q=B_{0,T}\int_{0}^{\infty}\mathrm{d}s_{1}\,\bar{f}(s_{1}\mid s_{0})\bar{g}(s_{% 0},s_{1})1_{(L_{1},U_{1})}(s_{1})\cdots\int_{0}^{\infty}\mathrm{d}s_{N}\,\bar{% f}(s_{N}\mid s_{N-1})\bar{g}(s_{N-1},s_{N})h(s_{N})1_{(L_{U},U_{N})}(s_{N}),$ (5.1)

with the potential functions

 $\bar{g}(\bar{S}_{n-1},\bar{S}_{n})={g}(\bar{S}_{n-1},\bar{S}_{n})\frac{f(s_{n}% \mid s_{n-1})}{\bar{f}(s_{n}\mid s_{n-1})}.$ (5.2)

This yields the Feynman–Kac representation

 $Q=B_{0,T}\mathbb{E}\bigg(h(\bar{S}_{N})\prod_{n=1}^{N}\bar{G}_{n}(\bar{S}_{n-1% },\bar{S}_{n})\bigg)$ (5.3)

in terms of the potential functions

 $\bar{G}_{n}(\bar{S}_{n-1},\bar{S}_{n})=1_{(L_{n},U_{n})}(\bar{S}_{n})\bar{g}(% \bar{S}_{n-1},\bar{S}_{n})$ (5.4)

and the Markov chain $\smash{(\bar{S}_{n})_{n\geq 0}}$, with

 $\Pr(\bar{S}_{n}\in\mathrm{d}s_{n}\mid\bar{S}_{n-1})=\bar{f}(s_{n}\mid\bar{S}_{% n-1})\,\mathrm{d}s_{n}.$ (5.5)

The importance-sampling formula (5.3) is rather well known. The corresponding $M$-particle model consists of $M$ particles evolving, between the selection times, as independent copies of the twisted Markov chain model $\bar{S}_{n}$; the selection–recycling procedure favors transitions $\smash{\bar{S}_{n-1}}\leadsto\smash{\bar{S}_{n}}$ that increase the density ratio $\smash{f(\bar{S}_{n}\mid\bar{S}_{n-1})/\bar{f}(\bar{S}_{n}\mid\bar{S}_{n-1})}$.

We end this section with a more sophisticated change of measure related to the payoff functions.

For any sequence of positive potential functions $\smash{(h_{n})_{0\leq n\leq N}}$ with $\smash{h_{N}=h}$, and using the fact that

 $h(\bar{S}_{N})=\frac{h_{N}(\bar{S}_{N})}{h_{N-1}(\bar{S}_{N-1})}\frac{h_{N-1}(% \bar{S}_{N-1})}{h_{N-2}(\bar{S}_{N-2})}\cdots\frac{h_{1}(\bar{S}_{1})}{h_{0}(% \bar{S}_{0})}h_{0}(\bar{S}_{0}),$ (5.6)

we have that

 $\displaystyle Q_{0}$ $\displaystyle=B_{0,T}h_{0}(\bar{s}_{0})\mathbb{E}\bigg(\prod_{n=1}^{N}\bigg(% \frac{h_{n}(\bar{S}_{n})}{h_{n-1}(\bar{S}_{n-1})}\bar{G}_{n}(\bar{S}_{n-1},% \bar{S}_{n})\bigg)\bigg)$ $\displaystyle=B_{0,T}h_{0}(\bar{s}_{0})\mathbb{E}\bigg(\prod_{n=1}^{N}\check{G% }_{n}(\bar{S}_{n-1},\bar{S}_{n})\bigg),$ (5.7)

with

 $\check{G}_{n}(\bar{S}_{n-1},\bar{S}_{n})=\bar{G}_{n}(\bar{S}_{n-1},\bar{S}_{n}% )\frac{h_{n}(\bar{S}_{n})}{h_{n-1}(\bar{S}_{n-1})}.$ (5.8)

For example, for the payoff functions discussed in the option pricing model (2.2), we can choose

 $h_{N}(x)=h(x)=\max(K-x,0)\quad\text{and}\quad h_{n}(x)=h(x)+1\quad\forall n (5.9)

Note that the $M$-particle model associated with the potential functions $\check{G}_{n}$ consists of $M$ particles evolving, between the selection times, as independent copies of the Markov chain $\bar{S}_{n}$; and the selection–recycling procedure favors transitions $\smash{\bar{S}_{n-1}}\leadsto\smash{\bar{S}_{n}}$ that increase the ratio $h_{n}(\bar{S}_{n})/h_{n-1}(\bar{S}_{n-1})$. For instance, in the example suggested in (5.9), the transitions $\smash{\bar{S}_{n-1}\leadsto\bar{S}_{n}}$ exploring regions far from the strike $K$ are more likely to duplicate.

The choice of potential functions (5.2) allows us to choose the reference Markov chain in order to randomly explore the state space during the mutation transitions. The importance-sampling Feynman–Kac model (5) is less intrusive. More precisely, without changing the reference Markov chain, the choice of potential functions (5.8) allows us to favor transitions that sequentially increase the payoff function. The importance-sampling models (5.2) and (5.8) can obviously be combined so as to change the reference Markov chain and favor the transitions that increase the payoff function.

## 6 Numerical results

Consider a simple knock-out barrier call option with constant lower and upper barriers $L=90$ and $U=110$, strike $K=100$ and maturity $T=0.5$ for the following market data: spot $S_{0}=100$, interest rate $r=0.1$, volatility $\sigma=0.3$ and zero dividends $q=0$. Tables 1 and 2 and Figures 2 and 3 show the exact closed-form solutions, SMC and standard Monte Carlo estimators, standard errors of the estimators and estimator efficiencies for this option for continuously and discretely monitored barrier cases. We perform $M=100\,000$ simulations for Monte Carlo estimators and $M=100\,000$ particles for SMC estimators, which are repeated fifty times (using independent random numbers) to calculate the final option price estimates and their standard errors.

Our calculations are based on sampling at equally spaced time slices $t_{1},\dots,t_{N}$. Note that we present results for $N=(1,2,4,8,16,32,64,128)$ not to demonstrate the convergence of a discretely monitored barrier to the continuous case or to address time discretization errors, but rather to illustrate and explain the behavior of the SMC, which improves the accuracy of the option price sampling estimator for a given time discretization. In the case of a real barrier option, the time discretization will be dictated by the stochastic process, window barrier structure, barrier monitoring type (eg, continuous, daily) and market data term structures.

For the Monte Carlo estimator (in the case of a continuously monitored barrier), we need to calculate the conditional probability of a barrier hit (2.1) between sampled dates only for asset simulated paths that do not breach the barrier condition during the option’s life and result in a nonzero payoff at maturity. For SMC estimators, these probabilities should be calculated for all time steps but only for particles that appear between the barriers. Thus, the direct calculation of computational effort is not straightforward. Instead, we can use the actual computing time to compare the methods using the following facts.

1. (1)

The CPU time $t_{\text{cpu}}$ is proportional to the number of simulations $M$ in the Monte Carlo method (or the number of particles $M$ in SMC).

2. (2)

Both Monte Carlo and SMC estimators are unbiased. Their standard errors are proportional to $1/\smash{\sqrt{M}}$, with different proportionality coefficients for SMC and Monte Carlo (for theoretical results about the variance of SMC estimators, see Cérou et al (2011)). While, for Monte Carlo, this coefficient is easily calculated as the standard deviation of asset path payoffs, for SMC there is no simple expression, and we must run independent calculations many times (eg, fifty times in our numerical example) to estimate the standard errors of the SMC estimators.

Thus, the squared standard error $s^{2}$ of an estimator is

 $s^{2}=\frac{\alpha}{t_{\text{cpu}}},$ (6.1)

where $\alpha$ depends on the method, ie, $\alpha=\alpha_{\text{MC}}$ for Monte Carlo and $\alpha=\alpha_{\text{SMC}}$ for SMC; these are easily found from the numerical results for $s^{2}$ and $t_{\text{cpu}}$ of the corresponding estimators. To compare the efficiency of the estimators, we calculate

 $\kappa=\frac{\alpha_{\text{MC}}}{\alpha_{\text{SMC}}}.$ (6.2)

Interpretation of $\kappa$ is straightforward: if the computing time for the SMC estimator is $t_{\text{SMC}}$, then the computing time taken for the Monte Carlo estimator to achieve the same accuracy as the SMC estimator is $\kappa\times t_{\text{SMC}}$, ie, $\kappa>1$ indicates that SMC is faster than Monte Carlo, and $\kappa<1$ indicates otherwise.

For our specific numerical example, the computing time for SMC is about only 10–20% longer than for Monte Carlo in the case of the discretely monitored barrier. In the case of the continuously monitored barrier, the SMC’s time is about double that of the Monte Carlo’s time, mainly because we need to calculate the conditional probability of barrier hit (2.1) between sampled dates, which is computationally expensive in the case of a double barrier. However, the standard error for the SMC estimator is always smaller than that for the Monte Carlo estimator (except in the limiting case $N=1$, where the barrier is monitored at maturity only when both standard errors are about the same). It is easy to see from our results that the SMC is superior to the Monte Carlo (except for when $N=1$). For both the discrete and continuous barrier cases we observe that the SMC efficiency coefficient $\kappa$ monotonically increases as the number of time steps $N$ increases. The accuracy (standard error) of the SMC estimator does not change much as $N$ increases because barrier-rejected asset-sampled values (particles) are resampled from particles between the barriers, and thus at maturity we still have $M$ particles between the barriers, regardless of $N$. The standard error of Monte Carlo estimator grows with $N$ because the number of simulated paths that will reach maturity without breaching the barrier condition will reduce as $N$ increases.

It is easy to see from Table 2 that in the case of a discretely monitored barrier the SMC efficiency $\kappa$ is roughly proportional to $1/\psi$, where $\psi$ is the probability of the underlying asset not hitting the barrier during the option’s life (ie, in this case it is the probability of the asset path reaching maturity without breaching the barrier). Note that in the case of a discretely monitored barrier, $\psi$ decreases as the number of time steps $N$ increases (ie, fewer paths will reach maturity without breaching the barrier as $N$ increases). In the case of a continuously monitored barrier, $\psi$ does not change with $N$ (it is about 0.5% for the option calculated in our numerical example; see Table 1). However, note that the Monte Carlo estimator for the continuously monitored barrier case is calculated by sampling asset paths over $N$ dates and multiplying the path payoff at maturity by the conditional probabilities of not hitting the barrier between sampled dates (2.1). Thus, the probability of the asset paths reaching maturity without breaching the barrier is the same as for the discrete barrier case. As a result, the standard error of the Monte Carlo estimator (for both discrete and continuous barriers) grows as $N$ increases.

Other numerical experiments not reported here show that the efficiency of the SMC over the Monte Carlo improves when the barriers become closer, ie, when the probability of the asset path hitting the barrier increases; this is also easy to see from the results in Table 2. If the probability of the asset path not hitting the barrier is high, then the performance of the SMC is about the same as – or slightly worse than – the Monte Carlo. Note that our implementation does not include any standard variance reduction techniques such as antithetics, importance sampling or control variates, or any parallel/vector computations. The algorithm was implemented using Fortran 90 and executed on a standard laptop (Windows 7, Intel i7-2640M CPU at 2.8 GHz, 4 GB RAM). While computing time is somewhat subjective (ie, it depends on specifics of our implementation), the ratio of standard errors (or ratio of squared standard errors) of the Monte Carlo and SMC estimators from Tables 1 and 2 strongly indicates the superiority of the SMC over Monte Carlo, bearing in mind that computational effort for the SMC is only about 10–100% greater than for the Monte Carlo.

## 7 Conclusion and discussion

In this paper, we presented the SMC method for pricing knock-out barrier options. General observations include the following.

• The standard error of the SMC estimator does not grow as the number of time steps increases, while the standard error of the Monte Carlo estimator can increase significantly. This is because, in the SMC, sampled asset values (particles) rejected by barrier condition are resampled from asset values between the barriers, and thus the number of particles between the barriers will not change; however, in Monte Carlo, the number of simulated paths not breaching the barrier will reduce as the number of time steps increases.

• The efficiency of the SMC versus the standard Monte Carlo improves when the probability of the asset path hitting the barrier increases (eg, the upper and lower barriers get closer or the number of time steps increases). Typically, the most significant benefit of the SMC is achieved for cases when the probability of not hitting the barrier is very small. Otherwise, its efficiency is comparable with the standard Monte Carlo.

• Implementation of the SMC requires little extra effort compared with the standard Monte Carlo method.

• Both SMC and Monte Carlo estimators are unbiased with standard errors proportional to $1/\smash{\sqrt{M}}$, where $M$ is the number of simulated asset paths for Monte Carlo and the number of particles for SMC; the proportionality coefficient for the SMC is different from that for Monte Carlo.

Further research may consider the development of SMC and Monte Carlo for the alternative solution presented in Section 2.2. Also note that it is straightforward to calculate a knock-in option as the difference between a vanilla option (ie, without a barrier) and a knock-out barrier option. However, it is not obvious how to develop the SMC estimator to calculate a knock-in option directly (ie, how to write a knock-in option price expectation via the Feynman–Kac representation formula (3.1)). This is a subject for future research. Note also that in this paper we focused on the case of one underlying asset for ease of illustration, but the SMC algorithm presented can easily be adapted for the case with many underlying assets and with additional stochastic factors, such as stochastic volatility.

## Declaration of interest

The authors report no conflicts of interest. The authors alone are responsible for the content and writing of the paper.

## References

Andersen, L., and Brotherton-Racliffe, R. (2006). Exact exotics. Risk 9(10), 85–89.

Bartoli, N., and Del Moral, P. (2001). Simulation et Algorithmes Stochastiques. Éditions Cépaduès, Toulouse.

Beaglehole, D. R., Dybvig, P. H., and Zhou, G. (1997). Going to extremes: correcting simulation bias in exotic option valuation. Financial Analysts Journal 55(1), 62–68. URL: http://bit.ly/2b5InMI (http://doi.org/fsjxm4).

Borodin, A., and Salminen, P. (1996). Handbook of Brownian Motion: Facts and Formulae. Birkhäuser, Basel (http://doi.org/bq2j).

Broadie, M., Glasserman, P., and Kou, S. (1997). A continuity correction for discrete barrier options. Mathematical Finance 7, 325–349 (http://doi.org/ddwz9f)

Carmona, R., Fouque, J.-P., and Vestal, D. (2009). Interacting particle systems for the computation of rare credit portfolio losses. Finance and Stochastics 13(4), 613–633 (http://doi.org/c2dc2n).

Carmona, R., Del Moral, P., Hu, P., and Oudjane, N. (2012). An introduction to particle methods with financial applications. In Numerical Methods in Finance, Carmona, R., Del Moral, P., Hu, P., and Oudjane, N. (eds), pp. 3–49. Springer (http://doi.org/bq2k).

Cérou, F., Del Moral, P., and Guyader, A. (2011). A nonasymptotic theorem for unnormalized Feynman–Kac particle models. Annales de l’Institut Henri Poincaré B: Probabilités et Statistiques 47(3), 629–649 (http://doi.org/ffh8bx).

Daley, D. J., and Vere-Jones, D. (2003). An Introduction to the Theory of Point Processes, Volume I: Elementary Theory and Methods, 2nd edn. Springer.

Del Moral, P. (2004). Feynman–Kac Formulae: Genealogical and Interacting Particle Systems with Applications. Probability and Its Applications. Springer (http://doi.org/bq2m).

Del Moral, P. (2013). Mean Field Simulation for Monte Carlo Integration. Monographs on Statistics and Applied Probability. Chapman and Hall/CRC.

Del Moral, P., and Patras, F. (2011). Interacting path systems for credit risk. In Credit Risk Frontiers, Brigo, D., Bielecki, T., and Patras, F. (eds), pp. 649–674. Wiley–Bloomberg Press (http://doi.org/bq2n).

Del Moral, P., Hu, P., Oudjane, N., and Rémillard, B. (2011). On the robustness of the Snell envelope. SIAM Journal on Financial Mathematics 2(1), 587–626 (http://doi.org/btbk4s).

Del Moral, P., Hu, P., and Oudjane, N. (2012a). Snell envelope with small probability criteria. Applied Mathematics and Optimization 66(3), 309–330 (http://doi.org/bq2p).

Del Moral, P., Rémillard, B., and Rubenthaler, S. (2012b). Monte Carlo approximations of American options that preserve monotonicity and convexity. In Numerical Methods in Finance, Carmona, R., Del Moral, P., Hu, P., and Oudjane, N. (eds), pp. 115–143. Springer (http://doi.org/bq2q).

Dewynne, J., and Wilmott, P. (1994). Partial to exotic. Risk 7(12), 53–57.

Giles, M. (2008a). Improved multilevel Monte Carlo convergence using the Milstein scheme. In Monte Carlo and Quasi-Monte Carlo Methods 2006, pp. 343–358. Springer (http://doi.org/d52wrs).

Giles, M. (2008b). Multilevel Monte Carlo path simulation. Operations Research 56(3), 607–617 (http://doi.org/fmrpfp).

Glasserman, P. (2004). Monte Carlo Methods in Financial Engineering. Springer.

Glasserman, P., and Staum, J. (2001). Conditioning on one-step survival for barrier option simulations. Operations Research 49(6), 923–937 (http://doi.org/ctg84z).

Gobet, E. (2009). Advanced Monte Carlo methods for barrier and related exotic options. Handbook of Numerical Analysis 15, 497–528 (http://doi.org/bsddw9).

Gobet, E., and Menozzi, S. (2010). Stopped diffusion processes: boundary corrections and overshoot. Stochastic Processes and Their Applications 120(2), 130–162 (http://doi.org/bb82kk).

He, H., Keirstead, W. P., and Rebholz, J. (1998). Double lookbacks. Mathematical Finance 8(3), 201–228 (http://doi.org/c8fm26).

Heynen, R., and Kat, H. (1994a). Crossing barriers. Risk 7(6), 46–51.

Heynen, R., and Kat, H. (1994b). Partial barrier options. Journal of Financial Engineering 3(3), 253–274.

Hull, J. C., and White, A. D. (1993). Efficient procedures for valuing European and American path-dependent options. Journal of Derivatives 1(1), 21–31 (http://doi.org/ctjp93).

Jasra, A., and Del Moral, P. (2011). Sequential Monte Carlo methods for option pricing. Stochastic Analysis and Applications 29(2), 292–316 (http://doi.org/dtk6ts).

Johannes, M. S., Polson, N. G., and Stroud, J. R. (2009). Optimal filtering of jump diffusions: extracting latent states from asset prices. Review of Financial Studies 22(7), 2759–2799 (http://doi.org/dz9mt9).

Karatzas, I., and Shreve, S. (1991). Brownian Motion and Stochastic Calculus. Springer.

Kat, H. M., and Verdonk, L. T. (1995). Tree surgery. Risk 8(2), 53–56.

Kunitomo, N., and Ikeda, M. (1992). Pricing options with curved boundaries. Mathematical Finance 2(4), 275–298 (http://doi.org/ftbn2s).

Peters, G. W., Brier, M., Shevchenko, P., and Doucet, A. (2013). Calibration and filtering for multi factor commodity models with seasonality: incorporating panel data from futures contracts. Methodology and Computing in Applied Probability 15(4), 841–874 (http://doi.org/bq2r).

Rubinstein, M., and Reiner, E. (1991). Breaking down the barriers. Risk 4(8), 28–35.

Shevchenko, P. V. (2003). Addressing the bias in Monte Carlo pricing of multi-asset options with multiple barriers through discrete sampling. The Journal of Computational Finance 6(3), 1–20 (http://doi.org/bq2s).

Shevchenko, P. V. (2011). Closed-form transition densities to price barrier options with one or two underlying assets. Technical Report EP11204, CSIRO.

Targino, R. S., Peters, G. W., and Shevchenko, P. V. (2015). Sequential Monte Carlo samplers for capital allocation under copula-dependent risk models. Insurance: Mathematics and Economics 61, 206–226 (http://doi.org/bq2t).

Only users who have a paid subscription or are part of a corporate subscription are able to print or copy content.

To access these options, along with all other subscription benefits, please contact [email protected] or view our subscription options here: http://subscriptions.risk.net/subscribe