Quasi-Monte Carlo Methods in Option Pricing

Quasi-Monte Carlo methods are finding a wider use in recent applications from computational finance. Their use is applied to a problem from option pricing.

by Giray Ökten

Monte Carlo and quasi-Monte Carlo methods describe a family of numerical tools that are used in numerous disciplines including economic modeling and finance. In general, quasi-Monte Carlo methods provide more accurate estimates than Monte Carlo methods. Their use  involves simulating the model under consideration using the so called "low-discrepancy" sequences, instead of pseudorandom numbers. We will examine some of these low-discrepancy sequences and write Mathematica codes to generate them. We will then illustrate their use in a problem from option pricing.

Monte Carlo and Quasi-Monte Carlo Methods

The general problem that Monte Carlo and quasi-Monte Carlo methods are applied to is to estimate the expected value E[X]=(∫_Ω)X dμ of a random variable X, defined on a probability space (Ω,μ). To elaborate on the theory of Monte Carlo and quasi-Monte Carlo methods, let us consider the case when for a particular choice for (Ω,μ) and X, the above problem reduces to estimate the definite integral of a function defined on the s-dimensional unit cube, i.e., to estimate I=(∫_((([0, 1])^s)))f(x)dx. 

In Monte Carlo methods, I is estimated by selecting N independent vectors (Subscript[x, 1]),...,(Subscript[x, N]) at random from the uniform distribution on (([0, 1])^s) and then evaluating the average \(TraditionalForm\`\(\(I\&^\)\(=\)\)\)(1/N) (∑_(n = 1)^N) f((Subscript[x, n])). In a practical implementation, the vectors (Subscript[x, 1]),...,(Subscript[x, N]) are usually generated from a deterministic algorithm. It is hoped that the numbers produced in such a deterministic manner imitate the true behavior of "random" numbers. We call these "imitations" pseudorandom numbers and their success is determined by their overall performance on a variety of statistical tests. From the strong law of large numbers, \(TraditionalForm\`I\&^ → I\) as N→∞, with probability one. In addition, from the central limit theorem, the error of the approximation, \(TraditionalForm\`\(\(I\&^ - I\)\(,\)\)\) is a random variable which is approximately a normal random variable with mean 0 and standard deviation σ/(Sqrt[N]), where (σ^2)=(∫_((([0, 1])^s)))(((f(x) - I))^2)dx. Therefore, Monte Carlo estimates satisfy the probabilistic error bound O(1/(Sqrt[N])). 

Quasi-Monte Carlo methods are often described as the deterministic versions of Monte Carlo methods. They emphasize the "uniformity" of the numbers we choose from (([0, 1])^s), rather than their randomness. Let {(Subscript[q, n])} be a sequence in (([0, 1])^s). The uniformity, or the evenness, of {(Subscript[q, n])} is usually measured by its star discrepancy \(TraditionalForm\`\(D\_N\%*\)(q\_1,  ... , q\_N)\) (or, \(TraditionalForm\`\(D\_N\%*\)(q\_n)\), in short) which is defined by: \(TraditionalForm\`\(D\_N\%*\)(q\_1,  ... , q\_N) = 
        sup\+\(αϵ[0, 1]\^s\) | \(\(1\/N\) \(∑\+\(n = 1\)\%N\( 1\_\(\([0, α\)\()\)\)\) \((q\_n)\)\) - \(\(α\_1\)\(·\)\) ... \)·α\_s | \).
Here, \(TraditionalForm\`\(\([0, α\)\()\)\) = \(\(\(\([0, α\_1\)\()\)\)\(×\)\) ... \)×\(\([0, α\_s\)\()\)\)\) and (Subscript[1, A]) (x) is the characteristic function of the set A which takes the value 1 if xϵA and 0 otherwise. A sequence {(Subscript[q, n])} is called a uniformly distributed sequence if \(TraditionalForm\`\(D\_N\%*\)(q\_1,  ... , q\_N)\)→0 as N→∞. If we construct the estimate \(TraditionalForm\`\(\(I\&^\)\(=\)\)\)(1/N) (∑_(n = 1)^N) f((Subscript[q, n])) using a uniformly distributed sequence {(Subscript[q, n])}, then \(TraditionalForm\`I\&^ → I\) as N→∞, from a well known result of the theory of uniform distribution (see [1]). An upper bound for the error of this approximation is given by the Koksma-Hlawka inequality, \(TraditionalForm\`\(\(|\)\(I - \(1\/N\) \(∑\+\(n = 1\)\%N f(
                q\_n)\)\)\(|\)\(\(≤\)\(\(V(f)\) \(\(D\_N\%*\)(
              q\_n)\)\)\)\)\), where V(f) is the variation of the function in the sense of Hardy and Krause ([1]), which is assumed to be finite. Therefore, we expect to obtain smaller error when using sequences with smaller star discrepancy.

The star discrepancy of the best known uniformly distributed sequences satisfies O((((log N))^s)/N). These sequences are called "low-discrepancy" sequences. Then, from the Koksma-Hlawka inequality, estimates obtained using a low-discrepancy sequence satisfy the deterministic error bound O((((log N))^s)/N). The classical low-discrepancy sequences are the van der Corput, Halton, Sobol, Faure and Niederreiter sequences, in the order of decreasing bounds for their discrepancy.

In the following section, we will study some of the low-discrepancy sequences in detail and write Mathematica codes to generate them. Then we will show how these sequences, as well as the pseudorandom sequence, can be used to estimate the value of an Asian option.

Low-Discrepancy Sequences

Van der Corput and Halton Sequences

One of the classical low-discrepancy sequences is the van der Corput sequence in base b, where b is any integer greater than one. To find the nth term, (Subscript[x, n]), of the sequence, we first write the unique digit expansion of n in base b: n = (∑_(j = 0)^∞) (Subscript[a, j])(n) 
            (b^j). This is the unique expansion that has only finitely many non-zero coefficients (Subscript[a, j])(n). The next step is to evaluate (Subscript[ϕ, b])(n), the radical inverse function in base b, which is defined as : (Subscript[ϕ, b])(
          n) = (∑_(j = 0)^∞) (Subscript[a, j])(n) 
            (b^(-j - 1)). The van der Corput sequence in base b is then given by (Subscript[x, n]) = (Subscript[ϕ, b])(n), for all n ≥ 0.

Let us illustrate in detail how we can use Mathematica to compute (Subscript[ϕ, 2])(6), a particular term of the van der Corput sequence in base 2. First we need to write 6 in its binary expansion. This can be done using the built-in Mathematica function, IntegerDigits. Here is the digit expansion of 6 in base 2: IntegerDigits[6,2] = {1, 1, 0}. The next step is to calculate the sum (Subscript[a, 2])(n) (2^(-3)) + (Subscript[a, 1])(n) 
          (2^(-2)) + (Subscript[a, 0])(n) (2^(-1)), where  (Subscript[a, 0]) = 0,  (Subscript[a, 1]) = 1, (Subscript[a, 2]) = 1.  We apply the Mathematica function Dot and use the listability of the exponent "^".

Our "list1" is Reverse[IntegerDigits[6,2]] = {0, 1, 1} and "list2" is Range[3] = {1,2,3}. We note that the argument of the function Range, the number 3, is the length of list1. We want to multiply each entry of list1 by 2 raised to the opposite of the corresponding entry of list2, and then add the products. That is, we want to compute:  0× (2^(-1)) + 
        1× (2^(-2)) + 1×(2^(-3)). This can be done by the pure function (#1*2^-#2)&. 
Generalizing the above expressions to any number n and base b, we obtain:

\(\(vanderCorput[n_, \ b_]\  := \ 
          IntegerDigits[n, 
              b] . \((b^
                Range[\(-Floor[
                      Log[b, n] + 1]\), \(-1\)])\);\)


      \(Attributes[vanderCorput]\  = \ Listable;\)\)

These are the first 10 van der Corput numbers in base 2, starting from n=1.

Table[vanderCorput[n, 2], {n, 10}]

Van der Corput numbers in pairwise relatively prime bases can be used to construct low-discrepancy sequences in higher dimensions. The s-dimensional Halton sequence in bases (Subscript[b, 1]),...,(Subscript[b, s]) is defined by, (Subscript[x, n]) = ((Subscript[ϕ, (Subscript[b, 1])])(n), (Subscript[ϕ, (Subscript[b, 2])])(n),  ... , (Subscript[ϕ, (Subscript[b, s])])(n)), where (Subscript[b, 1]),...,(Subscript[b, s]) are pairwise relatively prime numbers. Our choice for (Subscript[b, i]) is the ith prime number, which can be generated in Mathematica using the Prime function. The following code generates s-dimensional Halton vectors in bases (Subscript[p, 1]),...,(Subscript[p, s]), where (Subscript[p, i]) is the ith prime number.

Here are the first six three-dimensional Halton vectors:

Table[halton[n, 3], {n, 6}]

We now plot the first 1000 two-dimensional Halton vectors:

ListPlot[Table[halton[n, 2], {n, 1000}], Frame -> True];

Compare this graph with the plot of 1000 pseudorandom numbers:

ListPlot[Table[Random[], {1000}], Frame -> True];

It can be seen from the previous graphs that the Halton vectors are distributed more evenly than the pseudorandom vectors. It is indeed this evenness, or uniformity, that makes the low-discrepancy sequences work better than the pseudorandom numbers in general.

Scrambled Van der Corput and Halton Sequences

The uniformity of the van der Corput numbers can be further improved by permuting the coefficients in the digit expansion of n in base b. The scrambled van der Corput sequence in base b is defined by  
      (Subscript[s, b])(n) = 
        (∑_(j = 0)^∞) σ((Subscript[a, j])(n)) 
            (b^(-j - 1)), where σ is a permutation on the digits {0,1,...,b-1}. The following permutations, defined recursively, are given by Faure [2]:

Using these permutations, we now construct the scrambled van der Corput sequence in base b.

The scrambled Halton sequence is constructed using the scrambled van der Corput numbers.

These are the first six three-dimensional scrambled Halton vectors:

Table[scrhalton[n, 3], {n, 6}]

Faure Sequence

Another low-discrepancy sequence which satisfies smaller discrepancy bounds than the Halton sequence is the Faure sequence. Let x be a b-adic rational in [0,1] with expansion x = (∑_(j = 0)^∞) (Subscript[x, j]) 
            (b^(-j - 1)). Define a function P on b-adic rationals by P(
          x) = (∑_(j = 0)^∞) ξ((Subscript[x, j])) 
            (b^(-j - 1)). Here ξ((Subscript[x, j]))=(∑_(i≥j)) (({{i}, {j}}))(Subscript[x, i])  (mod b), where  (({{i}, {j}})) is the binomial coefficient. Then, the s-dimensional Faure sequence is defined by (Subscript[f, n]) = ((Subscript[ϕ, p])(n),P((Subscript[ϕ, p])(n)),  ... , (P^(s - 1))((Subscript[ϕ, p])(n))), where p is the smallest prime number greater than or equal to s.

To find the smallest prime number greater than or equal to a given positive integer s, we will use the function NextPrime, from the "NumberTheoryFunctions" package. This loads the package.

<< NumberTheory`NumberTheoryFunctions`

NextPrime[s] gives the smallest prime number larger than s. Then, the smallest prime number greater than or equal to s will be given by NextPrime[s - 0.1].

Here is the code that generates the s-dimensional Faure sequence. The function pascal produces the list ξ((Subscript[x, 1])),  ... , ξ((Subscript[x, N])) as the output when given  the list  (Subscript[x, 1]),...,(Subscript[x, N]) as the input. The powers of the function P, which is described above, are implemented using the built-in Nest function.

Here are some three dimensional Faure vectors:

Table[faure[n, 3], {n, 1, 10}]

And here is the plot of the first 1000 two-dimensional Faure vectors:

ListPlot[Table[faure[n, 2], {n, 1000}], Frame -> True];

Applications to Option Pricing: Evaluation of Asian Options

In this section we consider the pricing of average rate options. These options are Asian options whose payoff depends on the arithmetic average price of the underlying asset. The exercising features are assumed to be European. The asset price follows the Ito process dS = μSdt + σSdX, where t is time, μ is the drift, σ is the volatility of the asset and X=X(t) is a standard Brownian motion. The Ito process is discretized as: S(t+Δt) = S(t) + μ S(t)Δt + σS(t)(Sqrt[Δ])tZ , where Z is a random variable with the standard normal distribution.

For a finite sequence of asset prices (Subscript[S, 0]), (Subscript[S, 1]),  ... , (Subscript[S, K]), the payoff function for an average rate option is h((Subscript[S, 0]), (Subscript[S, 1]),  ... , (Subscript[S, K])) = 
        max((1/(K + 1)) (∑_(i = 0)^K) (Subscript[S, i]) - E, 0), where E is the exercise price. The value of the option at t = 0 is the expected value (1) E[(e^(-rT)) h((Subscript[S, 0]), (Subscript[S, 1]),  ... , (Subscript[S, K]))], conditioned on S(0) = (Subscript[S, 0]) and μ = r. Here T is the expiry and r is the risk free interest rate. A detailed analysis of these options can be found in Kemna & Vorst [3].

Valuation of Asian options by simulation is done by generating N random walks (Subscript[w, 1]),  ... , (Subscript[w, N]), where a random walk (Subscript[w, j]) determines K asset prices (Subscript[w, j]) = (Subscript[S, 1]^j),  ... , (Subscript[S, K]^j), at discretized time steps (Subscript[t, 1]),  ... , (Subscript[t, K]) = T (at (Subscript[t, 0]) = 0, the stock price is (Subscript[S, 0]), which is fixed for all random walks). The asset prices are generated using (2)  (Subscript[S, i]^j)=(Subscript[S, i - 1]^j) + 
              r (Subscript[S, i - 1]^j) Δt + σ (Subscript[S, i - 
                      1]^j) (Sqrt[Δ])t (Subscript[z, i]) , i = 1,2,...,K, where Δt = (T/K) is the time interval of the discretization. The number (Subscript[z, i]) is the ith number sampled from a standard normal distribution. The expected value (1) is then estimated by generating N random walks (Subscript[w, j]) = (Subscript[S, 1]^j),  ... , (Subscript[S, K]^j),  j = 1,...,N and evaluating the average (1/N) (∑_(j = 1)^N) (e^(-rT)) h(
                (Subscript[w, j])).

In the following, numberofrandomwalks, sini, exercise, sigma and deltaT stand for N, (Subscript[S, 0]), E, σ and Δt, respectively. We put:

\(K = 30; 
numberofrandomwalks = 5000; 
r = 0.1; 
sini = 40; 
      
exercise = 35; 
sigma = 0.4; 
T = 0.2; 
deltaT = T/K; \)

We will use the K-dimensional pseudorandom, Halton, scrambled Halton and Faure vectors to sample from the standard normal distribution in generating the random walks (Subscript[w, j]) from (2). To easily switch from one sequence to another, let us define:

\(\(	\(vector[i_] := 
      Which[seq == 1, Table[Random[], {K}], seq == 2, halton[i, K], seq == 3, 
        scrhalton[i, K], seq == 4, faure[i, K]]; \)\)\)

To sample from the standard normal distribution, we will use the following approximation to the inverse normal distribution [4].

\(\(zinv = 
      Compile[{x}, 
        Module[{c0, c1, c2, d1, d2, d3, p, t, output}, 
		c0 = 2.515517; 

          		c1 =  .802853; 
		c2 =  .010328; 
		d1 = 1.432788; 
		
          d2 =  .189269; 
		d3 =  .001308; 
		
          p = If[x >  .5, 1 - x, x]; 
		t = Sqrt[Log[1/p^2]]; 
		
          output = 
            t - \((c0 + c1*t + c2*t^2)\)/\((1 + d1*t + d2*t^2 + d3*t^3)\); 

          		If[x < 0.5, \(-output\), output]\ \ \ \ \ ]]; \)\)

We have used a compiled function to increase the speed of the calculations. Without compiling, the function zinv generates 10,000 numbers in 34.11 seconds. After compiling, it takes only 0.55 seconds to do the same operation! 

Here is the function that gives the next asset price (Subscript[S, i]^j) from  (Subscript[S, i - 1]^j) using (2). The second argument "innumber" will be selected from numbers sampled from the normal distribution.

nextprice[s_, innumber_] := 
      s + r*s*deltaT + sigma*s*Sqrt[deltaT]*zinv[innumber];

The expression, FoldList[nextprice,sini,vector[j]], gives the jth random walk. We then need to compute the average price and evaluate h((Subscript[S, 0]), (Subscript[S, 1]),  ... , (Subscript[S, K])) = h((Subscript[w, j])) = 
          max((1/(K + 1)) (∑_(i = 0)^K) (Subscript[S, i]) - E, 0). This is done by:

\(h[j_Integer] := 
        Max[Apply[Plus, FoldList[nextprice, sini, vector[j]]]
			
              \((numberofprices + 1)\) - exercise, 0]; \)

The final step is to calculate the average of these numbes as j runs from 1 to N. Here is the expression that evaluates (1/N) (∑_(j = 1)^N) 
                (e^(-rT)) h((Subscript[w, j])), an estimate for the value of the option using N random walks: 
Sum[h[j],{j,numberofsamples}]/numberofrandomwalks

It is often desirable to plot the estimates (1/n) (∑_(j = 1)^n) (e^(-rT)) h(
              (Subscript[w, j])) against the number of random walks n, as n = 1,...,N.  Such a graph will show how the estimates converge to the exact solution as n increases. It is also easier to compare the estimates obtained by the use of different sequences by comparing the graphs they produce. Using the FoldList function, we can generate the data efficiently:

Rest[FoldList[Plus, 0, Map[h, Range[numberofrandomwalks]]]]/
      Range[numberofrandomwalks];

The last step will be to plot the data obtained above, using the ListPlot built-in function. Now it is time to put together all our small pieces and write a code that will plot the estimates by using any sequence we specify as the input.

asian[seqinput_Integer]:=Module[{seq, data},
		seq = seqinput;
		vector[i_] := 
              Which[seq == 1, Table[Random[], {K}], seq == 2, halton[i, K], 
                seq == 3, scrhalton[i, K], seq == 4, faure[i, K]];
		nextprice[s_, innumber_] := 
              s + r*s*deltaT + sigma*s*Sqrt[deltaT]*zinv[innumber];
		h[j_Integer] := 
                Max[Apply[Plus, FoldList[nextprice, sini, vector[j]]]/
                      (K + 1) - exercise, 0];
		data = 
              Rest[FoldList[Plus, 0, Map[h, Range[numberofrandomwalks]]]]/
                Range[numberofrandomwalks];
		ListPlot[data, PlotJoined -> True, Frame -> True, 
              AxesOrigin -> {0, 0},  
              FrameLabel -> {"Number of random walks", "Estimate"}, 
              PlotLabel -> 
                Which[seq == 1, "Pseudorandom sequence", seq == 2, 
                  "Halton sequence", seq == 3, 
                  "Scrambled Halton sequence", seq == 4, 
                  "Faure sequence"] , 
              PlotStyle -> 
                Which[seq == 1, Hue[ .1], seq == 2, Hue[ .4], seq == 3, 
                  Hue[ .7], seq == 4, Hue[1]]]   ]

asian[1] plots the estimates when pseudorandom numbers are used.

graph1 = asian[1];

Similarly, asian[2] plots the estimates when the Halton sequence is used to generate the random walks.

graph2 = asian[2];

In the following graphs, the scrambled Halton and Faure sequences are used.

graph3 = asian[3];

graph4 = asian[4];

And here is the combined plot:

Show[graph1, graph2, graph3, graph4, 
      PlotLabel -> "Combined plot"];

There are no known analytic solutions for this problem. From the graphs it seems that the value of the option is between 5.2 and 5.6. To get a better understanding of the accuracy of the above estimates, we generated ten million random walks using pseudorandom numbers and obtained 5.429 as the estimate (this computation is done on a Digital Alpha 7640 using the Fortran programming language). We now plot this estimate as a horizontal line superimposed to the previous graph.

approx = 
        Plot[y = 5.429, {x, 1, 5000}, Frame -> True, 
          DisplayFunction -> Identity];

Show[graph1, graph2, graph3, graph4, approx, 
      PlotLabel -> "Combined plot", PlotRange -> {5, 5.8}];

The scrambled Halton sequence (blue) seems to give the best estimates, when the estimate of the ten million run is taken as a reference point. However, after 2000 random walks, the scrambled Halton and Faure (red) estimates become, relatively speaking, very close. They both provide better estimates, when compared to the pseudorandom (yellow) implementation. After a "poor" start, the Halton sequence (green) becomes competitive around 4000 random walks. [Editor's note: colors appear in the electronic version of this article.]

Notes

A detailed exposition of the general theory of Monte Carlo and quasi-Monte Carlo methods is given by Niederreiter [5]. Although low-discrepancy sequences provide more accurate estimates than Monte Carlo estimates in general, they have their own limitations. In recent research, several methods have been designed, in an attempt to combine the best features of Monte Carlo and quasi-Monte Carlo methods. These methods are usually refered as "hybrid-Monte Carlo methods". For applications of some hybrid-Monte Carlo methods to finance, see Ökten [6,7,8] and Owen [9].

References

[1] L. Kuipers and H. Niederreiter, Uniform Distribution of Sequences. New York: John Wiley & Sons (1974).
[2] H. Faure, Suggestions for Quasi-Monte Carlo Users, Lecture Notes in Pure and Applied Mathematics, Vol. 141, Marcel Dekker, Inc. (1993).
[3] A.G.Z. Kemna and A.C.F. Vorst, A Pricing Method for Options Based on Average Asset Values, Journal of Banking and Finance, 14, 113-129 (1990).
[4] M. Abramowitz and I.A. Stegun, Handbook of Mathematical Functions. New York:  Dover (1965).
[5] H. Niederreiter, Random Number Generation and Quasi-Monte Carlo Methods. Philadelphia CBMS-SIAM (1992).
[6] G. Ökten, Error Estimation for Quasi-Monte Carlo Methods, In Monte Carlo and Quasi-Monte Carlo Methods 1996 (Edited by H. Niederreiter et al.), pp. 353-368, Lecture Notes in Statistics #127, Springer-Verlag, New York (1998).
[7] G. Ökten, Applications of a Hybrid-Monte Carlo Sequence to Option Pricing, In Monte Carlo and Quasi-Monte Carlo Methods 1998 (Edited by J. Spanier and H. Niederreiter), Lecture Notes in Computational Science and Engineering, Springer-Verlag, New York, to appear.
[8]G. Ökten, Contributions to the Theory of Monte Carlo and Quasi-Monte Carlo Methods, Dissertation.com (1999) (www.dissertation.com/library/1120419a.htm).
[9] A.B. Owen and D.A. Tavella, Scrambled Nets for Value at Risk Calculations, In Understanding and Applying Value-at-Risk (Edited by S. Grayling), pp. 257-273, Risk Publications, London (1997).

ABOUT THE AUTHOR

ELECTRONIC SUBSCRIPTIONS

Included in the distribution for each electronic subscription is the file QuasiMonteCarlo.nb, containing Mathematica code for the material described in this article.