Improving Monte Carlo: Control Variates

I’ve already discussed quite a lot about Monte Carlo in quantitative finance. MC can be used to value products for which an analytical price is not available in a given model, which includes most exotic derivatives in most models. However, two big problems are the time that it takes to run and the ‘Monte Carlo error’ in results.

One technique for improving MC is to use a ‘control variate’. The idea is to find a product whose price is strongly correlated to the product that we’re trying to price, but which is more easy to calculate or which we already know the price of. When we simulate a path in MC, it will almost surely give either an under-estimate or an over-estimate of the true price, but we don’t know which, and averaging all of these errors is what leads to the Monte Carlo error in the final result. The insight in the control variate technique is to use the knowledge given to us by the control variate to reduce this error. If the two prices are strongly correlated and a path produces an over-estimate of product price, it most likely also produces an over-estimate of the control variate and visa versa, which will allow us to improve our estimate of the product we’re trying to price.

The textbook example is the Asian Option. Although the arithmetic version of the asian option discussed in previous posts has no analytic expression in BS, a similar Geometric asian option does have an analytic price. So, for a given set of model parameters, we can calculate the price of the option. As a reminder, the payoff of an arithmetic asian option at expiry is

    \[C_{\rm arit}(T) = \Bigl({1\over N}\sum_{i=0}^{N-1} S(t_i) - K \Bigr)^+\]

and the payoff of the related geometric averaging asian is

    \[C_{\rm geo}(T) = \Bigl( \bigl(\prod_{i=0}^{N-1} S(t_i)\bigr)^{1\over N} - K \Bigr)^+\]

Denoting the price of the arithmetic option as X and the geometric option as Y, the traditional monte carlo approach is to generate N paths, and for each one calculate X_i (the realisation of the payoff along the path) and take the average over all paths, so that

    \[{\mathbb E}[X] = {1 \over N} \sum_{i=0}^{N-1} X_i\]

which will get closer to the true price as N \to \infty.

Using Y as a control variate, we instead calculate

    \[{\mathbb E}[X] = {1\over N} \sum_{i=0}^{N-1}\bigl( X_i - \lambda( Y_i - {\mathbb E}[Y] ) \bigr)\]

where {\mathbb E}[Y] is the price of the geometric option known from the analytical expression, and \lambda is a constant (in this case we will set it to 1).

What do we gain from this? Well, consider the variance of X_i - \lambda (Y_i - {\mathbb E } [Y])

    \[{\rm Var}\bigl( X_i - ( Y_i - {\mathbb E}[Y] ) \bigr) = {\rm Var}( X_i ) + \lambda^2 {\rm Var} ( Y_i ) - 2 \lambda {\rm Cov}( X_i Y_i )\]

(since {\mathbb E } [Y] is known, so has zero variance) which is minimal for \lambda = \sqrt{ {\rm Var}(X_i) \over {\rm Var}(Y_i) }\cdot \rho(X_i,Y_i) in which case

    \[{ {\rm Var}\bigl( X_i - ( Y_i - {\rm E}[Y] ) \bigr) \over {\rm Var}( X_i )} = 1 - \rho(X_i,Y_i)^2\]

that is, if the two prices are strongly correlated, the variance of the price calculated using the control variate will be a significantly smaller. I’ve plotted a sketch of the prices of the two types of average for 100 paths – the correlation is about 99.98%. Consequently, we expect to see a reduction in variance of about 2000 times for a given number of paths (although we now have to do a little more work on each path, as we need to calculate the geometric average as well as the arithmetic average of spots). This is roughly 45 times smaller standard error on pricing – well over an extra decimal place, which isn’t bad – and this is certainly much easier than running 2000 times as many paths to achieve the same result.

The relationship between payout for a geometric and an arithmetic asian option, which here demonstrate a 99.98% sample correlation
The relationship between payout for a geometric and an arithmetic asian option, which here demonstrate a 99.98% sample correlation. Parameters used were: r: 3%; vol: 10%; K: 105; S(0): 100; averaging dates: monthly intervals for a year

Asian Options II: Monte Carlo

In a recent post, I introduced Asian options, an exotic option which pay the average of the underlying over a sequence of fixing dates.

The payoff of an Asian call is

C(T) = \Big({1 \over N}\sum^N_{i=0} S_{t_i} - K \Big)^+

and to price the option at an earlier time via risk-neutral valuation requires us to solve the integral

C(t) = P_{t,T}\ \big( \iint\cdots\int \big) \Big({1 \over N}\sum^N_{i=0} S_{t_i} - K \Big)^+ \phi(S_{t_1},S_{t_1},\cdots,S_{t_N}) dS_{t_1}dS_{t_2}\cdots dS_{t_N}

A simple way to do this is using Monte Carlo – we simulate a spot path drawn from the distribution \inline \phi(S_{t_1},S_{t_1},\cdots,S_{t_N}) and evaluate the payoff at expiry of that path, then average over a large number of paths to get an estimate of C(t).

It’s particularly easy to do this in Black-Scholes, all we have to do is simulate N gaussian variables, and evolve the underlying according to a geometric brownian motion.

I’ve updated the C++ and Excel Monte Carlo pricer on the blog to be able to price asian puts and calls by Monte Carlo – have a go and see how they behave relative to vanilla options. One subtlety is that we can no longer input a single expiry date, we now need an array of dates as our function input. If you have a look in the code, you’ll see that I’ve adjusted the function optionPricer to take a datatype MyArray, which is the default array used by XLW (it won’t be happy if you tell it to take std::vector). This array can be entered as a row or column of excel cells, and should be the dates, expressed as years from the present date, to average over.

Excel Monte-Carlo Pricer

I’ve updated the Monte-Carlo pricer to run in excel. Give it a go and the benefits should be obvious – no more messing around with consoles or repeated typing, excel provides a handy interface for running C++ files.

What you’ll need to do is:

  • Download the Excel pricing sheet
  • Download the .xll library (large!) – this is the equivalent of the old .exe file
  • Open the Excel file, and drag-and-drop the .xll file inside the pricing sheet
  • Click on the option to “enable this add-on for this session only”
  • The Option Price function should now re-calculate option price as you make changes to the parameters!

Try using the option pricer function yourself – it now works just like a native excel function, so you can drag it into many cells and it will adjust target cells as usual allowing for easy variation of parameters and plotting of results.

If you’ve successfully set up XLW and Code::Blocks as I described in a previous post, I’ve also made the source code available in a .zip file which you can load in Code::Blocks and fiddle with yourself. I’ll be doing a bit of updating soon to make the code slightly cleaner and allow pricing of more interesting options.

The old console project source code is still available in compiled form and with the source code in a zip file at the bottom of the MC pricers page.

Digital Options

Today I’m going to talk about the valuation of another type of option, the digital (or binary) option. This can be seen as a bit of a case study, as I’ll present the option payoff and the analytical price and greeks under BS assumptions, and give add-ons to allow pricing with the MONTE CARLO pricer. I’ve also updated the ANALYTICAL pricer to calculate the price and greeks of these options.

Digital options are very straight-forward, they are written on an underlying S and expire at a particular date t, at which point digital calls pay $1 if S is greater than a certain strike K or $0 if it is below that, and digital puts pay the reverse – ie. the payoff is

P_{\rm dig\ call} = \Big\{ \ \begin{matrix} \$1 \quad S \geq K\\ \$0 \quad S < K\end{matrix}

We can calculate the price exactly in the BS approximation using the same method that I used to calculate vanilla option prices by risk-neutral valuation as follows

C_{\rm dig\ call}(0) = \delta(0,t)\ {\mathbb E}[ C_{\rm dig\ call}(t) ]

where \inline C_{\rm dig\ call}(t') is the price of the option at time \inline t', and we know that this must converge to the payoff as \inline t \to t', so \inline C_{\rm dig\ call}(t) = P_{\rm dig\ call}

= \delta(0,t)\ \int^{\infty}_{-\infty} P_{\rm dig\ call} \cdot e^{-{1\over 2}x^2} dx

\inline P_{\rm dig\ call} is zero if \inline S = S_0 e^{(r - {1\over 2}\sigma^2)t + \sigma \sqrt{t} x} < K which corresponds to \inline x < {\ln{K \over S_0} - (r-{1\over 2}\sigma^2)t \over \sigma \sqrt{t}} = -d_2

= \delta(0,t)\ \cdot \$1 \cdot \int^{\infty}_{-d_2} e^{-{1\over 2}x^2} dx

= \$ \delta(0,t)\cdot \Phi(d_2)


d_1 = {\ln{S\over K}+(r+{1\over 2}\sigma^2)t \over \sigma \sqrt{t}} \quad ;\quad d_2 = {\ln{S\over K}+(r-{1\over 2}\sigma^2)t \over \sigma \sqrt{t}}

Since we have an analytical price, we can also calculate an expression for the GREEKS of this option by differentiating by the various parameters that appear in the price. Analytical expressions for a digital call’s greeks are:

\Delta = {\partial C \over \partial S} = e^{-rt}\cdot {\phi(d_2)\over S\sigma \sqrt{t}}

\nu = {\partial C \over \partial \sigma} = -e^{-rt}\cdot {d_1 \ \phi(d_2)\over \sigma}

\gamma = {\partial^2 C\over \partial S^2} = -e^{-rt}\cdot {d_1 \ \phi(d_2)\over S^2 \sigma^2 t} = -e^{-rt}\cdot {d_1\ \phi(d_1)\over S K \sigma^2 t}

{\rm Vanna} = {\partial^2 C \over \partial S \partial \sigma} =e^{-rt}\cdot {\phi(d_2)\over S \sigma^2 \sqrt{t}}\Big( d_1 d_2 - 1 \Big)

{\rm Volga} = {\partial^2 C \over \partial \sigma^2} = e^{-rt}\cdot {\phi(d_2)\over \sigma^2}\cdot \Big( d_1 + d_2 - d_1^2 d_2 \Big)


\phi(d_1) = {1\over \sqrt{2 \pi}}\ {e^{-{1\over 2}d_1^2} }

Holding a binary put and a binary call with the same strike is just the same as holding a zero-coupon bond, since we are guaranteed to receive $1 wherever the spot ends up, so the price of a binary put must be

e^{-rt} = e^{-rt}\cdot \Phi(d_2) + C_{\rm dig\ put}(t'=0)\quad \quad \quad \quad [1]

C_{\rm dig\ put}(0) = e^{-rt} \cdot \Phi(-d_2)

Moreover, differentiating equation [1] above shows that the greeks of a digital put are simply the negative of the greeks of a digital call with the same strike.

Graphs of these are shown for a typical binary option in the following graphs. Unlike vanilla options, these option prices aren’t monotonic in volatility: if they’re in-the-money, increasing vol will actually DECREASE the price since it makes them more likely to end out-of-the-money!

Price and first-order greeks for a digital call option.
Price and first-order greeks for a digital call option.
Second-order greeks for a digital call option. Greeks for digital puts are simply the negative of these values

One final point on pricing, note that the payoff of a digital call is the negative of the derivative of a vanilla call payoff wrt. strike; and the payoff of a digital put is the positive of the derivative of a vanilla put payoff wrt. strike. This means that any binary greek can be calculated from the corresponding vanilla greek as follows

\Omega_{\rm dig\ call} = -{\partial \Omega_{\rm call}\over \partial K}

\Omega_{\rm dig\ put} = {\partial \Omega_{\rm put}\over \partial K}

where here \inline \Omega represents a general greek.

If you haven’t yet installed the MONTE CARLO pricer, you can find some instructions for doing so in a previous post. The following links give the header and source files for binary calls and puts which can be dropped in to the project in your C++ development environment

These will register the option types with the option factory and allow monte carlo pricing of the options (so far, all of the options in the factory also have analytical expressions, but I’ll soon present some options that can only be priced by Monte Carlo).

New Section: Monte Carlo in C++

Good news, this blog has a new section – Monte Carlo pricers! I’m going to experiment to find the most convenient way of making it available for download, and would like it to be available both as a compiled .exe and also as uncompiled source code to allow readers to examine and alter it.

The first iteration can be found online now in the MONTE CARLO section. It is code for a straight-forward Monte Carlo pricer that will price calls and puts, with uniform random variates selected by park-miller and converted into gaussian variates using one of the box-muller processes (for more information on these, refer to the post on Random Numbers).

I’ve used a factory pattern so in principle it should be very straight forward to add new options and also RNG processes to the code, although there are a few more improvements that I will be making in the near future, which will extend the code to basic path-dependent options and allow easier data input by the user. I’ll be going over various aspects of the code in detail in future posts, and also putting together some instructions on how to compile the code yourself using Dev C++ (a free, open-source compiler) and how to add options and processes via the factory.

As always, if you have any problems getting the code to work, please tell me!