Valuation of Barrier Options using
Sequential Monte Carlo

Pavel V. Shevchenko

CSIRO Computational Informatics, 11 Julius Ave, North Ryde, NSW 2113, Australia

School of Mathematics and Statistics UNSW, Australia

Pavel.Shevchenko@csiro.au

Pierre Del Moral

School of Mathematics and Statistics UNSW, NSW 2052, Australia

p.del-moral@unsw.edu.au

Final version 24 July 2015, 1st version 17 May 2014

Abstract

Sequential Monte Carlo (SMC) methods have successfully been used in many ap-
plications in engineering, statistics and physics. However, these are seldom used
in financial option pricing literature and practice. This paper presents SMC
method for pricing barrier options with continuous and discrete monitoring of
the barrier condition. Under the SMC method, simulated asset values rejected
due to barrier condition are re-sampled from asset samples that do not breach
the barrier condition improving the efficiency of the option price estimator; while
under the standard Monte Carlo many simulated asset paths can be rejected by
the barrier condition making it harder to estimate option price accurately. We
compare SMC with the standard Monte Carlo method and demonstrate that the
extra effort to implement SMC when compared with the standard Monte Carlo is
very little while improvement in price estimate 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
SMC estimator is smaller and does not grow with the number of time steps when
compared to the standard Monte Carlo. In this paper we demonstrate that SMC
can successfully be used for pricing barrier options. SMC can also be used for
pricing other exotic options and also for cases with many underlying assets and
additional stochastic factors such as stochastic volatility; we provide general for-
mulas and references.

**Keywords:** Sequential Monte Carlo, particle methods, Feynman-Kac representation, barrier options, Monte Carlo, option pricing
# 1 Introduction

Sequential Monte Carlo (SMC) methods (also referred to as particle methods) have successfully been used in many applications in engineering, statistics and physics for many years, especially in signal processing, state-space modelling and estimation of rare event probability. SMC method coincides with Quantum Monte Carlo method introduced as heuristic type scheme in physics by Enrico Fermi in 1948 while studying neutron diffusions. From 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 to the couple of books Del Moral (2004, 2013), and references therein. The applications of these particle methods in mathematical finance has been started recently. For instance, using the rare event interpretation of SMC, Carmona et al. (2009) and Del Moral & Patras (2011) proposed SMC algorithm for computation of the probabilities of simultaneous defaults in large credit portfolios, Targino et al. (2015) developed SMC for capital allocation problems. There are many articles utilizing SMC for estimation of stochastic volatility, jump diffusion and state-space price models, e.g. Johannes et al. (2009) and Peters et al. (2013) to name a few. The applications of SMC methods in option pricing has been started recently by the second author in the series of articles Carmona et al. (2012); Del Moral et al. (2011, 2012a,b); Jasra & Del Moral (2011). However, these methods are not widely known among option pricing practitioners and option pricing literature.

The purpose of this paper is to provide simple illustration and explanation of SMC method and its efficiency. It can be beneficial to use SMC for pricing many exotic options. For simplicity of illustration, we consider barrier options with a simple geometric Brownian motion for the underlying asset. 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). A lot of related more complex 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 options with payoff depending upon the path extrema 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 & Kat (1994b), Kunitomo & Ikeda (1992), Rubinstein & 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 & Kat (1994a) and He et al.
(1998).

In practice, however, numerical methods are used to price the barrier options for a number of reasons, for example, if the assumptions of constant volatility and drift are relaxed or payoff is too complicated. Numerical schemes such as binomial and trinomial lattices (Hull & White, 1993; Kat & Verdonk, 1995) or finite difference schemes (Dewynne & 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.

Monte Carlo (MC) simulation method is a good general pricing tool for such instruments; for review of advanced Monte Carlo methods for barrier options, see Gobet (2009). Many studies have been done to address finding the extrema of the continuously monitored assets by sampling assets at discrete dates. The standard discrete-time MC 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/\sqrt{N}$ for $N \gg 1$, where N is the number of equally spaced sampling dates (see Broadie et al. 1997, that also shows how to approximately calculate discretely monitored barrier option via continuous barrier case with some shift applied to the barrier). Also, extrapolation of the 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 & 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. The 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 simulated option payoff by the conditional probability of the path not crossing the barrier between the sampled dates); also see (Glasserman, 2004, pp. 368-370). The technique is very efficient in the case of underlying asset following 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 the time region. Closely related method of sampling underlying asset conditional on not crossing a barrier is studied in Glasserman & Staum (2001). The method of Brownian bridge simulation can also be applied in the case of multiple underlying assets as studied in Shevchenko (2003). Importance sampling and control variates methods can be applied to reduce
the variance of the barrier option price MC 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,b) has introduced a multi-level Monte Carlo path simulation method for the pricing of financial options including barrier options that improves the computational efficiency of MC path simulation by combining results using different numbers of time steps. Gobet & Menozzi (2010) developed a procedure for multidimensional stopped diffusion processes accounting for boundary correction through shifting the boundary that can be used to improve the barrier option MC estimates in the case of multi-asset and multi-barrier options with more general underlying processes.

However, the coefficient of variation of the MC estimator grows when the number of asset paths rejected by the barrier condition increases (i.e. probability of asset path to reach maturity without breaching the barrier decreases; for example, when barriers are getting closer to the asset spot). This can be improved by SMC method that re-samples 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 MC estimators are unbiased and are converging to the true value as $1/\sqrt{M}$, where M is the number of simulations (asset paths) but SMC has smaller variance.

This paper presents SMC algorithm and provides comparison between SMC and MC 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 improve the accuracy of the option price sampling estimator for a given time discretization.

The organisation 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 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 in the final section.

## 2 Model

Assume that underlying asset $S_t$ follows risk neutral process

$$dS_t = S_t \mu dt + S_t \sigma dW_t, \quad (1)$$
where $\mu = r - q$ is the drift, $r$ is risk free interest rate, $q$ is continuous dividend rate (it corresponds to the foreign interest rate if $S_t$ is exchange rate or continuous dividends if $S_t$ is stock), $\sigma$ is volatility and $W_t$ is the standard Brownian motion. The interest rate can be function of time, and drift and volatility can be functions of time and underlying asset. In this paper, we do not consider time discretization errors; for simplicity, hereafter, we assume that model parameters are piece-wise constant functions of time.

## 2.1 Pricing Barrier Option

The today's fair price of continuously monitored knock-out barrier option with the lower barrier $L_t$ and upper barrier $U_t$ can be calculated as expectation with respect to risk neutral process (1), given information today at $t_0 = 0$ (i.e. conditional on $S_0 = s_0$)

$$Q_C = B_{0,T} E \left( h(S_T) 1_{A_t}(S_t)_{t \in [0,T]} \right), \quad B_{0,T} = e^{-\int_0^T r(\tau)d\tau}, \qquad (2)$$

where $B_{0,T}$ is the discounting factor from maturity $T$ to $t_0 = 0$; $1_A(x)$ is indicator function equals 1 if $x \in A$ and 0 otherwise; $h(x)$ is payoff function, i.e. $h(x) = \max(x - K, 0)$ for call option and $h(x) = \max(K - x, 0)$ for put option, where $K$ is strike price; and $A_t = (L_t, U_t)$. All standard barrier structures such as lower barrier only, upper barrier only or several window barriers can be obtained by setting $L_t = 0$ or $U_t = \infty$ for corresponding time periods.

Assume that drift, volatility and barriers are piecewise constant functions of time for time discretization $0 = t_0 < t_1 < \dots < t_N = T$. Denote corresponding asset values as $S_0, S_1, \dots, S_N$; the lower and upper barriers as $L_1, \dots, L_N$ and $U_1, \dots, U_N$ respectively; and drift and volatility as $\mu_1, \dots, \mu_N$ and $\sigma_1, \dots, \sigma_N$. That is, $L_1$ is the lower barrier for time period $[t_0, t_1]$; $L_2$ is for $[t_1, t_2]$, etc. and similar 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$ or $U_n = \infty$ respectively.

Denote the transition density from $S_n$ to $S_{n+1}$ as $f(S_{n+1}|S_n)$ which is just a lognormal density in the case of process (1) with solution

$$S_n = S_{n-1} \exp \left( \left( \mu_n - \frac{1}{2}\sigma_n^2 \right) \delta t_n + \sigma_n \sqrt{\delta t_n} Z_n \right), \quad n = 1, \dots, N, \qquad (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 barrier monitored at $t_0, t_1, \dots, t_N$ (discretely monitored barrier), the option price (2) simplifies to

$$Q_D = B_{0,T} E \left( h(S_N) \prod_{n=1}^{N} 1_{(L_n, U_n)}(S_n) \right). \quad (4)$$

It is a biased estimate of continuously monitored barrier option $Q_C$ such that $Q_D \to Q_C$ for $\delta t_n \to 0$; see Broadie et al. (1997) that also shows how to approximately calculate discretely monitored barrier option via continuous barrier option price with some shift applied to the barrier in the case of one-dimension Brownian motion (for high dimensional case and more general processes, see Gobet & Menozzi 2010).

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

$$Q_C = B_{0,T} \int_{L_1}^{U_1} ds_1 f(s_1|s_0) g(s_0, s_1) \dots \int_{L_N}^{U_N} ds_N f(s_N|s_{N-1}) g(s_{N-1}, s_N) h(s_N), \quad (5)$$

where $g(S_{n-1}, S_n)$ is 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 \left( -2 \frac{\ln(S_n/B_n) \ln(S_{n-1}/B_n)}{\sigma_n^2 \delta t_n} \right); \quad (6)$$

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

$$g(S_{n-1}, S_n) = 1 - \sum_{m=1}^{\infty} [R_n(\alpha_n m - \gamma_n, x_n) + R_n(-\alpha_n m + \beta_n, x_n)] \\ + \sum_{m=1}^{\infty} [R_n(\alpha_n m, x_n) + R_n(-\alpha_n m, x_n)], \quad (7)$$

where

$$x_n = \ln \frac{S_n}{S_{n-1}}, \alpha_n = 2 \ln \frac{U_n}{L_n}, \beta_n = 2 \ln \frac{U_n}{S_{n-1}}, \gamma_n = 2 \ln \frac{S_{n-1}}{L_n}, R_n(z, x) = \exp \left( - \frac{z(z-2x)}{2\sigma_n^2 \delta t_n} \right).$$

Typically few terms in the above summations are enough to obtain a good accuracy (in the actual implementation the number of terms can be adaptive to achieve the required accuracy; the smaller time step $\delta t_n$ the less number of terms is needed). Formulas (6) and (7) can easily be obtained from the well known distribution of maximum and minimum of a Brownian motion (see e.g. Borodin & Salminen, 1996; Karatzas & Shreve, 1991); also can be found in Shevchenko (2011).
The integral (5) can be rewritten as

$$
\begin{align}
Q_C &= B_{0,T} \int_0^\infty ds_1 f(s_1|s_0) g(s_0, s_1) 1_{(L_1,U_1)}(s_1) \cdots \nonumber \\
    &\qquad \int_0^\infty ds_N f(s_N|s_{N-1}) g(s_{N-1}, s_N) h(s_N) 1_{(L_U,U_N)}(s_N) \nonumber \\
    &= B_{0,T} \times E \left( h(S_N) \prod_{n=1}^N \left( 1_{(L_n,U_n)}(S_n) g(S_{n-1}, S_n) \right) \right). \tag{8}
\end{align}
$$

Alternative expression for the barrier option that might provide more efficient nu-
merical estimate is presented by formula (11) in the next section. It is not analysed in
this paper and subject of further study.

2.2 Alternative Solution for Barrier Option

The integral for barrier option price (8) can also be rewritten in terms of the Markov
chain $\hat{S}_n$, starting at $\hat{S}_0 = S_0$, with elementary transitions

$$
\mathrm{Pr} \left( \hat{S}_n \in ds_n \mid \hat{S}_{n-1} = s_{n-1} \right) := \frac{\mathrm{Pr} \left( S_n \in ds_n \mid S_{n-1} = s_{n-1} \right) 1_{(L_n, U_n)}(s_n)}{\mathrm{Pr} \left( S_n \in (L_n, U_n) \mid S_{n-1} = s_{n-1} \right)}. \quad (9)
$$

We readily check that

$$
\hat{S}_n = \hat{S}_{n-1} \exp \left( a_n + b_n \hat{Z}_n \right) \quad (10)
$$

with

$$
a_n := (\mu_n - \frac{1}{2}\sigma_n^2)\delta t_n \quad \text{and} \quad b_n := \sigma_n\sqrt{\delta t_n}.
$$

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

$$
A_n(\hat{S}_{n-1}) := \left[ \ln \left( \frac{L_n}{\hat{S}_{n-1}} \right) - a_n \right] / b_n \quad \text{and} \quad B_n(\hat{S}_{n-1}) := \left[ \ln \left( \frac{U_n}{\hat{S}_{n-1}} \right) - a_n \right] / b_n.
$$

Let $\Phi(x) := \int_{-\infty}^{x} \frac{1}{\sqrt{2\pi}} e^{-y^2/2} dy$ be the standard Normal (Gaussian) distribution function and its inverse function is $\Phi^{-1}(\cdot)$. In this notation, we have that

$$
\begin{align*}
\mathrm{Pr} (S_n \in (L_n, U_n) | S_{n-1} = s_{n-1}) &= \mathrm{Pr} (Z_n \in (A_n(s_{n-1}), B_n(s_{n-1})) | S_{n-1} = s_{n-1}) \\
&= \Phi(B_n(s_{n-1})) - \Phi(A_n(s_{n-1})).
\end{align*}
$$

We can also simulate the transition $\hat{S}_{n-1} \rightsquigarrow \hat{S}_n$ by sampling a uniform random variable $\mathcal{U}_n$ by taking in (10)

$$
\hat{Z}_n := \Phi^{-1} \left[ \Phi \left( A_n(\hat{S}_{n-1}) \right) + \mathcal{U}_n \left( \Phi \left( B_n(\hat{S}_{n-1}) \right) - \Phi \left( A_n(\hat{S}_{n-1}) \right) \right) \right].
$$
If we set

$$ \varphi_{k-1}(s_{k-1}) := \Pr(S_k \in (L_k, U_k) | S_{k-1} = s_{k-1}) = \Phi(B_k(s_{k-1})) - \Phi(A_k(s_{k-1})), $$

then we have that

$$ \left\{ \prod_{k=1}^{n} 1_{(L_k, U_k)}(s_k) \right\} \left\{ \prod_{k=1}^{n} \Pr(S_k \in ds_k \mid S_{k-1} = s_{k-1}) \right\} \\ = \left\{ \prod_{k=1}^{n} \Pr(S_k \in (L_k, U_k) \mid S_{k-1} = s_{k-1}) \right\} \left\{ \prod_{k=1}^{n} \Pr(\hat{S}_k \in ds_k \mid \hat{S}_{k-1} = s_{k-1}) \right\} \\ = \left\{ \prod_{k=1}^{n} \varphi_{k-1}(s_{k-1}) \right\} \left\{ \prod_{k=1}^{n} \Pr(\hat{S}_k \in ds_k \mid \hat{S}_{k-1} = s_{k-1}) \right\} $$

from which we conclude that

$$ Q_C = B_{0,T} \times E \left( h(\hat{S}_N) \prod_{n=1}^{N} \hat{G}_{n-1}(\hat{S}_{n-1}, \hat{S}_n) \right) \quad (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). \quad (12) $$

Explicitly, the option price integral becomes

$$ \begin{aligned} Q_C &= B_{0,T} \int_0^1 dw_1 (\Phi(\tilde{U}_1) - \Phi(\tilde{L}_1)) g(s_0, s_1) \cdots \int_0^1 dw_N (\Phi(\tilde{U}_N) - \Phi(\tilde{L}_N)) g(s_{N-1}, s_N) h(s_N) \\ &= B_{0,T} \int_0^1 \cdots \int_0^1 dw_1 \cdots dw_N h(s_N) \prod_{n=1}^N (\Phi(\tilde{U}_n) - \Phi(\tilde{L}_n)) g(s_{N-1}, s_N), \end{aligned} \quad (13) $$

where

$$ \tilde{U}_n = (\ln(U_n/s_{n-1}) - (\mu_n - \frac{1}{2}\sigma_n^2)\delta t_n)/(\sigma_n\sqrt{\delta t_n}), $$

$$ \tilde{L}_n = (\ln(L_n/s_{n-1}) - (\mu_n - \frac{1}{2}\sigma_n^2)\delta t_n)/(\sigma_n\sqrt{\delta t_n}), $$

$$ z_n = \Phi^{-1}[\Phi(\tilde{L}_n) + w_n(\Phi(\tilde{U}_n) - \Phi(\tilde{L}_n))], $$

$$ s_n = s_{n-1} \exp((\mu_n - \frac{1}{2}\sigma_n^2)\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 more efficient numerical estimate but it is not analysed in this paper.
# 3 Feynman-Kac representations

In this section, we provide the basic option price formulas under Feynman-Kac representation underlying SMC method; for detailed introduction of 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 continuously monitored barrier (8) can be written as

$$Q_C = B_{0,T} \times E \left( H(X_N) \prod_{n=0}^{N-1} G_n(X_n) \right) \quad (14)$$

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}) \times 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 (14) is the Feynman-Kac formula for discrete time models (see Carmona et al. 2012) which is used to develop SMC option price estimator.

In this notation, the discretely monitored barrier option expectation (4) also takes the following form

$$Q_D = B_{0,T} \times E \left( H(X_N) \prod_{n=0}^{N-1} \tilde{G}_n(X_n) \right) \quad (15)$$

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 formulae for barrier option expectation presented in Section 2.2 by formula (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 modified underlying asset process $\hat{S}_n$ given by (10), then we can rewrite the formula (11) as follows

$$Q_C = B_{0,T} \times E \left( H(\hat{X}_N) \prod_{n=0}^{N-1} \hat{G}_n(\hat{X}_n) \right) \quad (16)$$

with the potential function $\hat{G}_n$ defined in (12). We observe that the above expression has exactly the same form as (14) by replacing $(X_n, G_n)$ by $(\hat{X}_n, \hat{G}_n)$.

Once the option price expectation is written in 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 formulae 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 formulae describing the normalizing constants in terms of normalized Feynman-Kac measures. These mathematical objects are essential to define and to analyze particle approximation models. For instance, the particle approximation of normalizing constants are defined mimicking the multiplicative formula discussed above, by 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 the monographs (Del Moral, 2004, Section 2.7.1) and (Del Moral, 2013, Section 3.2.2).

Firstly, we observe that (14) can be written in the following form

$$Q_C = B_{0,T} \gamma_N(H) = B_{0,T} \gamma_N(1) \eta_N(H) \quad (17)$$

with the Feynman-Kac unnormalized $\gamma_N$ and normalized $\eta_N$ measures given for any function $\varphi$ by the formulae

$$\gamma_N(\varphi) = E \left( \varphi(X_N) \prod_{n=0}^{N-1} G_n(X_n) \right) \quad \text{and} \quad \eta_N(\varphi) = \gamma_N(\varphi)/\gamma_N(1). \quad (18)$$

Notice that the sequence of non negative measures $(\gamma_n)_{n \ge 0}$ satisfies for any bounded measurable function $\varphi$ the recursive linear equation

$$\gamma_n(\varphi) = \gamma_{n-1}(Q_n(\varphi)) \quad (19)$$
with the integral operator

$$Q_n(\varphi)(x) = G_{n-1}(x) K_n(\varphi)(x), \qquad (20)$$

where

$$K_n(\varphi)(x) = E(\varphi(X_n) | X_{n-1} = x) = \int K_n(x, dy) \varphi(y) \quad (21)$$

and $K_n(X_{n-1}, dx) := \text{Pr}(X_n \in dx | X_{n-1})$ is the Markov transition in the chain $X_n$. We prove this claim using the fact that

$$\begin{aligned}
\gamma_n(\varphi) &= E\left(E\left(\varphi(X_n) \prod_{p=0}^{n-1} G_p(X_p) \mid (X_0, \dots, X_{n-1})\right)\right) \\
&= E\left(E\left(\varphi(X_n) \mid (X_0, \dots, X_{n-1})\right) \prod_{p=0}^{n-1} G_p(X_p)\right) \\
&= E\left(E\left(\varphi(X_n) \mid X_{n-1}\right) \prod_{p=0}^{n-1} G_p(X_p)\right) \\
&= E\left(G_{n-1}(X_{n-1}) K_n(\varphi)(X_{n-1}) \prod_{p=0}^{n-2} G_p(X_p)\right) \\
&= \gamma_{n-1}(G_{n-1} \times K_n(\varphi)). \qquad (22)
\end{aligned}$$

By construction, we also have that

$$\begin{aligned}
\gamma_N(1) &= E\left(\prod_{n=0}^{N-1} G_n(X_n)\right) \\
&= E\left(G_{N-1}(X_{N-1}) \times \prod_{n=0}^{N-2} G_n(X_n)\right) = \gamma_{N-1}(G_{N-1}). \qquad (23)
\end{aligned}$$

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}) \qquad (24)$$

from which we conclude that

$$\gamma_N(1) = \prod_{0 \le n < N} \eta_n(G_n) \qquad (25)$$

and therefore

$$Q_C = B_{0,T} \times \left[ \prod_{0 \le n < N} \eta_n(G_n) \right] \times \eta_N(H), \qquad (26)$$

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 MC and SMC estimators and corresponding algorithms to calculate option price in the case of continuously and discretely monitored barrier conditions.

## 4.1 Standard Monte Carlo

Using process (3), simulate independent asset path realizations $S^{(m)} = (S_1^{(m)}, \dots, S_N^{(m)})$, $m = 1, \dots, M$. Then, the unbiased estimator for continuously monitored barrier option price integral (8) is a standard average of option price payoff realisations over simulated paths

$$ 
\begin{aligned}
\hat{Q}_C^{MC} &= B_{0,T} \frac{1}{M} \sum_{m=1}^{M} \left( h(S_N^{(m)}) \prod_{n=1}^{N} \{g(S_{n-1}^{(m)}, S_n^{(m)}) 1_{[L_n, U_n]}(S_n^{(m)})\} \right) \\
&= B_{0,T} \frac{1}{M} \sum_{m=1}^{M} \left( H(X_N^{(m)}) \prod_{n=0}^{N-1} G_n(X_n^{(m)}) \right)
\end{aligned}
\quad (27) $$

with $X_n^{(m)} = (S_n^{(m)}, S_{n+1}^{(m)})$ and the unbiased estimator for discretely monitored barrier option (4) is

$$ \hat{Q}_D^{MC} = B_{0,T} \frac{1}{M} \sum_{m=1}^{M} \left( H(X_N^{(m)}) \prod_{n=0}^{N-1} \tilde{G}_n(X_n^{(m)}) \right) \quad (28) $$

with $X_n^{(m)} = (S_n^{(m)}, S_{n+1}^{(m)})$.

## 4.2 Sequential Monte Carlo

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

### • Initial step

1. (proposition step) 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 (3); these are referred to as M (transition type) particles. Set

$$ G_0(X_0^{(m)}) = 1_{(L_1, U_1)}(S_1^{(m)}) \times g(S_0^{(m)}, S_1^{(m)}) $$

for each $1 \le m \le M$.
2. (acceptance-rejection step) Sample M random and [0, 1]-valued uniform variables $U_0^{(m)}$. The rejected transition type particles $X_0^{(m)}$ are those for which $G_0(X_0^{(m)}) < U_0^{(m)}$. The particles $X_0^{(m)}$ for which $U_0^{(m)} \le G_0(X_0^{(m)})$ are accepted. Notice that a transition type particle $X_0^{(m)}$ s.t. $S_1^{(m)} \notin (L_1, U_1)$ is instantly rejected (since its weight $G_0(X_0^{(m)}) = 0$ is null); and a transition type particle $X_0^{(m)}$ s.t. $S_1^{(m)} \in (L_1, U_1)$ is rejected with a probability $1 - G_0(X_0^{(m)})$.

3. (recycling-selection step) Resample each rejected transition type particle $X_0^{(m)}$ by resampling its $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)}), \quad (29)$$

where $\delta(y - y_0)$ is a point mass function centered at $y_0$ (i.e. 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 one). In other words, when a transition type particle, say $X_0^{(r)}$, is rejected for some index $r$ we replace it by one of the particle $X_0^{(m)}$ randomly chosen w.r.t. 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 (29) can be accomplished by Algorithm 4.1 in Section 4.3.

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

$$\tilde{X}_0^{(m)} = (S_0^{(m)}, \tilde{S}_1^{(m)}) \quad m = 1, \dots, M.$$

**Remark 4.1** By definition of (29) we notice that transition type particles $X_0^{(m)}$ s.t. $S_1^{(m)} \notin (L_1, U_1)$ have a null weight. Therefore, they cannot be selected in replacement of the rejected ones. Moreover, the transition type particles $X_0^{(m)}$ s.t. $S_1^{(m)} \in (L_1, U_1)$ with a large probability $G_0(X_0^{(m)})$ of non hitting the barrier within $[t_0, t_1]$ are more likely to be selected (in replacement of the rejected ones).

• Step 0 $\sim 1$

a) (proposition) For the 2nd 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 $\tilde{S}_1^{(m)}$ of the selected transitions $\tilde{X}_0^{(m)}$ at the previous step, using the process evolution (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_2^{(m)}) \times g(\tilde{S}_1^{(m)}, S_2^{(m)})$$

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

b) (acceptance-rejection) Sample $M$ random and $[0, 1]$-valued uniform variables $U_1^{(m)}$. The rejected transition type particles $X_1^{(m)}$ are those for which $G_1(X_1^{(m)}) < U_1^{(m)}$ and the particles $X_1^{(m)}$ for which $U_1^{(m)} \le G_1(X_1^{(m)})$ are accepted. That is a transition type particle $X_1^{(m)}$ s.t. $S_2^{(m)} \notin (L_2, U_2)$ is instantly rejected and $X_1^{(m)}$ s.t. $S_2^{(m)} \in (L_2, U_2)$ is rejected with a probability $1 - G_1(X_1^{(m)})$.

c) (recycling-selection) Resample each rejected transition type particles $X_1^{(m)}$ by resampling its $S_2^{(m)}$ component from the discrete distribution with the 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)}) \quad (30)$$

using e.g. efficient and simple Algorithm 4.1 in Section 4.3.

At the end of the acceptance-rejection-recycling scheme, we have $M$ (transition-type) particles denoted as $\tilde{X}_1^{(m)} = (\tilde{S}_1^{(m)}, \tilde{S}_2^{(m)})$, $1 \le m \le M$. A remark similar to Remark 4.1 is also applied here: transition type particles $X_1^{(m)}$ s.t. $S_2^{(m)} \notin (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 $X_1^{(m)}$ s.t. $S_2^{(m)} \in (L_2, U_2)$ with a large probability $g(\tilde{S}_1^{(m)}, S_2^{(m)})$ of non hitting the barrier within $[t_0, t_1]$ are more likely to be selected (in replacement of the rejected ones).

• Repeat steps a) to c) in **Step 0** $\rightsquigarrow 1$ for time steps $[t_2, t_3], \dots, [t_{N-1}, t_N]$.

Calculate the final unbiased option price estimator as

$$\hat{Q}_C^{SMC} = B_{0,T} \times \left[ \prod_{n=0}^{N-1} \frac{1}{M} \sum_{m=1}^{M} G_n(X_n^{(m)}) \right] \times \frac{1}{M} \sum_{m=1}^{M} H(X_N^{(m)}). \quad (31)$$

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 formula (26).
Note that $H(X_N^{(m)}) = h(\tilde{S}_N^{(m)})$, i.e. payoff at maturity is calculated using particles $\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_D$ defined in (15) is given by

$$ \hat{Q}_D^{SMC} = B_{0,T} \times \left[ \prod_{n=0}^{N-1} \frac{1}{M} \sum_{m=1}^{M} \tilde{G}_n(X_n^{(m)}) \right] \times \frac{1}{M} \sum_{m=1}^{M} H(X_N^{(m)}), \quad (32) $$

where $(X_n^{(m)})_{0 \le n \le N}$, $1 \le m \le M$ is obtained by the above algorithm with potential functions $(G_n)_{0 \le n \le N}$ replaced by the indicator potential functions $(\tilde{G}_n)_{0 \le n \le N}$.

In both cases, it may happen that all the particles 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 (16) for alternative option price expression (11) presented in Section 2.2. In this context, an unbiased estimator of $Q_C$ is given by

$$ \hat{Q}_C^{SMC} = B_{0,T} \times \left[ \prod_{n=0}^{N-1} \frac{1}{M} \sum_{m=1}^{M} \hat{G}_n(\hat{X}_n^{(m)}) \right] \times \frac{1}{M} \sum_{m=1}^{M} H(\hat{X}_N^{(m)}), \quad (33) $$

where $(\hat{X}_n^{(m)})_{0 \le n \le N}$, $1 \le m \le M$, is obtained by the above algorithm for $(X_n^{(m)})_{0 \le n \le N}$ with potential functions $(G_n)_{0 \le n \le N}$ replaced by the potential functions $(\hat{G}_n)_{0 \le n \le N}$ and process for $S_n$ is replaced by process $\hat{S}_n$ as described in Section 2.2.

**Remark 4.2** As we mentioned in the introduction of Section 3.2, the particle estimate in (31) is defined as in (26) by replacing the normalized Feynman-Kac measures $\eta_n$ by the particle empirical approximations. Formulae (32), and respectively (33), follow the same line of arguments based on the Feynman-Kac model (15), and respectively (16).

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 $S_1^{(4)}$ is rejected and resampled (moved to position $S_1^{(1)}$), particle $S_1^{(6)}$ is rejected and moved to position $S_1^{(3)}$. Then two particles located at $S_1^{(3)}$ will generate two particles at $t_2$, two particles located at $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 $S_1^{(1)}$, $S_1^{(2)}$, $S_1^{(3)}$, $S_1^{(5)}$ are also rejected in the case of continuously monitored barrier.
<figure><img src="image_1.png" /><figcaption>Figure 1: Illustration of Sequential Monte Carlo algorithm to calculate barrier option with the lower barrier at level L. Particle $S_1^{(4)}$ is rejected and moved to position $S_1^{(1)}$ (resampled), particle $S_1^{(6)}$ is rejected and moved to position $S_1^{(3)}$, etc. Note that it is possible that $S_1^{(1)}$, $S_1^{(2)}$, $S_1^{(3)}$, $S_1^{(5)}$ are also rejected in the case of continuously monitored barrier.</figcaption></figure>

## 4.3 Sampling from discrete distribution

For the benefit of the reader, in this section we present efficient and simple algorithm for sampling of the rejected particles from the discrete density required during *recycle-selection* step of SMC algorithm described in previous section, i.e. sampling from discrete densities (29) and (30).

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

$$f(x) = \sum_{m=1}^{M} p_m \delta(x - x_m) \qquad (34)$$

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 (34) and $X = F^{-1}(U)$ is a sample from $F(x)$ if $U$ is from uniform (0,1) distribution. It is important
to use computationally efficient method for sampling of R variables. If the order of samples is not important (as in the case of recycling-selection steps of SMC algorithm in Section 4.2) then, for example, one can sample (R + 1) independent exponential random variables ($\mathcal{E}_r$)$_{1\le r\le(R+1)}$ with unit parameter and set

$$ \mathcal{T}_r = \sum_{1 \le s \le r} \mathcal{E}_s \quad \text{and} \quad \mathcal{V}_r = \mathcal{T}_r / \mathcal{T}_{R+1}, \quad r = 1, 2, \dots, R+1. \qquad (35) $$

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 Poisson process, see e.g. (Bartoli & Del Moral, 2001, Example 3.6.9 and Section 2.6.2) or (Daley & Vere-Jones, 2003, Exercise 2.1.2). Then sampling of ($Y^{(r)}$)$_{1\le r\le R}$, by calculating $Y^{(r)} = F^{-1}(\mathcal{V}_r)$, can be accomplished using the following synthetic pseudo code.

### Algorithm 4.1

1. $k = 1$ and $r = 1$

2. While $r \le R$

• While $\mathcal{V}_r < p_1 + \dots + p_k$
  - $Y^{(r)} = x_k$
  - $r = r + 1$

• End while

• $k = k + 1$

3. End while

The computational cost of this sampling scheme is linear with respect to R. In particular, to simulate from the probability density (29) set $p_m = \frac{G_0(X_0^{(m)})}{\sum_{k=1}^{M} G_0(X_0^{(k)})}$ and $x_m = S_1^{(m)}$, and to simulate from the discrete distribution (30) set $p_m = \frac{G_1(X_1^{(m)})}{\sum_{k=1}^{M} G_1(X_1^{(k)})}$ and $x_m = S_2^{(m)}$ in (34).
## 4.4 Unbiasedness properties

SMC estimator for option price (31) can be written as

$$ \hat{Q}^{SMC} = B_{0,T} \times \gamma_{N}^{M}(1) \times \eta_{N}^{M}(H) \quad (36) $$

with the empirical measures $\eta_N^M$ given by

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

and normalizing constants

$$ \gamma_{N}^{M}(1) = \prod_{p=0}^{N-1} \frac{1}{M} \sum_{m=1}^{M} G_{p}(X_{p}^{(m)}) = \prod_{p=0}^{N-1} \eta_{p}^{M}(G_{p}). \quad (38) $$

In this notation, the M-particle approximations of the Feynman-Kac measures $\gamma_N$ for any function $\varphi$ are given by

$$ \gamma_N^M(\varphi) := \gamma_N^M(1) \times \eta_N^M(\varphi) \quad \Rightarrow \quad \hat{Q}^{SMC} = B_{0,T} \times \gamma_N^M(H). \quad (39) $$

Here, $\eta_N^M$ and $\gamma_N^M$ are particle empirical approximations of Feynman-Kac measures $\eta_N$ and $\gamma_N$ in the option price formula (26).

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

$$ \sup_{\|\varphi\| \le 1} \| E(\eta_N^M(\varphi)) - \eta_N(\varphi) \| \le c(N)/M \quad (40) $$

for some finite positive constant $c(N)$ whose values only depend on the time horizon $N$. That is, $\eta_N^M(\varphi)$ converges to $\eta_N(\varphi)$ as $M$ increases. The unnormalized particle measures $\gamma_N^M$ in (31), (32), and (33) are unbiased. On the other hand, the empirical measures $\eta_N^M(\varphi)$ can be expressed in terms of the ratio of two unnormalized quantities $\gamma_N^M(\varphi)$ and $\gamma_N^M(1)$. Taking into considerations the fluctuation of these unnormalized particle models, the estimate of the bias (40) is obtained using an elementary Taylor type expansion at the first order of this ratio.
To prove that $\gamma_N^M(H)$ is unbiased, i.e. $\hat{Q}^{SMC}$ is unbiased, recall that the particles evolve sequentially using a selection and a mutation transition. Thus we have the conditional expectation formula

$$ 
\begin{aligned}
& E \left( \eta_N^M(H) \mid \left( X_0^{(m)}, \dots, X_{N-1}^{(m)} \right)_{1 \le m \le M} \right) \\
&= E \left( H \left( X_N^{(1)} \right) \mid \left( X_0^{(m)}, \dots, X_{N-1}^{(m)} \right)_{1 \le m \le M} \right) \\
&= \sum_{1 \le m \le M} \frac{G_{N-1}(X_{N-1}^{(m)})}{\sum_{1 \le k \le M} G_{N-1}(X_{N-1}^{(k)})} K_N(H)(X_{N-1}^{(m)}), \quad (41)
\end{aligned}
$$

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

$$ 
\begin{aligned}
& E \left( \gamma_N^M(H) \mid \left( X_0^{(m)}, \dots, X_{N-1}^{(m)} \right)_{1 \le m \le M} \right) \\
&= \left[ \prod_{p=0}^{N-1} \eta_p^M(G_p) \right] \frac{1}{N} \sum_{1 \le m \le M} \frac{G_{N-1}(X_{N-1}^{(m)})}{\frac{1}{N} \sum_{1 \le k \le M} G_{N-1}(X_{N-1}^{(k)})} K_N(H)(X_{N-1}^{(m)}) \\
&= \left[ \prod_{p=0}^{N-2} \eta_p^M(G_p) \right] \times \eta_{N-1}^M(\mathcal{Q}_N(H)) \quad (42)
\end{aligned}
$$

with the one step Feynman-Kac semigroup $\mathcal{Q}_N$ introduced in (20). That is

$$ E \left( \gamma_N^M(H) \mid \left( X_0^{(m)}, \dots, X_{N-1}^{(m)} \right)_{1 \le m \le M} \right) = \gamma_{N-1}^M(\mathcal{Q}_N(H)) \quad (43) $$

and therefore

$$ E(\gamma_N^M(H)) = E(\gamma_{N-1}^M(\mathcal{Q}_N(H))). \quad (44) $$

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

$$ \gamma_0^M = \eta_0^M \Rightarrow E(\gamma_0^M(\varphi)) = E(\eta_0^M(\varphi)) = \eta_0(\varphi) = \gamma_0(\varphi) $$

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

$$ E(\gamma_n^M(\varphi)) = \gamma_n(\varphi) \quad (45) $$
at some rank $n$, for any $M \ge 1$ and any $\varphi$. In this case, arguing as above we have

$$E(\gamma_{n+1}^M(\varphi)) = E(\gamma_n^M(Q_{n+1}(\varphi))). \quad (46)$$

Under the induction hypothesis, this implies that

$$E(\gamma_{n+1}^M(\varphi)) = \gamma_n(Q_{n+1}(\varphi)) = \gamma_{n+1}(\varphi). \quad (47)$$

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

# 5 Importance sampling models

The Feynman-Kac representation formulae (14) and their particle interpretations discussed in Section 4.2 are far from being unique. For instance, using (8), for any non negative probability density functions $\bar{f}(s_n|s_{n-1})$, we also have that

$$Q = B_{0,T} \int_0^\infty ds_1 \bar{f}(s_1|s_0) \bar{g}(s_0, s_1) 1_{(L_1, U_1)}(s_1) \cdots \\ \int_0^\infty ds_N \bar{f}(s_N|s_{N-1}) \bar{g}(s_{N-1}, s_N) h(s_N) 1_{(L_U, U_N)}(s_N) \quad (48)$$

with the potential functions

$$\bar{g}(\bar{S}_{n-1}, \bar{S}_n) = g(\bar{S}_{n-1}, \bar{S}_n) \times \frac{f(s_n|s_{n-1})}{f(s_n|s_{n-1})}. \quad (49)$$

This yields the Feynman-Kac representation

$$Q = B_{0,T} \times E \left( h(\bar{S}_N) \prod_{n=1}^{N} \bar{G}_n(\bar{S}_{n-1}, \bar{S}_n) \right) \quad (50)$$

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) \quad (51)$$
and the Markov chain $(\overline{S}_n)_{n \ge 0}$, with

$$ \Pr(\overline{S}_n \in ds_n | \overline{S}_{n-1}) = \overline{f}(s_n | \overline{S}_{n-1}) ds_n. \qquad (52) $$

The importance sampling formula (50) is rather well known. The corresponding M-particle consist with M particles evolving, between the selection times, as independent copies of the twisted Markov chain model $\overline{S}_n$; and the selection/recycling procedure favors transitions $\overline{S}_{n-1} \rightsquigarrow \overline{S}_n$ that increase density ratio $f(\overline{S}_n | \overline{S}_{n-1}) / \overline{f}(\overline{S}_n | \overline{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 $(h_n)_{0 \le n \le N}$ with $h_N = h$, using the fact that

$$ h(\overline{S}_N) = \frac{h_N(\overline{S}_N)}{h_{N-1}(\overline{S}_{N-1})} \times \frac{h_{N-1}(\overline{S}_{N-1})}{h_{N-2}(\overline{S}_{N-2})} \times \dots \times \frac{h_1(\overline{S}_1)}{h_0(\overline{S}_0)} \times h_0(\overline{S}_0), \qquad (53) $$

we also have that

$$ \begin{aligned} Q_0 &= B_{0,T} \times h_0(\bar{s}_0) \times E \left( \prod_{n=1}^{N} \left( \frac{h_n(\bar{S}_n)}{h_{n-1}(\bar{S}_{n-1})} \right) \bar{G}_n(\bar{S}_{n-1}, \bar{S}_n) \right) \\ &= B_{0,T} \times h_0(\bar{s}_0) \times E \left( \prod_{n=1}^{N} \tilde{G}_n(\bar{S}_{n-1}, \bar{S}_n) \right) \end{aligned} \qquad (54) $$

with

$$ \tilde{G}_n(\overline{S}_{n-1}, \overline{S}_n) = \overline{G}_n(\overline{S}_{n-1}, \overline{S}_n) \times \frac{h_n(\overline{S}_n)}{h_{n-1}(\overline{S}_{n-1})}. \qquad (55) $$

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

$$ h_N(x) = h(x) = \max(K - x, 0) \quad \text{and} \quad \forall n < N \quad h_n(x) = h(x) + 1 \qquad (56) $$

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

The choice of the potential functions (49) allows to choose the reference Markov chain to explore randomly the state space during the mutation transitions. The importance sampling Feynman-Kac model (54) is less intrusive. More precisely, without changing
the reference Markov chain, the choice of the potential functions (55) allows to favor
transitions that increase sequentially the payoff function. The importance sampling
models (49) and (55) can be combined in an obvious way so that 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 market data: spot $S_0 = 100$, interest rate $r = 0.1$, volatility $\sigma = 0.3$ and zero dividends $q = 0$. Exact closed form solution, SMC and standard MC estimators, standard errors of the estimators, and estimator efficiencies for this option are presented in Tables 1 and 2 and Figures 2 and 3 for continuously and discretely monitored barrier cases. We perform $M = 100,000$ simulations for MC estimators and $M = 100,000$ particles for SMC estimators that are repeated 50 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 convergence of discretely monitored barrier to the continuous case and not to address time discretization errors, but to illustrate and explain the behavior of SMC that improves the accuracy of option price sampling estimator for given time discretization. In the case of real barrier option, the time discretization will be dictated by the stochastic process, window barrier structure, barrier monitoring type (e.g. continuous, daily) and market data term-structures.

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

* Computing CPU time $t_{cpu}$ is proportional to the number of simulations $M$ in MC method (or the number of particles $M$ in SMC).

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

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

$$s^2 = \alpha/t_{cpu}, \qquad (57)$$

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

$$\kappa = \alpha_{MC}/\alpha_{SMC}. \qquad (58)$$

Interpretation of $\kappa$ is straightforward; if computing time for SMC estimator is $t_{SMC}$, then the computing time for MC estimator to achieve the same accuracy as SMC estimator is $\kappa \times t_{SMC}$, i.e. $\kappa > 1$ indicates that SMC is faster than MC and $\kappa < 1$ otherwise.

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

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

Other numerical experiments not reported here show that efficiency of SMC over MC improves when barriers become closer, i.e. probability for asset path to hit the barrier increases; it is also easy to see from results in Table 2. If probability of asset path not hitting the barrier is large then performance of SMC is about the same or slightly worse than MC. Note that our implementation does not include any standard variance reduction techniques such as antithetics, importance sampling and control variates or any parallel/vector computations. The algorithm was implemented using Fortran 90 and executed on a standard laptop (Windows 7, Intel(R) i7-2640M CPU @ 2.8GHz, RAM 4 GB). While computing time is somewhat subjective (i.e. depends on specifics of our implementation), the ratio of standard errors (or ratio of squared standard errors) of MC and SMC estimators from Tables 1 and 2 strongly indicates SMC superiority over MC having in mind that computational effort for SMC is only about 10%-100% larger than for MC.
<figure><img src="image_2.png" /><figcaption>Figure 2: Relative efficiency of SMC estimator versus MC estimator measured by coefficient $\kappa$ versus number of time steps $N$ in the case of discretely monitored and continuously monitored barrier. If computing time for SMC estimator is $t_{SMC}$, then the computing time for MC estimator to achieve the same accuracy as SMC estimator is $\kappa \times t_{SMC}$.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 3: Relative standard error (in percent) of SMC and MC estimators in the case of continuously monitored barrier.</figcaption></figure>
<table><caption>Table 1: Comparison MC, $\hat{Q}_C^{MC}$, and SMC, $\hat{Q}_C^{SMC}$, option price estimators for continuously monitored barriers as the number of time steps $N$ increases. Exact price is 0.008061. Probability of the underlying asset not hitting the barrier during option life is $\psi = 0.005$.</caption><thead><tr><th>N</th><th>MC(stderr)</th><th>SMC(stderr)</th><th>κ</th></tr></thead><tbody><tr><td>1</td><td>0.008069(0.10%)</td><td>0.008074(0.12%)</td><td>0.39</td></tr><tr><td>2</td><td>0.008059(0.19%)</td><td>0.008077(0.13%)</td><td>1.14</td></tr><tr><td>4</td><td>0.008059(0.29%)</td><td>0.008064(0.14%)</td><td>1.58</td></tr><tr><td>8</td><td>0.008033(0.43%)</td><td>0.008046(0.15%)</td><td>3.98</td></tr><tr><td>16</td><td>0.008027(0.58%)</td><td>0.008066(0.12%)</td><td>8.15</td></tr><tr><td>32</td><td>0.008098(0.67%)</td><td>0.008063(0.13%)</td><td>13.25</td></tr><tr><td>64</td><td>0.008001(0.77%)</td><td>0.008070(0.13%)</td><td>16.12</td></tr><tr><td>128</td><td>0.007953(1.01%)</td><td>0.008050 (0.14%)</td><td>23.84</td></tr></tbody></table>

<table><caption>Table 2: Comparison MC, $\hat{Q}_D^{MC}$, and SMC, $\hat{Q}_D^{SMC}$, option price estimators for discretely monitored barriers as the number of time steps $N$ increases. $\psi$ is probability of the underlying asset not hitting the barrier.</caption><thead><tr><th>N</th><th>MC(stderr)</th><th>SMC(stderr)</th><th>κ</th><th>ψ</th></tr></thead><tbody><tr><td>1</td><td>0.8225(0.11%)</td><td>0.8229(0.12%)</td><td>0.69</td><td>0.359</td></tr><tr><td>2</td><td>0.5146(0.16%)</td><td>0.5140(0.10%)</td><td>2.11</td><td>0.229</td></tr><tr><td>4</td><td>0.2985(0.16%)</td><td>0.2985(0.10%)</td><td>2.19</td><td>0.137</td></tr><tr><td>8</td><td>0.1675(0.27%)</td><td>0.1684(0.11%)</td><td>4.98</td><td>0.080</td></tr><tr><td>16</td><td>0.0952(0.33%)</td><td>0.0957(0.11%)</td><td>7.25</td><td>0.048</td></tr><tr><td>32</td><td>0.0568(0.44%)</td><td>0.0566(0.13%)</td><td>10.54</td><td>0.029</td></tr><tr><td>64</td><td>0.0358(0.57%)</td><td>0.0361(0.13%)</td><td>17.84</td><td>0.019</td></tr><tr><td>128</td><td>0.0246(0.66%)</td><td>0.0249(0.14%)</td><td>20.12</td><td>0.013</td></tr></tbody></table>
# 7 Conclusion and Discussion

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

* Standard error of SMC estimator does not grow as the number of time steps increases while standard error of MC estimator can increase significantly. This is because in SMC, sampled asset values (particles) rejected by barrier condition are re-sampled from asset values between the barriers and thus the number of particles between the barriers will not change while in MC the number of simulated paths not breaching the barrier will reduce as the number of time steps increases.

* Efficiency of SMC versus standard MC improves when probability of asset path to hit the barrier increases (e.g. upper and lower barrier are getting closer or number of time steps increases). Typically, most significant benefit of SMC is achieved for cases when probability of not hitting the barrier is very small. Otherwise its efficiency is comparable to standard MC.

* Implementation of SMC requires little extra effort when compared to the standard MC method.

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

Further research may consider development of SMC and MC for alternative solution presented in Section 2.2. Also note that it is straightforward to calculate knock-in option as the difference between vanilla option (i.e. without barrier) and knock-out barrier option, however it is not obvious how to develop SMC estimator to calculate knock-in option directly (i.e. how to write knock-in option price expectation via Feynman-Kac representation formula (14)) which is a subject of future research. It is also worth to note that in this paper we focused on the case of one underlying asset for easy illustration while presented SMC algorithm 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 conflict of interests. The authors alone are responsible for the
writing of this work.

**References**

Andersen, L., & Brotherton-Racliffe, R. 2006. Exact Exotics. *Risk*, **9**(10), 85–89.

Bartoli, N., & Del Moral, P. 2001. *Simulation & Algorithmes Stochastiques*. Cépaduès éditions.

Beaglehole, D. R., Dybvig, P. H., & Zhou, G. 1997. Going to extremes: Correcting Simulation Bias in Exotic Option Valuation. *Financial Analyst Journal*, January/February, 62–68.

Borodin, A., & Salminen, P. 1996. *Handbook of Brownian Motion-Facts and Formulae*. Basel: Birkhauser Verlag.

Broadie, M., Glasserman, P., & Kou, S. 1997. A continuity correction for discrete barrier options. *Mathematical Finance*, **7**, 325–349.

Carmona, R., Fouque, J.-P., & Vestal, D. 2009. Interacting Particle Systems for the Computation of Rare Credit Portfolio Losses. *Finance and Stochastics*, **13**(4), 613–633.

Carmona, René, Del Moral, Pierre, Hu, Peng, & Oudjane, Nadia. 2012. An introduction to particle methods with financial applications. *Pages 3–49 of: Carmona, René, Del Moral, Pierre, Hu, Peng, & Oudjane, Nadia (eds), Numerical methods in finance*. Springer.

Cérou, F., Del Moral, P., & Guyader, A. 2011. A nonasymptotic theorem for unnormalized Feynman-Kac particle models. *Ann. Inst. Henri Poincaré Probab. Stat*, **47**(3), 629–649.

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

Del Moral, P. 2004. *Feynman-Kac Formulae. Genealogical and interacting particle approximations*. Probability and Applications. Springer.
Del Moral, P. 2013. *Mean field simulation for Monte Carlo integration*. Monographs on Statistics and Applied Probability. Chapman and Hall/CRC.

Del Moral, P., & Patras, F. 2011. Interacting path systems for credit risk. *Pages 649–674 of: Brigo, D., Bielecki, T., & Patras, F. (eds), *Credit Risk Frontiers*. WileyBloomberg Press.

Del Moral, Pierre, Hu, Peng, Oudjane, Nadia, & Rémillard, Bruno. 2011. On the Robustness of the Snell envelope. *SIAM Journal on Financial Mathematics, 2*(1), 587–626.

Del Moral, Pierre, Hu, Peng, & Oudjane, Nadia. 2012a. Snell envelope with small probability criteria. *Applied Mathematics & Optimization, 66*(3), 309–330.

Del Moral, Pierre, Rémillard, Bruno, & Rubenthaler, Sylvain. 2012b. Monte Carlo approximations of American options that preserve monotonicity and convexity. *Pages 115–143 of: Carmona, René, Del Moral, Pierre, Hu, Peng, & Oudjane, Nadia (eds), Numerical methods in finance*. Springer.

Dewynne, J., & Wilmott, P. 1994. Partial to exotic. *Risk Magazine*, December, 53–57.

Giles, M. 2008a. Improved multilevel Monte Carlo convergence using the Milstein scheme. *Pages 343–358 of: Monte Carlo and quasi-Monte Carlo methods 2006*. Springer.

Giles, M. 2008b. Multilevel Monte Carlo path simulation. *Operations Research, 56*(3), 607–617.

Glasserman, P. 2004. *Monte Carlo methods in financial engineering*. Springer.

Glasserman, P., & Staum, J. 2001. Conditioning on one-step survival for barrier option simulations. *Operations Research, 49*(6), 923–937.

Gobet, E., & Menozzi, S. 2010. Stopped diffusion processes: boundary corrections and overshoot. *Stochastic Processes and their Applications, 120*(2), 130–162.

Gobet, Emmanuel. 2009. Advanced Monte Carlo methods for barrier and related exotic options. *Handbook of Numerical Analysis, 15*, 497–528.

He, Hua, Keirstead, William P, & Rebholz, Joachim. 1998. Double lookbacks. *Mathematical Finance, 8*(3), 201–228.
Heynen, Ronald, & Kat, Harry. 1994a. Crossing barriers. *Risk*, **7**(6), 46–51.

Heynen, Ronald, & Kat, Harry. 1994b. Partial barrier options. *The Journal of Financial Engineering*, **3**(3), 253–274.

Hull, John C, & White, Alan D. 1993. Efficient procedures for valuing European and American path-dependent options. *The Journal of Derivatives*, **1**(1), 21–31.

Jasra, A., & Del Moral, P. 2011. Sequential Monte Carlo methods for option pricing. *Stochastic analysis and applications*, **29**(2), 292–316.

Johannes, M. S., Polson, N. G., & Stroud, J. R. 2009. Optimal Filtering of Jump Diffusions: Extracting Latent States from Asset Prices. *Review of Financial Studies*, **22**(7), 2759–2799.

Karatzas, I., & Shreve, S. 1991. *Brownian Motion and Stochastic Calculus*. Springer.

Kat, Harry M., & Verdonk, Leen T. 1995. Tree surgery. *Risk Magazine*, **8**(2), 53–56.

Kunitomo, Naoto, & Ikeda, Masayuki. 1992. Pricing Options With Curved Boundaries. **1. Mathematical finance**, **2**(4), 275–298.

Peters, G. W., Brier, M., Shevchenko, P., & 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.

Rubinstein, Mark, & Reiner, Eric. 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.

Shevchenko, P. V. 2011. Closed-form transition densities to price barrier options with one or two underlying assets. *CSIRO technical report EP11204*.

Targino, R. S., Peters, G. W., & Shevchenko, P. V. 2015. Sequential Monte Carlo Samplers for capital allocation under copula-dependent risk models. *Insurance: Mathematics and Economics*, **61**, 206–226.