Journal of Risk

Risk.net

The importance of being scrambled: supercharged quasi-Monte Carlo

Sergei Kucherenko and Julien Hok

  • We show that randomized quasi-Monte Carlo and standard quasi-Monte Carlo methods outperform the Monte Carlo method in finance.
  • Randomized quasi-Monte Carlo outperforms standard quasi-Monte Carlo showing faster and more stable convergence.
  • Randomized quasi-Monte Carlo allows us to compute confidence intervals, providing a practical error bound.
  • We find that Owen’s scrambling outperforms other randomized quasi-Monte Carlo methods such as digital shift.

In many financial applications, the quasi-Monte Carlo (QMC) method based on Sobolʹ low-discrepancy sequences (LDSs) outperforms the Monte Carlo (MC) method, showing faster and more stable convergence. However, unlike Monte Carlo, QMC lacks a practical error estimate. The randomized QMC (RQMC) method combines the best of both methods. The application of scrambled LDSs allows us to compute the confidence intervals around the estimated value, providing a practical error bound. The randomization of Sobolʹ LDSs is done via two methods: Owen’s scrambling and digital shift. These are then compared considering the computation of Asian options and Greeks using the hyperbolic local volatility model. RQMC demonstrates a superior performance to the standard QMC, showing increased convergence rates and providing practical error bounds around the estimated values. The efficiency of the RQMC method strongly depends on the scrambling methods used. We recommend using Sobolʹ LDSs with Owen’s scrambling. The application of effective dimension-reduction techniques such as the Brownian bridge or the principal component analysis is critical in order to dramatically improve the efficiency of QMC and RQMC methods based on Sobolʹ LDSs.

1 Introduction

Monte Carlo simulation is a unique universal method widely used in the valuation of complex financial instruments and risk management engines. Its O(1/N) convergence rate does not depend on the number of dimensions d, although it is rather slow. Here, N is the number of sampled points or the number of paths. The Monte Carlo method also provides practical error estimates through the computation of confidence intervals.

Unlike random numbers (known to have bad uniformity properties), on which Monte Carlo simulation is based, deterministic low-discrepancy sequences (LDSs) are designed to fill multidimensional space as uniformly as possible, resulting in the significantly improved convergence rate of the quasi-Monte Carlo (QMC) method based on these sequences. Asymptotically, the QMC method is O(1/N), which is much higher than the convergence rate of the Monte Carlo. However, this theoretical estimate depends on the dimensionality of the problem, and in practice, the number of paths N required to achieve a given standard error may not be any lower than that of the Monte Carlo. For problems in quantitative finance, the number of dimensions d can reach many thousands. This has led to the misconception that QMC is not efficient at high dimensions. In practice, the effectiveness of QMC depends not on the nominal dimension d but on the so-called effective dimensions. There are many problems for which the d-dimensional function f is dominated by the first few variables (the low effective dimension in the truncation sense) or can be well approximated by the sum of low-dimensional functions in the analysis of variance (ANOVA) decomposition (the low effective dimension in the superposition sense). Such functions are very efficiently integrated by QMC, achieving a convergence rate close to O(1/N). It has been shown that problems in finance typically have low effective dimensions, and that effective dimensions can be reduced by applying special sampling schemes.

One of the main drawbacks of the QMC method is that since LDSs are deterministic, there is no statistical method for computing the standard error of the estimate. Specifically, this means that there is no clear termination criterion for stopping simulation after the required tolerance has been reached. There are techniques, known as randomized QMC (RQMC) methods, that introduce appropriate randomizations into the construction of LDSs. This allows for the measurement of integration errors through a confidence interval, in a similar manner to that for the Monte Carlo method, while preserving and often improving the convergence rate of the QMC. While the superior performance of the QMC methods based on Sobol LDSs has been widely studied and made use of, there are only a handful of purely academic papers in which RQMC methods are applied to problems in finance. It has been shown that RQMC offers both enhanced efficiency in comparison with pure QMC and the ability to produce confidence intervals; however, it is still not widely used in finance by practitioners. This paper aims to bridge this gap. We consider two popular methods of LDS randomization: random digital shift (DS) and Owen’s scrambling.

The remainder of this paper is organized as follows. Section 2 provides a brief review of the Monte Carlo, QMC and RQMC methods. Section 3 introduces the ANOVA decomposition and the concept of effective dimensions. Section 4 presents the time-homogeneous hyperbolic local volatility (HLV) model. Section 5 considers the two different time-discretization schemes. Section 6 discusses the Monte Carlo simulation of option pricing and Greeks. Section 7 presents our numerical results. Section 8 states our conclusions.

2 Monte Carlo, quasi-Monte Carlo and randomized quasi-Monte Carlo methods

The Monte Carlo method solves a problem by simulating the underlying process and then calculating the average result. It can be formulated as the computation of the multidimensional integral

  I[f]=Hdf(X)dX.   (2.1)

Here, function f(X) is integrable in the d-dimensional unit hypercube Hd. The Monte Carlo quadrature formula is based on the probabilistic interpretation of an integral as an expectation. The standard Monte Carlo estimator of the expectation is

  μ^N=1Ni=1Nf(Xi),   (2.2)

where {Xi} is a sequence of length N of random points uniformly distributed in Hd. The approximation μ^N converges to I[f] with probability 1. An integration error according to the central limit theorem is σ(f)/N, where σ2(f) is the function variance. Although σ2(f) is typically unknown, an unbiased estimate of it can be obtained along with the confidence intervals (Table 1). The convergence rate of the Monte Carlo method does not depend on the number of variables d, but it is rather slow. It is well known that random number sampling is prone to clustering. As new points are added randomly, they do not necessarily fill the gaps between the already sampled points.

In the classical QMC method, independent random points {Xi} are replaced by a deterministic set of points, such as the LDS, which are designed to cover the unit hypercube more uniformly than random points. Successive LDS points “know” about the position of previously sampled points and “fill” the gaps between them. The QMC algorithm for the evaluation of the integral (2.1) has a form similar to (2.2), but the LDS points {Qi}, QiHd, i=1,,N, are used instead of the random points {Xi}. The Sobol LDSs, also known as the digital (t,d) sequences in base 2, are the most well known and widely used LDSs in finance due to their efficiency (Glasserman 2003).

The efficiency of a particular Sobol LDS generator depends on the so-called direction numbers. In this paper, we use BRODA’s SobolSeq generator (BRODA 2023). The Sobol sequences produced by BRODA’s SobolSeq satisfy two additional uniformity properties: property A for all dimensions (currently the maximum dimension d=131 072) and property A for adjacent dimensions. It has been shown by Sobol et al (2011) on a number of different tests that BRODA’s SobolSeq generators generally outperform the other LDS generators considered in their paper. These results were corroborated in other publications (see, for example, Renzitti et al 2020).

A major drawback of the QMC method is the lack of practical estimates of the integration error. A classical worst-case error bound for numerical integration by the QMC method is given by the Koksma–Hlawka inequality. Although this bound can be used to obtain asymptotic convergence rates, it is too conservative and complex for the computation of practical error estimates.

The RQMC method combines the accuracy of QMC with Monte Carlo-type error estimation. Consider a set of randomized replications {Vi} of {Qi}. In the RQMC method,

  1. (a)

    for a fixed i, each point Vi is uniformly distributed ViU[0,1]d; and

  2. (b)

    the point set {Vi}, i=1,,N, is an LDS with probability 1.

Consider a set of K-randomized replications {Vi}=Vik, k=1,,K. We define the kth RQMC estimator for (2.1) as

  μ^nk =1ni=1nf(Vik),   (2.3)
and the sample mean as
  μ¯n =1Kk=1Kμ^nk.   (2.4)

We note that μ^nk are independent and identically distributed; hence, the sample standard deviation of this estimator and the corresponding root mean square error (RMSE) can be computed in the same way as for the Monte Carlo simulation. Table 1 provides the RMSEs εMC and εRQMC for the Monte Carlo and RQMC, respectively. It also provides the expressions for confidence intervals. It is assumed that K is large enough that the sample mean μ¯n is normally distributed. Then, zδ denotes the (1-δ) quantile of the standard normal distribution with the cumulative distribution function F:F(zδ)=1-δ. For a 95% confidence interval, δ=0.05 and zδ/21.96. In the case of small values of K, the normal quantile should be replaced with the one from the Student t distribution on K-1 degrees of freedom.

Assuming f is of bounded variation, the RMSE is O(1/Kn(1-α)) with α>0. To obtain an accurate estimate of μ¯n, we have to take large n and small K to keep the cost nK at an acceptable level.

Table 1: MC and RQMC sample standard deviations σ, RMSEs ε and confidence intervals. [The total number of function evaluations N=nK.]
σMC=1(N-1)i=1N(f(Xi)-μ^N)2 σRQMC=1(K-1)k=1K(μ^nk-μ¯n)2
εMC=σMCN εRQMC=σRQMCK
μ^N±zδ/2εMC μ¯n±zδ/2εRQMC

Owen’s nested scrambling achieves maximum randomization of the LDS while retaining multidimensional stratification (Owen 1997). Consider a b-ary expansion of an LDS point in base b:

  Qij=p=1mqi,pjb-p.   (2.5)

Here, Qij is the jth dimensional component of Qi, j=1,,d, i=1,,N, N=bm, b2, and the coefficients qi,pj{0,1,,b-1}. Owen’s scrambled version Vij of Qij is obtained by permuting the digits qi,pj in the following way:

  vi,1j=πj(qi,1j),vi,2j=π(qi,1j)j(qi,2j),vi,3j=π(qi,1j,qi,2j)j(qi,3j),  

and so on. All uniform random permutations πj over the set of {0,1,,b-1} are mutually independent, but each of them depends on the previous leading digits of Qij. Let M be the number of digits used in the binary number representation (M=32 or M=64). Then, the permutation tree in the d-dimensional case would consist of d(bM-1)/(b-1) permutations. For Sobol LDSs with b=2, M=32 and a low-dimensional problem with d=100, scrambling would require storing in memory around 4.3×1011 permutations. One way to reduce the computational costs would be to carry out the permutations for the first k bits only and then generate the other bits randomly. In this paper we use a modification of Owen’s scrambling with additional permutations (Atanassov and Kucherenko 2021). It has reduced memory and CPU requirements. Owen (1997) showed that for sufficiently smooth functions, εRQMCO(1/n(3/2-α)). It is n times higher than the best achievable rate O(1/n(1-α)) for the standard (nonscrambled) nets. This reduction arises from random error cancellations.

The random DS method is simple to implement and, unlike Owen’s scrambling, it does not impose extra memory requirements. For simplicity we present it for the Sobol sequence. Consider a set of d-dimensional Sobol points {Qi} in base b=2 (see (2.5)). Generate a random vector UU[0,1]d and produce a randomized version Vi of Qi with the components vi,pj=(qi,pjupj), i=1,,N, j=1,,d, p=1,,m. Here, upj is the pth digit in the binary representation of Uj. The symbol denotes the digital addition operation (a bitwise XOR operator). We note that K-randomized replicas of {Qi} are obtained with the same set of {Qi} and different Uk.

We note that there are other types of LDS randomization that are less costly than nested Owen’s scrambling and more efficient than DS. A survey of these methods is given in L’Ecuyer (2018). They all satisfy properties (a) and (b) of the RQMC method above, but do not possess the increased rate of convergence of Owen’s scrambling.

3 Analysis of variance decomposition and effective dimension

ANOVA decomposition can be used to explain the efficiency of QMC and RQMC methods in finance. Consider an integrable function f(X) defined in the unit hypercube Hd. It can be decomposed as

  f(X)=f0+i=1dfi(Xi)+i=1di<jdfij(Xi,Xj)++f1,2,,d(X1,X2,,Xd).   (3.1)

Each of the components fi1,,is(Xi1,,Xis) is a function of a unique subset of variables from X. The components fi(Xi) are called first-order terms, fij(Xi,Xj) second-order terms, and so on. Under appropriate regularity conditions, the decomposition is unique if

  01fi1,,is(Xi1,,Xis)dXik=0,1ks.   (3.2)

In this case terms are orthogonal with respect to integrations (Bianchetti et al 2015). For square integrable functions, the total variance of f can be decomposed as

  σ2=i=1dσi2+i=1di<jdσij2++σ1,2,,d2.   (3.3)

Here,

  σi1,,is2=01fi1,,is2(Xi1,,Xis)dXi1,,dXis  

are called partial variances.

Let |u| be a cardinality of a set of variables u. Define the Sobol indexes as Su=σu2/σ2. The effective dimension of f(X) in the superposition sense is the smallest integer ds such that

  0<|u|<dsSu1-ϵ  

with small ϵ0. If ds is close to 1, it means that f is well approximated by a sum of ds (or fewer) dimensional functions. There are cases where the first few inputs are much more important than the others. If

  u1,,dtSu1-ϵ,  

then f has an effective dimension dt in the truncation sense. Low effective dimension in the truncation sense can sometimes be achieved by redesigning the sampling scheme in such a way that the first few ANOVA components account for most of the variance in f (see Section 5.2 for details).

4 Time-homogeneous hyperbolic local volatility model

It is well known that implied volatility (the volatility input to the Black–Scholes formula that generates the market European call or put price) generally depends on the option’s strike price K and maturity T. When the implied volatility is plotted against the strike price, the resulting graph is typically downward-sloping for equity markets, often called “volatility skew”. For other markets, such as foreign exchange options or equity index options, where the graph typically is upward-sloping at both ends, the more familiar term “volatility smile” is used (see, for example, Gatheral (2012) for more details). For our numerical analysis, we consider the time-homogeneous HLV model, which better captures the market skew. It corresponds to a parametric local volatility-type model in which the dynamic of the underlying under the risk-neutral measure is

  dS(t)=rS(t)dt+σ~(S(t))dW(t),S0=1,   (4.1)

where r is the risk-free interest rate and

  σ~(S)=ν{(1-β+β2)βS+(β-1)β(S2+β2(1-S)2-β)}.   (4.2)

Here, ν>0 is the level of volatility, β(0,1] is the skew parameter and W is the standard Brownian motion. This model, which was introduced by Jäckel (2008), corresponds to the Black–Scholes model for β=1 and exhibits a skew for the implied volatility surface when β1. We note that the skew increases significantly with decreasing β. For example, with ν=0.3, β=0.2, the difference in volatility between strikes at 50% and at 100% is about 15%.

5 Time discretization schemes

5.1 Euler discretization of the stochastic differential equation

We consider the pricing of an option on a single asset whose value S(t) is defined by the stochastic differential equation (SDE) (4.1). To guarantee a positive price in the simulation, we use the transformation Y(t)=ln(S(t)). Then, from (4.1), we obtain

  dY(t)=[r-12σ2(Y(t))]dt+σ(Y(t))dWt,Y(0)=log(S(0)),   (5.1)

where σ(Y)=σ~(eY)/eY.

For a general Monte Carlo pricing framework with SDE discretization, we use the Euler–Maruyama scheme (Glasserman 2003; Kloeden and Platen 1992). In a discrete case of d equally distributed time steps, it has the following form:

  Yd(ti+1) =Yd(ti)+[r-12σ2(Yd(ti))](ti+1-ti)  
      +σ(Yd(ti))ti+1-ti(W(ti+1)-W(ti)),   (5.2)

with Yd(0)=log(S(0)), Δt=T/d, ti=iΔt, i=0,,d.

In addition to the statistical noise, there is a discretization error. Theorem 10.2.2 in Kloeden and Platen (1992) provides the conditions for the Euler–Maruyama scheme to have a strong error convergence of order 12. Under stronger conditions, as in Kloeden and Platen (1992, Theorem 14.5.2), the scheme reaches a weak error convergence of order 1.

5.2 Discretization of the Wiener process

There are different algorithms for the discretization of the Brownian motion W in (5.1). The standard (incremental) discretization algorithm follows directly from the definition of W(t). It is defined by the relation

  W(ti)=W(ti-1)+ΔtZi,1id,   (5.3)

where (Zi) are independent standard normal variates obtained from random numbers or from the Sobol LDSs using the inverse normal cumulative distribution function.

The Brownian bridge (BB) discretization is based on conditional distributions: the value of W(ti) is generated from the values of W(tl), W(tm), lim, at earlier and later time steps. This discretization first generates the Brownian motion at the terminal point

  W(T)=TZ1  

and then fills the other points using already found values of W(ti). The generalized BB formula is given by

  W(ti)=(1-γ)W(tl)+γW(tm)+γ(1-γ)(m-l)ΔtZi,   (5.4)

where γ=(i-l)/(m-l). It can be seen from (5.4) that the variance of the stochastic part of the BB formula γ(1-γ)(m-l)Δt decreases rapidly over successive levels of refinement, and the first few points contain most of the variance. Moreover, the variance in the stochastic part of (5.4) is smaller than that in (5.3) for the same time steps.

For the Monte Carlo method, the BB scheme has the same efficiency as the standard one, but it does affect the efficiency of the QMC based on Sobol LDSs. Sobol defined the “Sobol sequence” as the LPτ sequence. The τ-value is a quality parameter that measures the uniformity of the point sets. The smaller the τ-value is, the more uniformly distributed the points are. This value is equal to zero for one- and two-dimensional Sobol sequences only. In higher dimensions, as d increases, the smallest possible values of τ increase as well. Hence, the initial coordinates of Sobol LDSs are much better distributed than the later high-dimensional coordinates.

The BB discretization uses low well-distributed coordinates from each d-dimensional LDS vector point to determine most of the structure of a path and reserves the later coordinates to fill in fine details. In other words, well-distributed coordinates are used for important variables, while higher not-so-well-distributed coordinates are used for far less important variables. Thus, the BB sampling reduces the effective dimension in the truncation sense, which leads to the much higher convergence rate of the QMC algorithm for the majority (but not all) of the payoffs (Bianchetti et al 2015).

6 Monte Carlo simulation of option pricing and Greeks

6.1 Option pricing

We consider a geometric average Asian call option whose payoff function is given by

  PA=max(S¯-K,0),   (6.1)

where S¯ is a geometric average at d equally spaced time points:

  S¯=(i=1dSi)1/d,   (6.2)

where Si is the asset price at time ti=iT/d, 1id.

In a risk-neutral setting, the value of a call option with maturity T and strike price K is the discounted value of its payoff,

  AC(T,K)=e-rT𝔼[PA].   (6.3)

There is no analytical formula for (6.3) in the HLV model and it is estimated by the Monte Carlo method. First, we approximate the asset price S(ti) with Sd(ti)=eYd(ti) by discretizing the SDE (5.1) (as described in Section 5.1). Second, the expectation of the Asian payoff (6.1) is computed with the Monte Carlo estimator as an arithmetic average of payoffs taken over a finite number N of simulated price paths:

  ACN(T,K)=e-rT[1Ni=1Nmax(S¯(i)-K,0)],   (6.4)

where S¯(i) is an approximation of S¯ using the simulated ith price path.

6.2 Sensitivity factors

Sensitivity factors, or Greeks, are derivatives of the price AC(T,K) with respect to specific parameters such as spot price or volatility. They are computed for hedging and risk management purposes. In this paper, we focus only on the Delta defined as Δ=AC(T,K)/S0, where S0 is the current spot price. In the dynamic hedging, Delta corresponds to the number of assets we should hold for each option shorted to maintain a Delta-neutral position. As there is no analytical formula for the value of ACN(T,K), Delta can be estimated by Monte Carlo simulation and the finite-difference method. In the case of the central difference scheme, Delta is computed as

  ΔACN(T,K,S0+ϵs)-ACN(T,K,S0-ϵs)2ϵs,   (6.5)

where ϵs=hS0 is the increment and h is a shift parameter.

The path recycling of both pseudorandom sequences and LDSs is used to minimize the variance (as suggested, for example, by Glasserman (2003)). We note that the error analysis for Greeks is more complex than that for prices, since the variance of the Monte Carlo simulation is mixed with the bias due to the approximation of derivatives with finite differences. For the sensitivity factor estimation not to be entirely hidden by the Monte Carlo noise, in our computations the shift is chosen to be large enough; that is, h=0.01 (see Glasserman (2003) for detailed discussions).

7 Numerical results

In this section we present the results from the simulations of prices and the sensitivity factor Δ for Asian call options on a single underlying. The following parameters were used in the simulations: S0=100, r=3%, T=1, ν=30%, β=0.5 and the number of discrete time steps d=256. In the single underlying case, d corresponds to the problem dimensionality. We consider in-the-money (ITM), at-the-money (ATM) and out-of-the-money (OTM) options with the strike prices 80, 100 and 120, respectively, to investigate the effect of moneyness.

Numerical simulations using the Monte Carlo, QMC and RQMC methods were performed in order to compare the convergence of each method. The standard (incremental) discretization of Brownian motion was used in the Monte Carlo method. The BB algorithm was used in the QMC and RQMC methods (Section 5.2). The Mersenne Twister generator was used for Monte Carlo simulations and RQMC with DS. BRODA’s Sobol sequence generator with additional uniformity properties (described in Section 2) was used for QMC simulations (BRODA 2023; Sobol et al 2011). For QMC and RQMC to achieve an optimal uniformity sampling, the number of points n was taken to be powers of two. The reference values of prices and Deltas were obtained by the Monte Carlo method by averaging over K=10 independent runs, with each run using n=218 paths.

7.1 Price and Delta convergence

Asian call prices (a) ITM, (b) ATM and (c) OTM, with respect to the number of paths N=n, K=1 (on log 2 scale).
Figure 1: Asian call prices (a) ITM, (b) ATM and (c) OTM, with respect to the number of paths N=n, K=1 (on log 2 scale).
Asian call Deltas (a) ITM, (b) ATM and (c) OTM, with respect to the number of paths N=n, K=1 (on log 2 scale).
Figure 2: Asian call Deltas (a) ITM, (b) ATM and (c) OTM, with respect to the number of paths N=n, K=1 (on log 2 scale).

First, we analyze the convergence plots, that is, the values of the prices (Figure 1) and Deltas (Figure 2) versus the number of paths. One trial (K=1) is used for both Monte Carlo and RQMC runs. For simulated prices, RQMC with Owen’s scrambling outperforms all other methods for the ITM and ATM cases, converging faster to the reference levels. Its efficiency is followed by that of QMC and RQMS with DS methods, which outperform MC for the ITM and ATM cases. These two methods show similar performance. Similar but less pronounced trends are present in the case of OTM.

For simulated Deltas, RQMC with Owen’s scrambling marginally outperforms all other methods for the ITM case. It shows a similar performance to QMC for the ATM case. Both the RQMC and QMC methods slightly outperform Monte Carlo in the case of OTM.

We note that the step-like behavior of the convergence patterns in some QMC and RQMC plots is likely to be a result of the design of the Sobol sequence generators, as can be seen in Figure 2 of Sobol et al (2011, p. 70).

We can conclude that in all considered cases (ITM, ATM and OTM), the RQMC method with Owen’s scrambling shows a faster convergence than the other methods. The RQMC with DS gives a similar performance to QMC.

7.2 Confidence intervals

Table 2: εMC and εRQMC at K=10, n=212 (N=40 960) of price estimations.
  ITM ATM OTM
εMC 3.26×10-2 2.2×10-2 1.24×10-3
εRQMC (Owen) 3.67×10-4 6.09×10-4 4.04×10-4
εRQMC (DS) 7.12×10-4 1.03×10-3 4.59×10-4
εMC/εRQMC (Owen) 89 36 3
Table 3: εMC and εRQMC at K=10, n=212 (N=40 960) of Delta estimations.
  ITM ATM OTM
εMC 3.29×10-4 2.43×10-3 4.09×10-4
εRQMC (Owen) 3.23×10-5 5.98×10-4 1.41×10-4
εRQMC (DS) 4.04×10-5 4.34×10-4 1.46×10-4
εMC/εRQMC (Owen) 10 4 3

Tables 2 and 3 show the RMSE (Table 1) for the Monte Carlo and RQMC methods of price and Delta estimation, respectively. The ratios of Monte Carlo to RQMC with Owen’s scrambling error estimates show a dramatic improvement when using RQMC, with the largest improvement ratio being for ITM. The RQMC with Owen’s scrambling method is on average better than the RQMC with the DS method, producing a smaller εRQMC.

Asian call Deltas with 95% confidence intervals computed at K=10 for ITM call options for RQMC with (a) Owen and (b) DS methods, with respect to the number of paths n (on log 2 scale).
Figure 3: Asian call Deltas with 95% confidence intervals computed at K=10 for ITM call options for RQMC with (a) Owen and (b) DS methods, with respect to the number of paths n (on log 2 scale).

Figure 3 shows the prices and Deltas with confidence intervals versus the number of simulation paths n for the ITM call. Visually, the results for the ATM and OTM cases are similar (therefore they are not shown). As expected, the confidence intervals shorten as the number of paths increases. RQMC with Owen’s scrambling offers the most accurate results by significantly reducing the bounds of the confidence intervals.

7.3 Performance analysis

We also analyze the relative performances of the considered methods in terms of their convergence rates. For all considered sampling schemes the following power law for the integration error is observed empirically in numerical tests:

  εnCnα.   (7.1)

For the Monte Carlo method α=0.5. For applications of the QMC and RQMC methods to financial problems it is quite common that α>0.5. Its value can be very close to 1 irrespective of the nominal dimension when the effective dimensions are low.

In Section 2 we discussed the fact that because LDSs are deterministic, there are no statistical measures such as variances associated with them. Hence, the constant C in (7.1) is not a variance and (7.1) does not have a probabilistic interpretation. In practice, the RMSE for the Monte Carlo, RQMC and QMC methods for any fixed n can be estimated by computing the following error averaged over K independent runs:

  εn=1Kk=1K(V-Vn(k))2,   (7.2)

where V is the exact, or estimated, value of the integral (an option price or Delta in our case) obtained at a very large n, and Vn(k) is the simulated value for the kth run, performed using n paths.

We note the difference between the definition of εRQMC given in Table 1 (computed with respect to μ¯n) and the definition of εn given in (7.2) (computed with respect to V).

For Monte Carlo and RQMC, runs based on different seed points are statistically independent. In the case of QMC, different runs are obtained using nonoverlapping sections of the LDS. In our computations K=10.

Table 4: Extracted α in the QMC and RQMC (Owen’s scrambling) methods.
  ITM ATM OTM
QMC (price) 1.0 0.95 0.82
RQMC (price) 0.7 0.79 0.79
QMC (Delta) 0.65 0.64 0.71
RQMC (Delta) 0.66 0.62 0.62
RMSE for ITM Asian call prices (a) ITM, (b) ATM and (c) OTM, with respect to the number of paths n, K=10. Owen's scrambling is used in RQMC.
Figure 4: RMSE for ITM Asian call prices (a) ITM, (b) ATM and (c) OTM, with respect to the number of paths n, K=10. Owen’s scrambling is used in RQMC.
RMSE for ITM Asian call Deltas (a) ITM, (b) ATM and (c) OTM, with respect to the number of paths n, K=10. Owen's scrambling is used in RQMC.
Figure 5: RMSE for ITM Asian call Deltas (a) ITM, (b) ATM and (c) OTM, with respect to the number of paths n, K=10. Owen’s scrambling is used in RQMC.

Figures 4 and 5 show the RMSE versus the number of paths n for the Monte Carlo, QMC and RQMC methods on a log2log2 scale. We fitted the regression lines that follow the power law (7.1) to extract the convergence rates α: they are the slopes of the regression lines (Table 4). We also extracted the intercepts of the regression lines (log2(C) (7.1)) (though these are not presented here). These intercepts provide useful information about the efficiency of the QMC, RQMC and Monte Carlo methods: lower intercepts mean that the simulated value starts closer to the exact value.

As expected, for the Monte Carlo method, α=0.5 for all cases, although these results are not presented in Table 4. For the QMC and RQMC methods α>0.5 for all cases. It is higher for the prices than it is for the Deltas. Although it is marginally higher for QMC, the intercepts C of the regression lines are always lower for RQMC than for other methods for considered ranges of n. This makes RQMC the most efficient method (although the efficiencies of RQMC and QMC are similar for the OTM case).

8 Conclusions

We presented and discussed the results of an application of the Monte Carlo, QMC and RQMC methods for derivative pricing and risk analysis based on the HLV model. The results presented for the Asian option show the superior performance of the QMC and RQMC methods. RQMC not only increases the rate of convergence of the QMC but also allows for the computation of the confidence intervals around the estimated value. The efficiency of the RQMC method strongly depends on the scrambling methods used. We advise using Sobol LDSs with Owen’s scrambling as it is the most efficient method.

Declaration of interest

Sergei Kucherenko serves as a director for BRODA, UK. However, this potential conflict of interest did not influence the design, conduct or interpretation of the study. The authors are committed to presenting unbiased and objective findings, and all data analysis and interpretation was conducted independently. The authors alone are responsible for the content and writing of the paper.

References

  • Atanassov, E., and Kucherenko, S. (2021). Implementation of Owen’s scrambling with additional permutations for Sobol sequences. Report, BRODA Ltd., UK. URL: http://www.broda.co.uk/.
  • Bianchetti, M., Kucherenko, S., and Scoleri, S. (2015). Pricing and risk management with high-dimensional quasi-Monte Carlo and global sensitivity analysis. Wilmott Magazine 2015(78), 46–70 (https://doi.org/10.1002/wilm.10434).
  • BRODA (2023). High-dimensional Sobol sequence generators. Report, BRODA Ltd., UK. URL: http://www.broda.co.uk/.
  • Gatheral, J. (2012). The Volatility Surface: A Practitioner’s Guide. Wiley (https://doi.org/10.1002/9781119202073).
  • Glasserman, P. (2003). Monte Carlo Methods in Financial Engineering. Stochastic Modelling and Applied Probability, Volume 53. Springer (https://doi.org/10.1007/978-0-387-21617-1).
  • Jäckel, P. (2008). Hyperbolic local volatility. Working Paper, OTC Analytics. URL: http://www.jaeckel.org/HyperbolicLocalVolatility.pdf.
  • Kloeden, P. E., and Platen, E. (1992). Numerical Solution of Stochastic Differential Equations. Stochastic Modelling and Applied Probabiltiy, Volume 23. Springer (https://doi.org/10.1007/978-3-662-12616-5).
  • L’Ecuyer, P. (2018). Randomized quasi-Monte Carlo: an introduction for practitioners. In Monte Carlo and Quasi-Monte Carlo Methods, Owen, A. B., and Glynn, P. W. (eds), pp. 29–52. Springer Proceedings in Mathematics and Statistics, Volume 241. Springer (https://doi.org/10.1007/978-3-319-91436-7_2).
  • Owen, A. B. (1997). Scrambled net variance for integrals of smooth functions. Annals of Statistics 25(4), 1541–1562 (https://doi.org/10.1214/aos/1031594731).
  • Renzitti, S., Bastani, P., and Sivorot, S. (2020). Accelerating CVA and CVA sensitivities using quasi-Monte Carlo methods. Wilmott Magazine 2020(108), 78–93 (https://doi.org/10.2139/ssrn.3193219).
  • Sobol, I. M., Asotsky, D., Kreinin, A., and Kucherenko, S. (2011). Construction and comparison of high-dimensional Sobol generators. Wilmott Magazine 2011(56), 64–79 (https://doi.org/10.1002/wilm.10056).

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

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

You are currently unable to copy this content. Please contact info@risk.net to find out more.

You need to sign in to use this feature. If you don’t have a Risk.net account, please register for a trial.

Sign in
You are currently on corporate access.

To use this feature you will need an individual account. If you have one already please sign in.

Sign in.

Alternatively you can request an individual account here