Efficient Monte Carlo Algorithm Using Antithetic Variate and Brownian Bridge Techniques for Pricing the Barrier Options with Rebate Payments

Email: hmzubaidi@uqu.edu.sa Abstract: The down-and-out call barrier option with rebate payment on dividend-paying stock is simulated using a new version of the Monte Carlo algorithm. The standard Monte Carlo method for simulating such an option suffers from two sources of errors: Hitting time error inherent from time stepping and the Monte Carlo statistical error. We present a modified version of Monte Carlo method that can reduce these errors efficiently using the Brownian bridge technique for the hitting time error and the antithetic variate approach for the statistical error. We found that the Brownian bridge technique is responsible for improving the order of convergence in hitting time from one half to one and the antithetic variate technique can speed up the Monte Carlo simulation by reducing the variance of the computed payoff, giving almost twice as much accuracy. The standard error and the coefficient of variation are applied in order to measure the effectiveness of the volatility of the underlying option.


Introduction
A European option is a contract giving the option holder the right, but not the obligation to buy ( called call option) or to sell (called put option) an asset at a prescribed time, known as expiry time or maturity T, for a certain price K, known as strike price (Hull, 2009). Throughout this paper, we restrict ourselves on the European call option and the other type can be dealt with, in the same manner.
Basically, evaluating the price of European option, (V c for the call option), is based on analysis of the classical Black-Scholes model (Black and Sholes, 1973), where the money market consists of a riskless cash bond whose value at time t is described by B(t) = e rt , where r is the continuously compounded annual risk-free interest rate, and a single risky asset whose price S(t) follows a Geometric Brownian motion defined by: where, W(t) is a standard Brownian motion and µ and σ are constants denoting the expected rate of return and the volatility of the asset price, respectively (Hull, 2009;Black and Sholes, 1973). Using Ito's formula, a closed form solution for the asset price defined by (1) in its Ito form is obtained as (Kloeden and Platen, 2011): At expiry time, T, the payoff of the call option is calculated as (Black and Sholes, 1973): The price of European call option then can be found by discounting the expected of its payoff C(T) at the risk-free rate with respect to the risk-neutral probability measure Q(i.e., µ = r in (2)). Thus: This class of option is known as standard or vanilla call European option on non dividend-paying stock.
However, for dividend-paying stock as we consider in this study, the stock price S in risk-neutral world at expiry time T can be calculated by putting µ = rq in (2), where q is dividend payable (Hull, 2009) Other type of European options that become increasingly popular in world markets, are options with barrier features, commonly called European barrier options. Unlike the standard European options, the payoff of European barrier options depends not only on the price at expiry time but also on whether or not the underlying asset price has reached a specified boundary during the lifetime of options (Hull, 2009).
Two basic types of European barrier options can be identified: Knock-out and knock-in options. Knock-out option is one where the option is worthless if the asset price touches or passes the barrier. In contrast, the knock-in option is activated if the asset price hits or crosses the predefined barrier level. The attractiveness of such options is that they can be cheaper than their counterparts of standard European options, since they risk either not being knocked in or being knocked out. In this study, we restrict our attention to single barrier knock-out options, particularly down-and-out options. Using out-in parity, the value of the corresponding knock-in can be found as the difference between the vanilla call option and the down-and-out call option.
Valuation formulas have been derived for the barrier options in the literature. Merton (1973) provided the first analytical formula for a down-and-out call option which was the oldest barrier option type. Reiner and Rubinstein (1991) provide closed form solutions for the prices of all four types of barrier on both call and put options with the assumption that the underlying asset price follows Geometric Brownian process. Later in (Haug, 1998) gave a generalization of these formulas that were provided by Reiner and Rubinstein (1991).
Some barrier options specify that a fixed rebate R is to be given to the holder of the options if the knock-out and knock-in options become worthless. This can make the barrier options more attractive to the potential purchasers by compensating them for the loss of the option when the knock-out option is knocked out or when knock-in option is never knocked in (Chriss, 1997). Here, we will consider the down-and-out barrier option with rebate payment, and this case is considered as the first exit time problem.
The analytical expressions are not available for many cases of barrier options such as options with multiple assets and path-dependent options. Numerical techniques, particularly Monte Carlo algorithms, play crucial role in such a situation. Monte Carlo algorithms are a class of algorithms based on simulating a large number of price paths of the underlying assets under the risk neutral probability measure, computing the payoffs, averaging, and discounting at the risk free rate (Glasserman, 2003). Boyle (1997) provides an early discussion of the use of Monte Carlo methods in pricing options, and since that time, this technique has become one of the standard methods for the calculation of option values due to its flexibility in handling the increasing complexity of options problems.
However, in the case of simulating barrier options, the standard Monte Carlo algorithm becomes inefficient since it yields high statistical and discretization errors, due to the knockout feature of such options. The method has slow convergence and produces statistical errors with order of convergence , where M is the number of simulations (Glasserman, 2003;Higham, 2004) and discretization error or hitting time error with order of convergence of , with N time steps (Gobet, 2009).
Basically, the 95% confidence interval of the underlying option is defined as (Glasserman, 2003;Higham, 2004): where, V M is the Monte Carlo estimation of the option value and Var represents the computed variance of the value of the barrier option. Observe that the quantity Var M represents the standard error in Monte Carlo Simulation. Obviously, one way to reduce such an error and to improve the Monte Carlo approximation is by taking more samples, which makes it an expensive work in terms of computational time. For example, improving the accuracy by getting an extra digit, that is, shrinking a confidence interval by a factor of 10, requires 100 times as many samples (Higham, 2004). Moreover, if we would reduce the variance by a half, we need to add four times of samples. The standard error is also proportional to the square root of the variance and this motivates an alternative approach by concentrating on reducing the variance other than adding more samples. Such techniques are known as variance reduction techniques and its most widely used technique is the antithetic variate technique. For background on this technique and more details about it, see as a brief sampling (Hull, 2009;Glasserman, 2003;Boyle, 1997;Higham, 2004;Boyle et al., 1997). There are many other variance reduction approaches such as control variate, stratified sampling, importance sampling and moment matching methods that have been discussed in literature; see for example (Glasserman, 2003) and the references given there. Recently, Giles and Szpruch (2013) used Giles' multilevel Monte Carlo technique for pricing many financial problems including barrier options in order to minimize the overall variance for a given computational cost.
For continuously monitored barrier options, the MC method also suffers from the hitting time error inherent from time stepping. As mentioned above, this error is proportional to 1 N , due to the fact that the standard MC algorithm overestimates the actual values of the option prices because the option may be knocked out between the discrete computational nodes (Gobet, 2009). Inspired by (Mannella, 1999), one can improve such a slow convergence by using Brownian interpolation to calculate the exit probability of the barrier. This technique is known as a Brownian bridge technique and according to (Gobet, 2000), this simple technique is responsible for improving the order of convergence from under some conditions on functional of the stock price. For more analysis on using this technique in barrier options, we refer the reader to (Glasserman, 2003;Gobet, 2009;Baldi, 1995;Moon, 2008) and for other applications, we refer to (Alzubaidi and Shardlow, 2014) for neuroscience models and (Mannella, 1999;Jansons and Lythe, 2003;Buchmann, 2005) for physical systems. The rest of the paper is arranged as follows. In section 2, the problem of the down-and-out call barrier option with rebate payment is formulated and necessary notations are introduced. In section 3, we compute the analytical expressions for the value of the underlying option. In section 4, a Monte Carlo algorithm for simulating the barrier option is discussed. In addition, we introduce a variance reduction technique called antithetic variate technique and an error reduction approach known as a Brownian bridge technique. Moreover, we present a modified version of Monte Carlo algorithm by combining these techniques with the standard MC method. In section 5, we include numerical experiments concerning the underlying option, in order to compare the standard MC for efficiency and accuracy with the modified version of MC algorithm. The resultant errors are analyzed and the effectiveness of the volatility is studied. Section 6 contains our conclusions and some ideas for future works.

Problem Formulation
Let H, where H < S 0 , denotes the level of the barrier of the down-and-out call barrier option with rebate payment R. The first hitting time of the stock price S t with barrier H, where 0 < t < T and T is expiry time, is thus defined by (Glasserman, 2003;Moon, 2008): The discounted payoff of such an option is thus given by (Moon, 2008): The value of the option price at time t = 0 can be formally given as (Moon, 2008): where, the expectation here is taken with respect to the risk-neutral probability measure and K is strike price of the barrier option. The rebate here is called non-deferred rebate since the rebate R is paid as soon as the barrier is hit.

Explicit Formulas
Theorem 3.1. (The Barrier Option with Zero Rebate) , then the theoretical value of V doc at time t = 0 of a down-and-out call option, with strike price K, volatility σ, dividend payable q, risk-free interest rate r, barrier H and expiration time T, is given by: Where: N denotes the standard normal distribution function defined as: Horfelt (2003).

Theorem 3.2. (The Pure Rebate Option)
The theoretical value of Rebate option in down-andout call option V doc at time t = 0 is given by: Where: Proof. For the proof, Horfelt (2003). Consequently, the value of the down-and-out call barrier option with rebate payment can be calculated as the sum of the values of the options given in (8) and (9). Thus: Which will be used to check the efficiency of the numerical simulation algorithms that will be discussed in next section.

Monte Carlo Algorithm with Error and Variance Reduction Techniques
We begin with discretization of the time interval [0, We next compute S n+1 = S(t n+1 ) at each time step for n = 0, 1, ···, N-1 by: where, ∆W n = W(t n+1 )-W(t n ) are increments in the Brownian motion, which are normally distributed, and can be simulated as n n W tη ∆ = ∆ where each η n is independent and identically distributed standard normal random variable (Glasserman, 2003;Higham, 2004).
The standard Monte Carlo approximation for the expected value of the discounted payoff of the downand-out call barrier option with rebate payment defined in (6) can be calculated as (Glasserman, 2003;Higham, 2004;Moon, 2008): where, M is the simulation number, f j is the discounted payoff at simulation step j and: is the approximation for the first hitting time τ H . Then the resultant error produced using the standard MC can be written as (Moon, 2008): , where, ξ H and ξ s represent the hitting time error and the statistical error, respectively. In the following subsections, we will discuss the ways to reduce these errors using the Brownian bridge technique for the hitting time error and antithetic variate approach for the statistical error.

Brownian Bridge Technique
Using time stepping methods to simulate the payoff of the barrier options may cause large errors in the calculation of the probability of first time the option price reaches the barrier. This is due to the possibility that the process may attain the barrier and come back, within the time step. Approximating the continuous sample paths of Brownian motion using discrete random walks gives the values only at the beginning and the end of the time step and so we have no information about the behaviour of the continuous process during the time step (Buchmann, 2005).
Inspired by (Mannella, 1999;Moon, 2008) dealt with this situation by applying a simple hitting test after each time step using the distribution of the Brownian bridge from S(t n ) = S n to S(t n+1 ) = S n+1 . To be more precise, the price S t for t ∈ [t n , t n+1 ], can be considered as a constantdrift Brownian bridge pinned at S(t n ) = S n to S(t n+1 ) = S n+1 and therefore, the distribution of the first hitting time through the barrier H with respect to this bridge can be approximated as: To check if the barrier H is crossed during the time step [t n , t n+1 ], we sample a uniformly distributed random variable u n ∼ U(0, 1) and then compare it to the probability exit P n defined by (12). If u n < P n , we assume that the exit event has occurred and in this case, a predefined rebate R need to be paid. We consider S τH ≤ H for some τ H ∈ (t n , t n+1 ), which means that the discounted payoff can be given as e -rτH R. Otherwise, the process S t is assumed to remain within the time step for all t ∈ [t n , t n+1 ]. For an approximation to the first hitting time τ H , one can choose the midpoint of the time interval.

Antithetic Variate Technique
The antithetic variate technique attempts to reduce the variance by exploiting the existence of negative correlations between two estimates of the underlying option (Merton, 1973). To illustrate how to use this technique, suppose our objective is to estimate the option value defined by (7). From (11), the risk neutral stock price S n , n = 0, 1, ··· , N are generated using the sequence of independent and identically distributed standard normal random variables η n , n = 0, 1, ··· , N. The idea behind such a technique is that if η n are used to simulate the increment of Brownian path, then their antithetic variates -η n , which are also normally distributed random variables, can be used to simulate the increment of the reflection of the path about the origin (Giles and Szpruch, 2013). The output computed using the increment of Brownian path will be negatively correlated with that obtained by its reflection, resulting in a reduction in variance.
Recall that the standard MC estimator for the option value defined in (7) is given as: where, f j is the discounted payoff of the underlying barrier option at the simulation step j.
The antithetic alternative estimator is simply the averages of all of 2M observations (Glasserman, 2003;Higham, 2004): where, the j f is the discounted payoff of the option at the simulation step j which is obtained by changing the sign of all random standard normal samples used for calculating f j . The pairs f j and j f are called the antithetic variates. The V A is thus the sample mean of M independent observations: converges in distribution to N(0, 1) (Glasserman, 2003;Higham, 2004). We now have: < in order to reduce the computed variance, and that will be valid when f is monotonic (Glasserman, 2003;Higham, 2004).
Based on these techniques, the modified version of MC algorithm is carried out as follows: Algorithm 1: Modified version of Monte Carlo algorithm for the down-and-out call barrier option with rebate payment 1: for j = 1 to M do 2: for n = 1 to N do 3: Generate an N(0, 1) samples η n : n = 1, 2, ···, N,

Results and Discussion
In our numerical experiments, we will simulate the down-and-out call option with rebate payment using the algorithms considered in previous section. The efficiency of such algorithms will be examined by comparing the simulated results to the analytical results obtained by (10). The statistical and hitting time errors will be analyzed. In addition, we will study the effectiveness of volatility on the pricing of underlying option using such algorithms. Figure 2 displays the CPU time against the error in computed discounted expected payoff of the underlying option for the two methods. To avoid any influence from statistical errors, M is to be chosen as 10 6 . The number of the time steps is adjusted in each method until we reach the desired accuracy. As seen from the plot, the modified version of MC has a significantly better computational time required to achieve a given accuracy level than the standard MC. For example, in order to achieve the accuracy level 0.03, we needed 30000 time steps in each path for the standard MC and only 40 time steps in each path for the modified method. The ratio between the two CPU time in this case is equal to 159.6783, which means that the modified version of MC is approximately 160 times more efficient than the standard MC.  The analysis of these results as follows. The approximations of down-and-out call option obtained by standard MC algorithm overestimate the actual values, even for small values of ∆t, because the exact path of the underlying option may attain the barrier and comeback, within the time step, causing the hitting time error. Inspired by (Mannella, 1999), one can use conditional exit probability after each time step to correctly check whether the barrier is being hit or not. As discussed above, this approach is known as a Brownian bridge technique and according to (Gobet, 2000), this simple technique is responsible for improving the order of convergence from half order to one order, which coincides with our numerical observations shown in Fig.  1. The overall accuracy depends not only on the hitting time error, but also on the sampling error arising from Monte Carlo approximation to the mean of the option, which decreases as more realizations are generated. Thus the sampling error is dominated and in order to estimate the option value with reasonable accuracy, a large value of M is required, which is computationally expensive.

Convergence Properties
As shown in Table1, using the antithetic variate technique can lead to saving in computation time by reducing the variance of option value. To be more precise, Table 1 shows the 95% confidence intervals and the ratio of their widths (1.96×standard errors) for time step ∆t = 0.0025 (N = 400) and M = 10 4 , 10 5 , 10 6 , 10 7 using the underlying numerical techniques. The standard MC is overestimating the option values, even for largest sample size M = 10 7 . Thus the confidence interval for M = 10 7 , which is [5.5937, 5.6073], does not even include the theoretical value of the down and-out call option (5.2835), whereas for each sample size M, this analytical value is always contained in the corresponding confidence interval when the antithetic variate technique is used. Moreover, the antithetic variate technique shrinks the confidence intervals by a factor of around 1.6. Observe that the width of the corresponding confidence interval scales with the square root of the variance of the underlying option and is inversely proportional to the square root of the number of simulations M, which makes it an expensive work to improve the approximation by taking more samples. Broadly speaking, the use of Brownian bridge technique and the antithetic variate approach significantly speeds the convergence and improves the computational efficiency compared to the standard MC.

Volatility
The stock price volatility σ is critically important to option price since it is considered as a measure of the uncertainty about stock returns. Figure 3 and 4 show the effectiveness of the volatility σ on the pricing of the down-and-out call option with rebate using the standard MC and the new version of MC algorithm. We chose ∆t = 0.0025, M = 10 6 , with keeping other parameters as in Fig. 1. Following (Hull, 2009), the stocks typically have a volatility between 15 and 60%. In Fig. 3, we plot the computed mean value of underlying option against these values of volatility using the two numerical methods. Standard MC results are represented using squares symbols whereas the circles represent the results obtained using the MC with the variance and error reduction techniques. The estimation of 95% confidence interval of each computation is shown as error bars. For reference values of our computations, the theoretical Black-Scholes prices are plotted as shaded circles.
There is excellent agreement between the analytical Black-Scholes values and those obtained using the MC algorithm combined with antithetic variate approach and Brownian bridge technique for given values of volatility. Moreover, all Black-Scholes values are contained in the corresponding 95% confidence intervals of simulation results obtained by such a technique. As σ increases, the results obtained using the modified MC remain stable and converge to the corresponding analytical results. In contrast, as volatility σ increases, the standard MC simulated results go up and diverge away from the analytical results, giving poor convergence with wider widths of confidence intervals.
To understand the behaviour of these results, we plot in Fig. 4 the variance of underlying option against the volatility σ for each method, with M = 10 6 and N = 400. As we increase the volatility more, the corresponding computed variance increases. Moreover, It is clearly seen that using the modified MC method greatly reduces the variance of the estimates by approximately above 50% for each volatility, persistent with the theoretical framework discussed in section 4.     Table 2 summarizes the results of evaluating the standard errors, for the underlying option with different values of volatility σ, chosen as in Fig. 3. First, we observe that the standard errors are increasing functions of volatility σ. The use of antithetic variate technique affects the estimation of the standard errors. For example, the standard errors have been reduced from ±0.0084 to ±0.0053 for σ = 0.15, and from ± 0.031 to ±0.0194 for σ = 0.60, giving almost twice as much accuracy as the standard MC. Moreover, our results indicate that the standard errors obtained using the modified version of MC algorithm usually much less than the standard MC results calculated using 2M samples (Hull, 2009).
In the investing world, the coefficient of variation CV or the ratio of standard deviation to the expected payoff can be considered as a measure of how much volatility we are assuming in comparison to the expected return from an investment. There for, an investment with lower value for the CV is preferred to a one with higher CV. Table 3 displays the coefficient of variation CV as a function of volatility σ that are chosen as in Fig. 3. The results indicate that the values of CV obtained by standard MC have been reduced by around a third for the low-level of volatility and around a quarter for high-level of volatility, when the modified MC algorithm is used.
To be more precise, the ratio of the CV new to the CV stand raises from 0.6575 to 0.7658 as we increase the volatility σ from 15% to 60%. Moreover, for the chosen parameters, the CV is always greater than one, which is the indication of a greater variability in the price of the underlying option.

Conclusion
A barrier option is a type of path-dependent option whose payoff at maturity depends on whether or not the underlying asset price has hit a specified barrier during the life of the option. In our work, we considered the down-and-out call barrier option with rebate payment. Explicit expressions for underlying option were derived. We applied the standard Monte Carlo algorithm for simulating such an option. The statistical and hitting time errors obtained by standard MC algorithm have been analyzed. We discussed the ways to reduce these errors using the antithetic variate technique for the statistical error and the Brownian bridge technique for the hitting time error. We found that the Brownian bridge technique is responsible for improving the order of convergence in hitting time from one half to one, coinciding with the theoretical work done previously; (Gobet, 2009;2000). We also found that the antithetic variate technique can speed up the Monte Carlo simulation by reducing the variance of the computed option price by around 50%, with the same number of realizations. In statistical problems, the parameters are usually estimated using the 95% confidence interval to guarantee capturing the explicit value of the underlying option, with high probability. Indeed, using the antithetic variate technique shrinks the 95% confidence interval of the simulated value of the underlying option efficiently, giving almost twice as much accuracy; see for more analysis (Glasserman, 2003;Higham, 2004). A modified version of Monte Carlo algorithm has been introduced using these reduction techniques. Finally, the standard error and the coefficient of variation have been applied to measure the effectiveness of the volatility computed using the standard and the modified version methods. Based on the Black-Scholes framework, the volatility is considered to be constant during the life of the option; however, this is usually inconsistent with real stock markets. By assuming that the volatility of the stock price is subject to fluctuations, it becomes possible to model the options more accurately. The models with this feature are known as stochastic volatility models either with jumps such as Bates model or without jumps such as Heston model. This would be an interesting area to consider in our future work.

Ethics
This article is original and contains unpublished material. The corresponding author confirms that all of the other authors have read and approved the manuscript and no ethical issues involved.