Importance Sampling for Rainbow Option Pricing

1. Introduction

Importance sampling is a general variance reduction technique in Monte Carlo simulation, and it has numerous applications in computational finance; see, for example, Glasserman et al. (1999, 2008), Glasserman (2004), Guasoni and Robertson (2008), and Wang (2012). In importance sampling, samples are simulated from alternative probability distributions, and an unbiased estimator is formed by multiplying each sample with an appropriate likelihood ratio. The literature on importance sampling is vast. Various methodologies and performance analyses can be found in this very partial list: Sigmund (1976), Heidelberger (1995), Sadowsky (1996), Dupuis and Wang (2004), Blanchet et al. (2007), Rubinstein and Kroese (2007), Bucklew (2010), and the references therein. Although most of the work on importance sampling is empirical, the goal of this paper is to build efficient importance sampling algorithms for the estimation of the price of high-dimensional rainbow options with rigorously provable asymptotically optimal schemes in the context of rare event simulation.

Almost all of the applications of importance sampling to computational finance are concerned with state-independent changes of measure. Within our context, state independence means that the change of measure can only depend on time. In other words, the alternative sampling distribution, which is usually obtained from solving an appropriate variational problem, remains fixed at any given time, regardless of the state of the system (Glasserman et al. 1999, Glasserman 2004, Guasoni and Robertson 2008). For such importance sampling schemes to attain efficiency, certain minimax conditions need to be imposed (see, e.g., Glasserman et al. 1999, condition 2.14 and Guasoni and Robertson 2008, condition 3.7). However, such minimax conditions will not hold in general, especially in the setting of rainbow options because many payoff functions are not convex. In Wang and Zhou (2015), a modified crossentropy scheme to select alternative sampling distributions from a class of mixtures for rainbow option simulation has been introduced. However, the scheme lacks theoretical performance analysis, and it requires pilot samples that may incur nontrivial computational overhead in higher dimensions.

In a series of papers, Dupuis and Wang (2004, 2007, 2009) have established that the change of measure needs to be state dependent (i.e., the change of measure can depend on both time and system state) in order to achieve efficiency in general. Such schemes, named dynamic importance sampling (dynamic IS) schemes, are built from suitable subsolutions to relevant Hamilton–Jacobi–Bellman (HJB) equations, and they have found much success in the setting of queueing networks. The goal of this paper is to systematically construct efficient state-dependent importance sampling schemes for rainbow options. We will demonstrate that the framework of dynamic importance sampling through subsolutions can be adapted to the simulation of rainbow options as well. However, this approach has limited success because it requires the option payoff function to satisfy certain conditions (Section 4.3). When applicable, the resulting dynamic importance sampling schemes are explicitly calculable and are not difficult to implement. Moreover, rigorous performance analysis can be carried out to establish the efficiency or asymptotic optimality of such schemes.

The difficulty of applying this dynamic importance sampling to general rainbow options lies in the construction of appropriate subsolutions to the associated HJB equations. This is largely because of the complicated and high-dimensional nature of a rainbow option’s payoff function. In order to circumvent this difficulty, we observe that all of these HJB equations are derived under the classical large deviation embedding, which has played an important role in the design of importance sampling schemes (Glasserman et al. 1999, Dupuis and Wang 2007, Guasoni and Robertson 2008). In this paper, we introduce a novel alternative large deviations embedding, which will reduce the complexity of the problem considerably. In short, it becomes sufficient to only investigate the estimation of probabilities instead of general expected values. Consequently, the construction of the important sampling schemes becomes much easier. The resulting asymptotically optimal (under this alternative embedding) importance sampling schemes are still state dependent, but the dependency is much simpler.

The paper is organized as follows. Section 2 establishes the model for the financial market. In Section 3, we give an overview of importance sampling and asymptotic optimality as well as a literature review. Section 4 considers the traditional large deviations embedding and establishes the connection between subsolutions and asymptotic optimality. We also discuss the construction of subsolutions when the terminal conditions take a certain form. The alternative large deviations embedding and the resulting importance sampling scheme are introduced in Section 5. Examples of numerical simulations are given in Section 6. We conclude the paper with further remarks and heuristic discussions in Sections 7 and 8.

2. Market Model and Rainbow Options

Because we are only interested in option pricing, we assume that the probability space $(\mathrm{\Omega },\mathcal{F})$ is equipped with the risk-neutral probability measure P. Let r be the constant risk-free interest rate. Consider a financial market with d risky assets whose prices at time t are denoted by a column vector $S(t)≔{[S_1(t),\dots ,S_d(t)]}^{\prime }$. We will assume that the financial market is complete and that stock prices are geometric Brownian motions. More precisely, let $B≔\{B(t)={[B_1(t),\dots ,B_d(t)]}^{\prime }:t\ge 0\}$ be a d-dimensional Brownian motion with strictly positive definite covariance matrix $\mathrm{\Sigma }≔[{\mathrm{\Sigma }}_{ij}]$ of size $d\times d$ and ${\mathrm{\Sigma }}_{ii}=1$ for every $i=1,2,\dots ,d$. Let $\sigma ≔{[{\sigma }_1,\dots ,{\sigma }_d]}^{\prime }$ be a column vector of strictly positive constants that denotes the volatilities of stock prices. For $i=1,2,\dots ,d$, the stock price $S_i$ is modeled by the geometric Brownian motion

We are interested in estimating the price of a rainbow option with expiration date T and payoff $h(S(T))≔h(S_1(T),\dots ,S_d(T))$, where $h:{\mathbb{R}}_{+}^d\to [0,\infty )$ is a nonnegative function defined on the open positive orthant

The option price is simply the expected value $v≔E[e^{-rT}h(S(T))]$. Prices of most rainbow options do not admit explicit formulae. The goal of our paper is to construct efficient importance sampling schemes for simulating such option prices.

For future analysis, it is convenient to introduce the log stock price process $X(t)≔{[X_1(t),\dots ,X_d(t)]}^{\prime }$, where for each $i=1,2,\dots ,d$,

In particular, letting $\mathrm{\Gamma }$ be a matrix such that $\mathrm{\Gamma }{\mathrm{\Gamma }}^{\prime }≔\mathrm{\Sigma }$ (see Remark 1), we can rewrite $X(T)=A+C\mathrm{\Gamma }Z$, where

C is a $d\times d$ diagonal matrix with $C_{ii}={\sigma }_i$ for $i=1,\dots ,d$, and $Z≔{[Z_1,\dots ,Z_d]}^{\prime }$ is a d-dimensional normal random vector with distribution $N(0,T{I}_d)$. Therefore, the discounted payoff can be expressed as

where we have adopted the notation that $\exp \{x\}≔{[\exp \{x_1\},\dots ,\exp \{x_d\}]}^{\prime }$ for any vector $x≔{[x_1,\dots ,x_d]}^{\prime }\in {\mathbb{R}}^d$.

Finally, we would like to impose some mild technical conditions on the payoff function h, which are satisfied by nearly all rainbow options. Throughout the paper, we will assume that h is of linear growth; that is, there exist some nonnegative constants a and b such that

for every $s\in {\mathbb{R}}_{+}^d$. Under this assumption, the option price v is finite. We will assume that the set $D≔\{h>0\}$ is either open or closed and that there is some continuous nonnegative function $\overline{h}:{\mathbb{R}}_{+}^d\to [0,\infty )$ such that $\overline{h}=h$ on $\{h>0\}$ or equivalently, $h=\overline{h}{1}_{\{h>0\}}$. The choice of $\overline{h}$ is not essential, and we can choose $\overline{h}$ so that it also satisfies the growth Condition (3); indeed, one can otherwise redefine $\overline{h}$ to be the minimum of itself with $a+b\Vert s\Vert$.

Remark 1.

When $\mathrm{\Sigma }$ is strictly positive definite, there is a particularly convenient choice of $\mathrm{\Gamma }$ based on Cholesky factorization (Glasserman 2004). In this case, $\mathrm{\Gamma }$ is lower triangular, and its components are explicitly available.

Remark 2.

Even though our analysis is largely based on market models where the stock prices are geometric Brownian motions, this restriction is not essential. The idea is applicable to more general jump diffusion models. We will discuss this extension with more details in Section 7.

3. Overview of Importance Sampling

Importance sampling is a general variance reduction technique for Monte Carlo simulation. It is particularly effective in the context of rare event simulation. The idea is to generate samples from a different probability distribution and multiply each outcome with an appropriate likelihood ratio so that the resulting estimator is unbiased.

Hereon, we will equip the sample space $(\mathrm{\Omega },\mathcal{F})$ with probability measures other than the original risk-neutral probability measure P. The expected value with respect to P will simply be denoted by $E[·]$, whereas the expected value with respect to some other probability measure, say Q, will be denoted by $E^Q[·]$. If those alternative probability measures are indexed by some parameter $\epsilon >0$, say $Q=P_{\epsilon }$, then the expected value with respect to those probabilities will be denoted by $E^{\epsilon }[·]$ for simplicity when no confusion is incurred.

Suppose one is interested in estimating the expected value $v≔E[X]$ of some random variable X. Let Q be an alternative probability measure such that $P\ll Q$. Then,

In other words, importance sampling produces an unbiased estimator by generating samples from an alternative probability measure Q and multiplying each sample of X by the likelihood ratio $dP/dQ$. The key in the design of importance sampling algorithms is the choice of the alternative sampling distribution Q. When appropriately designed, importance sampling can be very efficient, reducing the variance of the estimator by orders of magnitude. It is particularly useful in the estimation of probabilities or expected values related to rare events.

A commonly used criterion for selecting the alternative sampling distribution is the so-called “asymptotic optimality” (Sigmund 1976, Heidelberger 1995). Its definition requires that the present estimation problem be embedded into a sequence of estimation problems, where certain types of large deviations convergence results hold. To fix ideas, suppose that we wish to estimate the expected values of nonnegative random variables $\{X_{\epsilon }\}$ indexed by $\epsilon >0$ (the original problem corresponds to some $\epsilon$; i.e., $X=X_{\epsilon }$ for some $\epsilon$). Denote $v_{\epsilon }≔E[X_{\epsilon }]$. We assume that a large deviation type of asymptotics exists; that is, there exists a constant $\gamma \in \mathbb{R}$ such that

Let $P_{\epsilon }$ be a new probability measure such that $P\ll P_{\epsilon }$. The importance sampling estimator denoted by ${\widehat{v}}_{\epsilon }$ is given by

where $dP/d{P}_{\epsilon }$ denotes the Radon–Nikodým derivative or the likelihood ratio. Because the estimator is unbiased under $P_{\epsilon }$, minimizing its variance is equivalent to minimizing its second moment. However, by the Jensen inequality and the unbiasedness of ${\widehat{v}}_{\epsilon }$, we have the lower bound

Definition 1

(Asymptotic Optimality). An importance sampling scheme or change of measure $\{P_{\epsilon }\}$ is said to be asymptotically optimal if this lower bound is achieved: that is, if

$$\underset{\epsilon \to 0}{\overline{\lim }}\epsilon \log E^{\epsilon }[{\widehat{v}}_{\epsilon }^2]\le 2\gamma .$$

Remark 3.

To put everything into the framework of asymptotic optimality, we need to embed the original estimation problem into a sequence of estimation problems. Such embedding is not unique, which is a useful observation that we will exploit later on to deal with certain rainbow options.

4. Efficient Importance Sampling

It has been established in Glasserman and Wang (1997) and Dupuis and Wang (2004, 2007) that state-independent changes of measure will not lead to efficient importance sampling schemes in general, whereas asymptotically optimal state-dependent schemes can be constructed from classical subsolutions to related Hamilton–Jacobi–Bellman equations. For the financial market model under consideration, the stock prices are driven by Brownian motion B. Under an alternative sampling distribution that we consider, the Brownian motion B will become a Brownian motion with drift. Therefore, we expect that asymptotically optimal changes of measure correspond to state-dependent drifts in general.

4.1. The Classical Large Deviation Embedding

Recall the definition of G in (2). The goal is to estimate the option price $v≔E[G(Z)]$, where $Z≔{[Z_1,\dots ,Z_d]}^{\prime }$ is a d-dimensional normal random vector with distribution $N(0,T{I}_d)$. Following the setup in Glasserman et al. (1999), define $F(z)≔\log G(z)$ for any $z\in {\mathbb{R}}^d$. Note that F takes values in $\mathbb{R}\cup \{-\infty \}$ because G can be zero. Now, fix any $\epsilon >0$. Define

The option price of interest corresponds to $\epsilon =1$. Under very mild conditions, one can establish the large deviations or exponential growth/decay rate of $v_{\epsilon }$. The following result is a slight generalization of the Varadhan integral lemma, and its proof can be found in Glasserman et al. (1999, lemma 2.2).

Lemma 1.

Suppose that $f:{\mathbb{R}}^d\to \mathbb{R}\cup \{-\infty \}$ is continuous and satisfies the growth condition $f(z)\le c_1+c_2\Vert z{\Vert }^2$ for some $c_2<1/2$ and for all $z\in D$, where $D\subseteq {\mathbb{R}}^d$ is either open or closed. Let Z be a d-dimensional normal random vector with distribution $N(0,T{I}_d)$. Then,

$$\underset{\epsilon \to 0}{\lim }\epsilon \log E[e^{f(\sqrt{\epsilon }Z)/\epsilon }{1}_D(\sqrt{\epsilon }Z)]=\underset{z\in D}{\sup }[f(z)-\Vert z{\Vert }^2/2T].$$

Under our assumptions, $h=\overline{h}{1}_{\{h>0\}}$ for some nonnegative, continuous function $\overline{h}:{\mathbb{R}}_{+}^d\to [0,\infty )$ that satisfies the linear growth Condition (3). Now, define $\overline{G}$ according to (2) similarly with h replaced by $\overline{h}$. Finally let $\overline{F}≔\log \overline{G}$. It is not difficult to check that $\overline{F}$ is continuous and satisfies the growth condition of Lemma 1; indeed, $\overline{F}(z)\le c_1+c_2\Vert z\Vert$ for some $c_1,c_2$ and every $z\in {\mathbb{R}}^d$. Because h and $\overline{h}$ agree on the set $\{h>0\}$, it follows that $\overline{G}=G$ and thus, $\overline{F}=F$ on the set $D≔\{F>-\infty \}$. Furthermore, because $\{h>0\}$ is either open or closed and $z\mapsto \exp \{A+C\mathrm{\Gamma }z\}$ is a continuous bijection from ${\mathbb{R}}^d$ to ${\mathbb{R}}_{+}^d$, D is either open or closed. Thus, we may apply Lemma 1 to obtain that

Finally, we observe that the supremum is always attained at some $z^{*}\in \overline{D}$ when F is continuous on $\overline{D}$ because of the linear growth condition of F.

4.2. HJB Equation and Subsolution

Even though Dupuis and Wang (2004, 2007) are mostly concerned with discrete time dynamics, the technique therein can be naturally adapted to the diffusion setting. We should omit these repetitive heuristic details of deriving the HJB equation and directly state the definition of a subsolution. To this end, define $W≔{\mathrm{\Gamma }}^{-1}B$. Then, W is a d-dimensional standard Brownian motion. Define the scaled state process indexed by $\epsilon$:

Then, the quantity of interest can be written as

We say a continuously differentiable function $V:[0,T]\times {\mathbb{R}}^d\to \mathbb{R}$ is a classical subsolution with bounded gradient if $\nabla V$ is uniformly bounded and if V satisfies the partial differential inequality (see Remark 5)

for any $\forall (t,x)\in [0,T)\times {\mathbb{R}}^d$ with terminal condition

Note that the term $\Vert u{\Vert }^2/2$ is the large deviation rate function for the standard d-dimensional normal distribution $N(0,I_d)$ and is closely associated with the large deviations of d-dimensional standard Brownian motion $\{W^{\epsilon }\}$ (Dembo and Zeitouni 1998, theorem 5.2.3). The change of measure corresponding to a classical subsolution V is given by Dupuis and Wang (2007):

More precisely, the alternative sampling distribution is defined by the probability measure $P_{\epsilon }$, where

Because $u^{*}$ is assumed to be bounded, the right-hand side indeed defines a probability measure. By the Girsanov theorem (Karatzas and Shreve 1991), the process $\overline{W}$ defined by

is a standard Brownian motion under $P_{\epsilon }$. The corresponding importance sampling estimate for $v_{\epsilon }$ is simply

The performance of this estimate is characterized by the following theorem, which is essentially a continuous-time version of Dupuis and Wang (2007, theorem 8.1). This is the key result in the subsolution approach to importance sampling.

Theorem 1.

Let $V:[0,T]\times {\mathbb{R}}^d\to \mathbb{R}$ be a classical subsolution to the HJB Equations (6) and (7) with bounded gradient. If V is also twice continuously differentiable with respect to x with uniformly bounded Hessian matrix, then

$$\underset{\epsilon \to 0}{\overline{\lim }}\epsilon \log E^{\epsilon }[{\widehat{v}}_{\epsilon }^2]\le -2V(0,0).$$

Proof.

Let ${\nabla }^{2}V≔[{\nabla }_{ij}^{2}V]$ denote the $d\times d$ Hessian matrix of V and C be a constant such that $|{\nabla }^2{V}_{ij}(t,x)|\le C$ for every $i,j=1,\dots ,d$ and $(t,x)\in [0,T]\times {\mathbb{R}}^d$. Define the process $M≔\{M_t:t\in [0,T]\}$ with

$$\begin{array}{ll}{M}_t\hfill & ≔\exp \{-\frac{2}{\epsilon }V(t,W_{\epsilon }(t))-\frac{1}{\sqrt{\epsilon }}{\displaystyle {\int }_0^t{u}^{*}}(s,W_{\epsilon }(s))\hspace{0.17em}dW(s)+\frac{1}{2\epsilon }{\displaystyle {\int }_0^t}\Vert u^{*}(s,W_{\epsilon }(s)){\Vert }^2\hspace{0.17em}ds-\mathit{\text{Cdt}}\}.\hfill \end{array}$$(10)

Because of the boundedness of all first-order derivatives of V, V is of linear growth in both t and x. Therefore, it is not difficult to see that $E[M_t^2]\le K$ for some constant K and all $t\in [0,T]$. Applying the Itô formula and using (8), we have

$$\begin{array}{ll}\frac{d{M}_t}{M_t}\hfill & =\frac{2}{\epsilon }\left[-\frac{\partial V}{\partial t}+\Vert \nabla V{\Vert }^2+\frac{1}{2}\Vert u^{*}{\Vert }^2+〈\nabla V,u^{*}〉\right]\hspace{0.17em}dt-\left(Cd+\text{tr}{\nabla }^{2}V\right)\hspace{0.17em}dt-\frac{1}{\sqrt{\epsilon }}\left[2\nabla V+u^{*}\right]dW(t)\hfill \\ \hfill & =\frac{2}{\epsilon }\left[-\frac{\partial V}{\partial t}+\frac{1}{2}\Vert \nabla V{\Vert }^2\right]\hspace{0.17em}dt-\left(Cd+\text{tr}{\nabla }^{2}V\right)\hspace{0.17em}dt-\frac{1}{\sqrt{\epsilon }}\left[2\nabla V+u^{*}\right]dW(t),\hfill \end{array}$$

with the understanding that all of the coefficients on the right-hand side are taking values at $(t,W_{\epsilon }(t))$. It follows that M is indeed a local supermartingale because the coefficients of dt are negative by the subsolution property of V and that $\Vert {\nabla }_{ii}V{\Vert }_{\infty }\le C$ for all $i=1,\dots ,d$. Because M is nonnegative, it is thus a true supermartingale. In particular, $E[M_0]\ge E[M_T]$. Combined with the terminal Condition (7), it leads to

$$\begin{array}{lll}{e}^{-2V(0,0)/\epsilon }\hfill & \ge e^{-\mathit{\text{CdT}}}E\left[e^{-2V(T,W_{\epsilon }(T))/\epsilon }\frac{dP}{d{P}_{\epsilon }}\right]\hfill & \hfill \\ \hfill & \ge e^{-\mathit{\text{CdT}}}E\left[e^{2F(W_{\epsilon }(T))/\epsilon }\frac{dP}{d{P}_{\epsilon }}\right]\hfill & \hfill \\ \hfill & =e^{-\mathit{\text{CdT}}}{E}^{\epsilon }\left[e^{2F(W_{\epsilon }(T))/\epsilon }{\left(\frac{dP}{d{P}_{\epsilon }}\right)}^2\right]=e^{-\mathit{\text{CdT}}}{E}^{\epsilon }[{\widehat{v}}_{\epsilon }^2].\hfill \end{array}$$

Now, taking logarithm on both sides, multiplying with $\epsilon$, and then, sending $\epsilon$ to zero, we complete the proof. □

Remark 4.

Theorem 1 asserts that if $V(0,0)=-\gamma$, then the importance sampling estimate is asymptotically optimal. In many practical applications, $V(0,0)$ can be chosen to be arbitrarily close to $-\gamma$ for a classical subsolution V.

Remark 5.

The equation should be viewed as an Isaacs equation for interpretation. The interesting point of this Isaacs equation is that it is equivalent to the HJB equation by a scaling of two, which is the intrinsic reason why state-dependent importance sampling schemes can attain asymptotic optimality. In this paper, we simply use the form of the equivalent HJB equation. Details of the Isaacs equation may be found in Dupuis and Wang (2004, 2007).

4.3. Construction of Classical Subsolutions

Because of Theorem 1, building an efficient importance sampling scheme amounts to constructing a classical subsolution V of (6) and (7) with $V(0,0)$ as close to $\gamma$ as possible. In many applications, such subsolutions can be obtained as a mollification of piecewise affine weak subsolutions (Dupuis and Wang 2007).

Throughout this section, we assume that F is continuous on $\overline{D}$, where $D≔\{F>-\infty \}$. Recall the definition of $\gamma$ in (5), and denote by $z^{*}\in \overline{D}$ the maximizing point. We start with the simplest case where F is concave. Define $V:[0,T]\times {\mathbb{R}}^d\to \mathbb{R}$ by

We claim that V is a classical subsolution. Indeed, it satisfies (6) with equality. As for the terminal Inequality (7), we observe that $z^{*}/T\in \partial F$ because $z^{*}$ attains the supremum of (5). Therefore, by concavity,

which is exactly (7). Finally, $V(0,0)=-F(z^{*})+\Vert z^{*}{\Vert }^2/2T=-\gamma$, and thus, the corresponding state-independent importance sampling scheme is asymptotically optimal.

A more common situation is that F is the maximum of a concave function. For many options, the corresponding F can be written in this form, such as multistrike call options. Assume that $F≔F_1\vee \cdots \vee F_m$, where $F_i$ is a concave function for each $i=1,\dots ,m$. Define for each i

where $z^i\in {\mathbb{R}}^d$ is the maximizer for (5) with $F_i$ in place of F. Note that each $V_i$ satisfies the HJB Equation (6) with equality and $V_i(T,x)\le -F_i(x)$. Define $V≔V_1\wedge \cdots \wedge V_m$. Then, V satisfies (6) wherever it is differentiable with $V(T,x)\le -F(x)$. Moreover,

However, we cannot apply Theorem 1 directly because V is not differentiable everywhere. To do so, one can mollify the (weak) subsolution V to obtain a classical subsolution. To fix ideas, let $\delta$ be a small positive number, and define

Then, $V^{\delta }$ is a continuously differentiable function, and

with

It is easy to verify because of the convexity of the mapping $x\mapsto \Vert x{\Vert }^2$ for $x\in {\mathbb{R}}^d$ that $V^{\delta }$ satisfies the subsolution property (6). Furthermore, by the definitions of $V^{\delta }$ and V,

and

In other words, $V^{\delta }$ is a classical subsolution, where $V(0,0)$ is only $2\delta \log m$ away from asymptotic optimality (see Remark 6). In practice, $\delta$ is chosen to be small to guarantee good numerical performance.

Remark 6.

To achieve optimality, one can let $\delta$ depend on $\epsilon$, and a result analogous to Theorem 1 can be obtained. More precisely, if $\delta \to 0$ and $\delta /\epsilon \to \infty$, then the corresponding importance sampling scheme is asymptotically optimal. The proof is verbatim to that of Theorem 1 except that the trace of $V^{\delta }$ is now bounded by $C/\delta$ for some constant $\delta$, which is the reason behind condition $\delta /\epsilon \to \infty$.

5. An Alternative Embedding and Importance Sampling

In this section, we investigate an alternative large deviations embedding and the resulting importance sampling scheme.

5.1. Motivation and Overview

The classical large deviation embedding introduced by (4) and the subsequent subsolution approach yield a general framework for constructing asymptotically optimal importance sampling schemes. We have shown that if the payoff of a rainbow option can be expressed as the maximum of concave functions, the construction of subsolutions can be quite straightforward. However, many rainbow options (even simple ones, such as basket call) do not belong to this category. For those options, it is unclear, particularly in high dimensions, how to find systematic ways to construct subsolutions that lead to asymptotically optimal importance sampling schemes.

This prompts us to introduce a large deviation embedding other than (4). A distinct advantage of this alternative embedding is that the resulting importance sampling scheme is very easy to construct because it does not depend on the specific form of the payoff function. Rather, the construction will only rely on the distance from the origin to the region where the payoff function is strictly positive. For this reason, such importance sampling schemes are said to be “universal” (Dupuis and Wang 2004). They are still state-dependent importance sampling schemes, but the dependence is minimal; see (20).

5.2. An Alternative Large Deviation Embedding

To describe this new approach, we recall the definition of G in (2) and that of $D≔\{G>0\}$. Consider the following large deviations scaling. For any $\epsilon >0$, define

where Z is a d-dimensional normal random vector with distribution $N(0,T{I}_d)$. The quantity of interest corresponds to $\epsilon =1$. Here, we have abused the notation as $v_{\epsilon }$ is also used in the classical embedding. But, there should be no ambiguity given the context. We also define the set (which actually corresponds to a binary call)for every $\epsilon >0$. The heuristic is that a good importance sampling change of measure for estimating $P(A_{\epsilon })$ should also be good for estimating $v_{\epsilon }$. Even though $A_{\epsilon }$ may not be convex itself or the union of convex sets, one can design universal importance sampling (universal IS) schemes for it (Dupuis and Wang 2007). We adopt the following mild assumptions throughout this section.

Assumption 1

. D is nonempty, and $\overline{D^{°}}=\overline{D}$, where $D^{°}$ denotes the interior of D.

Assumption 2.

G is continuous on $\overline{D}$.

Assumption 3.

For some positive constants a and b, $|G(z)|\le a{e}^{b\Vert z\Vert }$ for all $z\in {\mathbb{R}}^d$.

Assumption 4.

$\overline{D}$ does not contain the origin.

We would like to remark that Assumption 4 is imposed only to avoid the trivial case where importance sampling is not really necessary. It will not affect the theoretical results in any way.

Lemma 2.

Let Z be a d-dimensional normal random vector with distribution $N(0,T{I}_d)$. Then, we have

$$\gamma ≔\underset{\epsilon \to 0}{\lim }\epsilon \log v_{\epsilon }=\underset{\epsilon \to 0}{\lim }\epsilon \log P(Z\in A_{\epsilon })=-\frac{1}{2T}\underset{z\in \overline{D}}{\inf }\Vert z{\Vert }^2,$$(17)

where the infimum is achieved at some $z^{*}\in \partial D$.

Proof.

Note that the last equality in (17) follows from classic large deviations and Assumption 1. It is trivial that the infimum is attained at some $z^{*}\in \partial D$ because D is nonempty.

To analyze $v_{\epsilon }$, we fix an arbitrarily small $\delta >0$ and define the set $D_{\delta }≔\{z\in {\mathbb{R}}^d:G(z)>\delta \}$. Note that $D_{\delta }^{°}$ is nonempty for $\delta$ small enough because $D^{°}$ is nonempty. It follows that $v_{\epsilon }\ge \delta P(\sqrt{\epsilon }Z\in D_{\delta })$, and classical large deviations theory yields

$$\underset{\epsilon \to 0}{\underset{¯}{\lim }}\epsilon \log v_{\epsilon }\ge \underset{\epsilon \to 0}{\underset{¯}{\lim }}[\epsilon \log \delta +\epsilon \log P(\sqrt{\epsilon }Z\in D_{\delta })]\ge -\frac{1}{2T}\underset{z\in D_{\delta }^{°}}{\inf }\Vert z{\Vert }^2.$$

Now, letting $\delta \downarrow 0$, it is straightforward to show by Assumption 1 that $D_{\delta }^{°}\uparrow D^{°}$, and thus,

$$\underset{\epsilon \to 0}{\underset{¯}{\lim }}\epsilon \log v_{\epsilon }\ge -\frac{1}{2T}\underset{z\in D^{°}}{\inf }\Vert z{\Vert }^2=-\frac{1}{2T}\underset{z\in \overline{D}}{\inf }\Vert z{\Vert }^2.$$

As for the other direction, take an arbitrarily large constant $M>0$, and observe that

$$G(\sqrt{\epsilon }Z)\le G(\sqrt{\epsilon }Z)1_{\{G(\sqrt{\epsilon }Z)>M\}}+M{1}_{A_{\epsilon }}(Z).$$(18)

It follows that

$$\begin{array}{ll}\underset{\epsilon \to 0}{\overline{\lim }}\epsilon \log v_{\epsilon }\hfill & \le \underset{\epsilon \to 0}{\overline{\lim }}\epsilon \log \left\{E[G(\sqrt{\epsilon }Z)1_{\{G(\sqrt{\epsilon }Z)>M\}}]+MP(Z\in A_{\epsilon })\right\}\hfill \\ \hfill & =\max \{\underset{\epsilon \to 0}{\overline{\lim }}\epsilon \log E[G(\sqrt{\epsilon }Z)1_{\{G(\sqrt{\epsilon }Z)>M\}}],\hfill \\ \hfill & \underset{\epsilon \to 0}{\overline{\lim }}\epsilon \log [MP(Z\in A_{\epsilon })]\}.\hfill \end{array}$$

Because the second term on the right-hand side equals ${\overline{\lim }}_{\epsilon \to 0}\epsilon \log P(Z\in A_{\epsilon })$, it remains to show that

$$\underset{M\to \infty }{\lim }\underset{\epsilon \to 0}{\overline{\lim }}\epsilon \log E[G(\sqrt{\epsilon }Z)1_{\{G(\sqrt{\epsilon }Z)>M\}}]=-\infty .$$

Because of the growth condition of G, it is sufficient to show that

$$\underset{M\to \infty }{\lim }\underset{\epsilon \to 0}{\overline{\lim }}\epsilon \log E[a{e}^{b\sqrt{\epsilon }\Vert Z\Vert }{1}_{\{a{e}^{b\sqrt{\epsilon }\Vert Z\Vert }>M\}}]=-\infty .$$

Observe that by the Hölder inequality,

$$E[a{e}^{b\sqrt{\epsilon }\Vert Z\Vert }{1}_{\{a{e}^{b\sqrt{\epsilon }\Vert Z\Vert }>M\}}]\le {\left\{E[a^2{e}^{2b\sqrt{\epsilon }\Vert Z\Vert }]P(a{e}^{b\sqrt{\epsilon }\Vert Z\Vert }>M)\right\}}^{1/2}.$$

By the dominated convergence theorem, $E[a^2{e}^{2b\sqrt{\epsilon }\Vert Z\Vert }]\to a^2{e}^{2b}$ as $\epsilon \to 0$. Furthermore, classical large deviations theory implies that

$$\begin{array}{ll}\underset{\epsilon \to 0}{\lim }\epsilon \log P(a{e}^{b\sqrt{\epsilon }\Vert Z\Vert }\ge M)\hfill & =-\frac{1}{2T}\inf \left\{\Vert z{\Vert }^2:a{e}^{b\Vert z\Vert }\ge M\right\}\hfill \\ \hfill & =-\frac{1}{2T{b}^2}{\left(\log \frac{M}{a}\right)}^2\to -\infty \hfill \end{array}$$

as $M\to \infty$. This completes the proof. □

5.3. Asymptotic Optimality

Below is the main result of the paper. It establishes the connection between the estimation of $v_{\epsilon }$ and that of $\{P(Z\in A_{\epsilon })\}$, and it confirms the heuristic that a good importance sampling scheme for the latter is also good for the former. In particular, if an importance sampling scheme is asymptotically optimal for estimating $\{P(Z\in A_{\epsilon })\}$, then it is also asymptotically optimal for estimating $\{v_{\epsilon }\}$ under some mild conditions.

Theorem 2.

Assume that $\{P_{\epsilon }\}$ is asymptotically optimal for the estimation of $\{P(Z\in A_{\epsilon })\}$. If for some $q>1$,

$$\underset{\epsilon \to 0}{\overline{\lim }}\epsilon \log E[1_{A_{\epsilon }}(Z){(dP/d{P}_{\epsilon })}^q]<\infty ,$$

then $\{P_{\epsilon }\}$ is asymptotically optimal for the estimation of $\{v_{\epsilon }\}$ as well.

Proof.

Denote the likelihood $L_{\epsilon }≔dP/d{P}_{\epsilon }$. The importance sampling estimator is just

$${\widehat{v}}_{\epsilon }=G(\sqrt{\epsilon }Z)L_{\epsilon },$$

where Z is sampled from $P_{\epsilon }$. Because of Lemma 2, it is sufficient to show that

$$\underset{\epsilon \to 0}{\overline{\lim }}\epsilon \log E^{\epsilon }[{\widehat{v}}_{\epsilon }^2]\le 2\gamma ,$$(19)

where $\gamma$ is defined to be the right-hand side of (17). Plugging in the form of ${\widehat{v}}_{\epsilon }$, the left-hand side (LHS) of (19) equals

$$\text{LHS}=\underset{\epsilon \to 0}{\lim }\epsilon \log E[{\widehat{v}}_{\epsilon }^2{L}_{\epsilon }^{-1}]=\underset{\epsilon \to 0}{\lim }\epsilon \log E[G^2(\sqrt{\epsilon }Z)L_{\epsilon }].$$

Take any large number M. It follows from (18) and ${(a+b)}^2\le 2a^2+2b^2$ that

$$G^2(\sqrt{\epsilon }Z)\le 2G^2(\sqrt{\epsilon }Z)1_{\{G(\sqrt{\epsilon }Z)>M\}}+2M^2{1}_{A_{\epsilon }}(Z).$$

By assumption, $\{P_{\epsilon }\}$ is asymptotically optimal for estimating $\{P(Z\in A_{\epsilon })\}$. It follows that

$$\begin{array}{lll}\underset{\epsilon \to 0}{\overline{\lim }}\epsilon \log E[2M^2{1}_{A_{\epsilon }}(Z)L_{\epsilon }]\hfill & =\underset{\epsilon \to 0}{\overline{\lim }}\epsilon \log E[1_{A_{\epsilon }}(Z)L_{\epsilon }]\hfill & \hfill \\ \hfill & =\underset{\epsilon \to 0}{\overline{\lim }}\epsilon \log E^{\epsilon }[1_{A_{\epsilon }}(Z)L_{\epsilon }^2] \hfill & =2\gamma .\hfill \end{array}$$

Now, for the other term, we can again use the Hölder inequality. Define $\alpha$ such that $1/\alpha +1/q=1$. Then,

$$\begin{array}{ll}\hfill & \underset{\epsilon \to 0}{\overline{\lim }}\epsilon \log E[2G^2(\sqrt{\epsilon }Z)1_{\{G(\sqrt{\epsilon }Z)>M\}}{L}_{\epsilon }]\hfill \\ \hfill & =\underset{\epsilon \to 0}{\overline{\lim }}\epsilon \log E[G^2(\sqrt{\epsilon }Z)1_{\{G(\sqrt{\epsilon }Z)>M\}}{1}_{A_{\epsilon }}(Z)L_{\epsilon }]\hfill \\ \hfill & \le \underset{\epsilon \to 0}{\overline{\lim }}\epsilon \log {\{E[G^{2\alpha }(\sqrt{\epsilon }Z)1_{\{G(\sqrt{\epsilon }Z)>M\}}]\}}^{1/\alpha }{\{E[1_{A_{\epsilon }}(Z)L_{\epsilon }^q]\}}^{1/q}\hfill \\ \hfill & \le \underset{\epsilon \to 0}{\overline{\lim }}\epsilon /\alpha \log E[G^{2\alpha }(\sqrt{\epsilon }Z)1_{\{G(\sqrt{\epsilon }Z)>M\}}]+\underset{\epsilon \to 0}{\overline{\lim }}\epsilon /q\log E[1_{A_{\epsilon }}(Z)L_{\epsilon }^q].\hfill \end{array}$$

Because of the assumption, the second term on the right-hand side cannot take a value of $\infty$. Thus, it remains to show that

$$\underset{M\to \infty }{\lim }\underset{\epsilon \to 0}{\overline{\lim }}\epsilon \log E[G^{2\alpha }(\sqrt{\epsilon }Z)1_{\{G(\sqrt{\epsilon }Z)>M\}}]=-\infty .$$

This is almost verbatim to the second half of the proof to Lemma 2, and we omit the details. Therefore, Inequality (19) follows readily. This completes the proof. □

5.4. Universal Importance Sampling

To apply Theorem 2, we need efficient importance sampling schemes for estimating $\{P(Z\in A_{\epsilon })\}$. Here, we investigate the universal schemes brought up in Dupuis and Wang (2007, section 10.7) and argue that it is indeed asymptotically optimal for the estimation of $\{P(Z\in A_{\epsilon })\}$ as well as $\{v_{\epsilon }\}$.

Recall $D≔\{G>0\}$ and $\gamma$ in (17). Fix any vector $\theta \in {\mathbb{R}}^d$ such that $\Vert \theta \Vert \le 1$, and define $u^{*}:[0,T)\times {\mathbb{R}}^d\to {\mathbb{R}}^d$ by

The importance sampling change of measure $\{P_{\epsilon }\}$ is again given by the Girsanov transform

Theorem 3.

$\{P_{\epsilon }\}$ is asymptotically optimal for the estimation of both $\{P(Z\in A_{\epsilon })\}$ and $\{v_{\epsilon }\}$.

Proof.

We first show that $\{P_{\epsilon }\}$ is asymptotically optimal for estimating $\{P(Z\in A_{\epsilon })\}$. This part is very similar to the proof of Theorem 1. Fix any $\delta >0$, and define $V:[0,T]\times {\mathbb{R}}^d\to \mathbb{R}$ by (see Dupuis and Wang 2007, section 10.7)

$$V(t,x)≔-\sqrt{\frac{-2\gamma }{T}(\Vert x{\Vert }^2+\delta )}-\gamma \left(1+\frac{t}{T}\right).$$

It is not difficult to show that the gradient and Hessian of V are uniformly bounded by some constant, say C. Analogous to the proof of Theorem 1, we define the process $M≔\{M_t:t\in [0,T]\}$ exactly as in (10). Because $u^{*}$ is bounded, we easily have $E[M_t^2]\le K$ for some constant K and all $t\in [0,T]$. Applying the Itô formula, we have

$$\begin{array}{ll}\frac{d{M}_t}{M_t}\hfill & =\frac{4\gamma }{\epsilon T}\left(\sqrt{\frac{\Vert \sqrt{\epsilon }{W}_{\epsilon }(t){\Vert }^2}{\Vert \sqrt{\epsilon }{W}_{\epsilon }(t){\Vert }^2+\delta }}-\frac{\Vert \sqrt{\epsilon }{W}_{\epsilon }(t){\Vert }^2}{\Vert \sqrt{\epsilon }{W}_{\epsilon }(t){\Vert }^2+\delta }\right)1_{\{W_{\epsilon }(t)\ne 0\}}\hspace{0.17em}dt\hfill \\ \hfill & +\frac{2\gamma }{\epsilon T}(1-\Vert \theta {\Vert }^2)1_{\{W_{\epsilon }(t)=0\}}\hspace{0.17em}dt-\left(Cd+\text{tr}{\nabla }^{2}V\right)\hspace{0.17em}dt-\frac{1}{\sqrt{\epsilon }}\left[2\nabla V+u^{*}\right]dW(t)\hfill \end{array}$$

with the understanding that all of the coefficients on the right-hand side are taking values at $(t,W_{\epsilon }(t))$. Clearly, the coefficients of dt are negative (note that $\gamma \le 0$). Define $F(x)≔\log 1_D(x)$. We observe that for $x\in \overline{D}$, $\Vert x{\Vert }^2\ge -2\gamma T$ by the definition of $\gamma$ in (17). Thus, for all $x\in \overline{D}$,

$$V(T,x)=-\sqrt{\frac{-2\gamma }{T}(\Vert x{\Vert }^2+\delta )}-2\gamma \le -\sqrt{\frac{-2\gamma }{T}}\Vert x\Vert -2\gamma \le 0.$$

In particular, we have $V(T,x)\le -F(x)$ for all $x\in {\mathbb{R}}^d$. The rest of the proof is almost verbatim to that of Theorem 1, and we arrive at

$$\underset{\epsilon \to 0}{\overline{\lim }}\epsilon \log E^{\epsilon }[1_{A_{\epsilon }}(Z){(dP/d{P}_{\epsilon })}^2]\le -2V(0,0)=2\gamma +2\sqrt{\frac{-2\gamma \delta }{T}}.$$

Letting $\delta \to 0$, we conclude that $\{P_{\epsilon }\}$ is asymptotically optimal for estimating $\{P(Z\in A_{\epsilon })\}$. To complete the proof, it is sufficient to show that for some $q>1$,

$$\underset{\epsilon \to 0}{\overline{\lim }}\epsilon \log E[1_{A_{\epsilon }}(Z){(dP/d{P}_{\epsilon })}^q]\le \underset{\epsilon \to 0}{\overline{\lim }}\epsilon \log E[{(dP/d{P}_{\epsilon })}^q]<\infty .$$(21)

We then invoke Theorem 2. This is trivial because $\Vert u^{*}\Vert$ is uniformly bounded. Thus,

$$\begin{array}{ll}E[{(dP/d{P}_{\epsilon })}^q]\hfill & =E\left[\exp \left\{\frac{-q}{\sqrt{\epsilon }}{\displaystyle {\int }_0^T{u}^{*}}(t,W_{\epsilon }(t))\hspace{0.17em}dW(t)+\frac{q}{2\epsilon }{\displaystyle {\int }_0^T}\Vert u^{*}(t,W_{\epsilon }(t)){\Vert }^2\hspace{0.17em}dt\right\}\right]\hfill \\ \hfill & \le E\left[\exp \left\{\frac{-q}{\sqrt{\epsilon }}{\displaystyle {\int }_0^T{u}^{*}}(t,W_{\epsilon }(t))\hspace{0.17em}dW(t)-\frac{q^2}{2\epsilon }{\displaystyle {\int }_0^T}\Vert u^{*}(t,W_{\epsilon }(t)){\Vert }^2\hspace{0.17em}dt\right\}\right]\hfill \\ \hfill & \times \exp \left\{\frac{(q^2+q)T\Vert u^{*}{\Vert }_{\infty }^2}{2\epsilon }\right\}=\exp \left\{\frac{(q^2+q)T\Vert u^{*}{\Vert }_{\infty }^2}{2\epsilon }\right\},\hfill \end{array}$$

and Inequality (21) follows readily. This completes the proof. □

Remark 7.

Both the dynamic importance sampling scheme outlined in Section 4 and the universal importance sampling schemes (20) are state dependent. However, the latter is much easier to construct because it only requires information on $\gamma$ or equivalently, the distance between the origin and the region D where the payoff function is strictly positive. Moreover, it is more broadly applicable because it does not require the payoff function to take any specific form. However, the price that we pay for this simplicity and generality is the efficiency. Even though the universal schemes are asymptotically optimal under the alternative large deviation scaling, they will not be as efficient as those schemes in Section 4 when the latter is applicable. This is quite obvious as the universal schemes do not utilize the details of the payoff functions. It is also demonstrated by the empirical evidence in the subsequent numerical experiments.

6. Numerical Experiments

In this section, we will present numerical examples on the pricing of rainbow options. For some options, we will use the dynamic importance sampling scheme outlined in Section 4, and for all options, we will also perform the universal importance sampling outlined in Section 5. In Tables 1–7, which show the numerical results, “R.E.” denotes the relative error of the importance sampling estimate (standard error divided by the estimate), and “var ratio” denotes the ratio of the variance of the plain Monte Carlo estimate to that of the importance sampling estimate. In all of the simulations, the time interval [0, T] is divided into $k=50$ equal-length subintervals (see Remark 8). Each estimate is based on n = 100,000 samples.

Table 1. Spread Call Option with Two Underlying Assets

All of the rainbow options that we consider satisfy the conditions in Section 2, including the linear growth Condition (3). Also, recall the definitions of vector A and square matrices C and $\mathrm{\Gamma }$ in (1) and thereafter. Note that $F(z)≔\log G(z)$, where G is defined in (2). When we apply universal IS, the set of interest is $D≔\{G>0\}$.

Example 1.

In this example, we consider a spread call option with $d=2$ underlying assets and payoff

$$h(S)≔{[S_1(T)-S_2(T)-K]}^{+},$$

where $K>0$ is a constant. In this case, for any $z\in {\mathbb{R}}^2$, the corresponding $F(z)$ and D are given by

$$\begin{array}{l}F(z)=\log {(\exp \{x_1\}-\exp \{x_2\}-K)}^{+},\hfill \\ D=\{z:\exp \{x_1\}-\exp \{x_2\}>K\},\hfill \end{array}$$

where $x≔{[x_1,x_2]}^{\prime }=A+C\mathrm{\Gamma }z$ for $z\in {\mathbb{R}}^2$ (Table 1).

1. Dynamic IS. One can show that F is a continuous concave function (we omit the technical details). Therefore, V in (11) is a classical subsolution with $z^{*}$ as the unique maximizer:

   $$z^{*}=\underset{z\in {\mathbb{R}}^d}{\text{arg\hspace{0.17em}max}}\left[F(z)-\frac{1}{2T}\Vert z{\Vert }^2\right].$$

   The corresponding change of measure is given by (8) and (9), where $u^{*}=-\nabla V(t,x)=z^{*}/T$. Note that $z^{*}$ must satisfy $\nabla F(z^{*})=z^{*}$.

2. Universal IS. The corresponding $\gamma$ is given by Lemma 2, where the minimizing z must lie on $\partial D=\{z:\exp \{x_1\}-\exp \{x_2\}=K\}$.

   We use the MATLAB function fminunc supplied with gradient (one can also use the iterative method in Glasserman et al. 1999) to solve for $z^{*}$. The equation of $\gamma$ can be converted into an equation of a single variable and solved by the bisection method. We omit the technical details and simply present the simulation results. For future reference, we also record $z^{*}$ and the maximizing $\widehat{z}$ for Lemma 2. The parameters are set at $r=0.05$, $T=1$, and

   $$S(0)=\left[\begin{array}{c}35\\ 30\end{array}\right], \left[\begin{array}{c}{\sigma }_1\\ {\sigma }_2\end{array}\right]=\left[\begin{array}{c}0.3\\ 0.4\end{array}\right], \mathrm{\Sigma }=\left[\begin{array}{cc}1& 0.2\\ 0.2& 1\end{array}\right].$$

Example 2.

In this example, we consider an outperformance digital call option with d underlying assets and payoff

$$h(S)=1_{\{S_1(T)\vee \cdots \vee S_d(T)\ge K\}},$$

where $K>0$ is a constant (Table 2). In this case, for any $z\in {\mathbb{R}}^d$, the corresponding $F(z)$ and D are given by

$$\begin{array}{l}F(z)=\log 1_{\{x_1\vee \cdots \vee x_d\ge \log K\}},\hfill \\ D=\{z:x_1\vee \cdots \vee x_d\ge \log K\},\hfill \end{array}$$

Table 2. Outperformance Digital Call Option with Three Underlying Assets

Table 2. Outperformance Digital Call Option with Three Underlying Assets

	$K=80$		$K=100$		$K=120$
	Dynamic	Universal	Dynamic	Universal	Dynamic	Universal
Estimate	3.41E-3	3.45E-3	2.41E-4	2.35E-4	2.05E-5	2.09E-5
R.E., %	0.64	2.97	0.65	3.40	0.68	5.44
Var ratio	71.2	3.2	1,012	38.3	9,256	138.9

---

Notes. R.E. denotes the relative error of the importance sampling estimate (standard error divided by the estimate). Var ratio denotes the ratio of the variance of the plain Monte Carlo estimate to that of the importance sampling estimate.

where $x≔{[x_1,\dots ,x_d]}^{\prime }=A+C\mathrm{\Gamma }z$.

1. Dynamic IS. $F≔F_1\vee \cdots \vee F_d$, where $F_i(z)≔\log 1_{\{x_i\ge \log K\}}$ is a concave function for each i. Therefore, $V≔V_1\wedge \cdots \wedge V_d$ in (11) is a (weak) subsolution, where for each i, $V_i$ is given by (12) with $z^i$ as the unique maximizer:

   $$z^i=\underset{z\in {\mathbb{R}}^d}{\text{arg\hspace{0.17em}max}}\left[F_i(z)-\frac{1}{2T}\Vert z{|}^2\right].$$

   Note that $F_i$ takes values of either zero or $-\infty$. Thus, the maximizer $z^i$ necessarily satisfies $F_i(z^i)=0$. From this observation, it is very straightforward to show that

   $$z^i=\frac{\log K-A_i}{\Vert {\theta }^i{\Vert }^2}·{\theta }^i,$$
   where ${\theta }^i$ denotes the ith column vector of the matrix ${(C\mathrm{\Gamma })}^{\prime }$. A classical subsolution $V^{\delta }$ can be obtained by mollification (13) with derivatives (14) and (15). The corresponding change of measure is given by (8) and (9), where $u^{*}=-\nabla V^{\delta }(t,x)$.

2. Universal IS. In this case, the corresponding $\gamma$ is given by Lemma 2, where obviously, the minimizing z must coincide with one of the $z^i$’s above, and it is easily identified as

   $$\gamma =-\frac{1}{2T}\underset{i=1,\dots ,d}{\min }\Vert z^i{\Vert }^2=-\frac{1}{2T}\underset{i=1,\dots ,d}{\min }\frac{{(\log K-A_i)}^2}{\Vert {\theta }^i{\Vert }^2}.$$

For simulation, we consider a digital call with $d=3$ underlying assets. The mollification parameter $\delta$ is set to be 0.1 in dynamic IS, and the other parameters are $r=0.05$, $T=1$, and

$$S(0)=\left[\begin{array}{c}40\\ 35\\ 40\end{array}\right], \left[\begin{array}{c}{\sigma }_1\\ {\sigma }_2\\ {\sigma }_3\end{array}\right]=\left[\begin{array}{c}0.2\\ 0.3\\ 0.1\end{array}\right], \mathrm{\Sigma }=\left[\begin{array}{ccc}1& 0.2& 0.3\\ 0.2& 1& -0.5\\ 0.3& -0.5& 1\end{array}\right].$$

Example 3.

In this example, we consider a multistrike call option with d underlying assets and payoff

$$h(S)≔\max \{S_1(T)-K_1,\dots ,S_d(T)-K_d,0\},$$

where $K>0$ is a constant (Table 3). In this case, for any $z\in {\mathbb{R}}^d$, the corresponding $F(z)$ and D are given by

$$\begin{array}{l}F(z)=\max \{\log {[\exp \{x_1\}-K_1]}^{+},\dots ,\log {[\exp \{x_d\}-K_d]}^{+}\},\hfill \\ D=\{z:\max \{x_1-\log K_1,\dots ,x_d-\log K_d\}>0\},\hfill \end{array}$$

Table 3. Multistrike Call Option with Four Underlying Assets

Table 3. Multistrike Call Option with Four Underlying Assets

	$\mathrm{\Delta }=20$		$\mathrm{\Delta }=30$		$\mathrm{\Delta }=40$
	Dynamic	Universal	Dynamic	Universal	Dynamic	Universal
Estimate	5.38E-2	5.45E-2	2.95E-3	3.01E-3	1.36E-4	1.38E-4
R.E., %	0.47	1.59	0.59	2.68	0.72	3.90
Var ratio	58.0	4.9	481.9	21.4	2,113	70.5

---

Notes. R.E. denotes the relative error of the importance sampling estimate (standard error divided by the estimate). Var ratio denotes the ratio of the variance of the plain Monte Carlo estimate to that of the importance sampling estimate.

where $x≔{[x_1,\dots ,x_d]}^{\prime }=A+C\mathrm{\Gamma }z$.

1. Dynamic IS. $F≔F_1\vee \cdots \vee F_d$, where $F_i(z)≔\log {[\exp \{x_i\}-K_i]}^{+}$ is a concave function for each i. Therefore, $V≔V_1\wedge \cdots \wedge V_d$ of (11) is a (weak) subsolution, where for each i, $V_i$ is given by (12) with $z^i$ as the unique maximizer:

   $$z^i=\underset{z\in {\mathbb{R}}^d}{\text{arg\hspace{0.17em}max}}\left[F_i(z)-\frac{1}{2T}\Vert z{\Vert }^2\right].$$

   It can be derived (we omit the technical details) that for each i, $z^i=x^i{\theta }^i$, where ${\theta }^i$ denotes the ith column vector of the matrix ${(C\mathrm{\Gamma })}^{\prime }$ and $x^i$ is the unique solution to the equation

   $$1+\frac{K_i}{e^{A_i+x\Vert {\theta }^i{\Vert }^2}-K_i}-\frac{x}{T}=0, x\ge \frac{\log K_i-A_i}{\Vert {\theta }^i{\Vert }^2}.$$

   Each $x^i$ can be obtained by the bisection method. A classical subsolution $V^{\delta }$ can be obtained by mollification (13) with derivatives (14) and (15). The corresponding change of measure is given by (8) and (9), where $u^{*}=-\nabla V^{\delta }(t,x)$.

2. Universal IS. The calculation of $\gamma$ is almost verbatim to that of Example 2 with K replaced by $K_i$ in each instance. Indeed,

   $$\gamma =\underset{i=1,\dots ,d}{\min }\frac{{(\log K_i-A_i)}^2}{\Vert {\theta }^i{\Vert }^2}.$$

   The mollification parameter $\delta$ is set to be 0.1 in dynamic IS. We consider $d=4$ underlying assets with parameters $r=0.05$, $T=1$, and

   $$\begin{array}{l}S(0)=\left[\begin{array}{c}40\\ 35\\ 30\\ 30\end{array}\right], \left[\begin{array}{c}{\sigma }_1\\ {\sigma }_2\\ {\sigma }_3\\ {\sigma }_4\end{array}\right]=\left[\begin{array}{c}0.1\\ 0.1\\ 0.2\\ 0.2\end{array}\right], \mathrm{\Sigma }=\left[\begin{array}{cccc}1& 0.2& 0.3& 0\\ 0.2& 1& 0.4& 0.2\\ 0.3& 0.4& 1& 0.3\\ 0& 0.2& 0.3& 1\end{array}\right],\\ \left[\begin{array}{c}{K}_1\\ K_2\\ K_3\\ K_4\end{array}\right]=S(0)+\mathrm{\Delta }\left[\begin{array}{c}1\\ 1\\ 1\\ 1\end{array}\right].\end{array}$$

   In order to investigate how the dimension of the problem affects the performance of these schemes, we also consider a group of $N=12$ underlying assets and multistrike call options associated with the first d assets, with $d=2,4,6,8,10,12$. The initial stock price $S(0)$, the volatility vector $\sigma ={[{\sigma }_1,{\sigma }_2,\dots ,{\sigma }_N]}^{\prime }$, and the covariance matrix $\mathrm{\Sigma }$ are randomly generated. Below is a realization of these parameters:

   $$\begin{array}{l}S(0)=\left[\begin{array}{c}18\\ 12\\ 22\\ 24\\ 18\\ 18\\ 12\\ 21\\ 30\\ 25\\ 15\\ 21\end{array}\right], \sigma =\left[\begin{array}{c}0.13\\ 0.15\\ 0.18\\ 0.13\\ 0.19\\ 0.17\\ 0.12\\ 0.21\\ 0.11\\ 0.14\\ 0.20\\ 0.26\end{array}\right], \left[\begin{array}{c}{K}_1\\ K_2\\ K_3\\ K_4\\ K_5\\ K_6\\ K_7\\ K_8\\ K_9\\ K_{10}\\ K_{11}\\ K_{12}\end{array}\right]=S(0)+30\left[\begin{array}{c}1\\ 1\\ 1\\ 1\\ 1\\ 1\\ 1\\ 1\\ 1\\ 1\\ 1\\ 1\end{array}\right],\\ \mathrm{\Sigma }=\left[\begin{array}{c}\begin{array}{cccccccccccc}\hfill 1.00\hfill & \hfill 0.08\hfill & \hfill 0.23\hfill & \hfill -0.15\hfill & \hfill -0.33\hfill & \hfill 0.32\hfill & \hfill 0.19\hfill & \hfill -0.31\hfill & \hfill 0.05\hfill & \hfill 0.12\hfill & \hfill 0.27\hfill & \hfill -0.13\hfill \\ \hfill 0.08\hfill & \hfill 1.00\hfill & \hfill 0.16\hfill & \hfill 0.11\hfill & \hfill 0.16\hfill & \hfill 0.24\hfill & \hfill 0.26\hfill & \hfill 0.16\hfill & \hfill 0.17\hfill & \hfill -0.31\hfill & \hfill 0.21\hfill & \hfill -0.12\hfill \\ \hfill 0.23\hfill & \hfill 0.16\hfill & \hfill 1.00\hfill & \hfill -0.35\hfill & \hfill -0.35\hfill & \hfill -0.18\hfill & \hfill 0.39\hfill & \hfill 0.03\hfill & \hfill -0.14\hfill & \hfill -0.36\hfill & \hfill -0.18\hfill & \hfill 0.27\hfill \\ \hfill -0.15\hfill & \hfill 0.11\hfill & \hfill -0.35\hfill & \hfill 1.00\hfill & \hfill 0.33\hfill & \hfill -0.20\hfill & \hfill -0.12\hfill & \hfill -0.07\hfill & \hfill 0.30\hfill & \hfill 0.22\hfill & \hfill -0.17\hfill & \hfill -0.35\hfill \\ \hfill -0.33\hfill & \hfill 0.16\hfill & \hfill -0.35\hfill & \hfill 0.33\hfill & \hfill 1.00\hfill & \hfill -0.30\hfill & \hfill -0.08\hfill & \hfill 0.28\hfill & \hfill 0.27\hfill & \hfill -0.09\hfill & \hfill 0.06\hfill & \hfill -0.14\hfill \\ \hfill 0.32\hfill & \hfill 0.24\hfill & \hfill -0.18\hfill & \hfill -0.20\hfill & \hfill -0.30\hfill & \hfill 1.00\hfill & \hfill -0.21\hfill & \hfill -0.38\hfill & \hfill 0.12\hfill & \hfill -0.01\hfill & \hfill 0.41\hfill & \hfill -0.17\hfill \\ \hfill 0.19\hfill & \hfill 0.26\hfill & \hfill 0.39\hfill & \hfill -0.12\hfill & \hfill -0.08\hfill & \hfill -0.21\hfill & \hfill 1.00\hfill & \hfill 0.06\hfill & \hfill -0.04\hfill & \hfill 0.05\hfill & \hfill -0.03\hfill & \hfill -0.13\hfill \end{array}\\ \begin{array}{cccccccccccc}\hfill -0.31\hfill & \hfill 0.16\hfill & \hfill 0.03\hfill & \hfill -0.07\hfill & \hfill 0.28\hfill & \hfill -0.38\hfill & \hfill 0.06\hfill & \hfill 1.00\hfill & \hfill 0.00\hfill & \hfill -0.40\hfill & \hfill -0.33\hfill & \hfill -0.01\hfill \\ \hfill 0.05\hfill & \hfill 0.17\hfill & \hfill -0.14\hfill & \hfill 0.30\hfill & \hfill 0.27\hfill & \hfill 0.12\hfill & \hfill -0.04\hfill & \hfill 0.00\hfill & \hfill 1.00\hfill & \hfill 0.01\hfill & \hfill -0.06\hfill & \hfill 0.02\hfill \\ \hfill 0.12\hfill & \hfill -0.31\hfill & \hfill -0.36\hfill & \hfill 0.22\hfill & \hfill -0.09\hfill & \hfill -0.01\hfill & \hfill 0.05\hfill & \hfill -0.40\hfill & \hfill 0.01\hfill & \hfill 1.00\hfill & \hfill -0.14\hfill & \hfill -0.37\hfill \\ \hfill 0.27\hfill & \hfill 0.21\hfill & \hfill -0.18\hfill & \hfill -0.17\hfill & \hfill 0.06\hfill & \hfill 0.41\hfill & \hfill -0.03\hfill & \hfill -0.33\hfill & \hfill -0.06\hfill & \hfill -0.14\hfill & \hfill 1.00\hfill & \hfill -0.01\hfill \\ \hfill -0.13\hfill & \hfill -0.12\hfill & \hfill 0.27\hfill & \hfill -0.35\hfill & \hfill -0.14\hfill & \hfill -0.17\hfill & \hfill -0.13\hfill & \hfill -0.01\hfill & \hfill 0.02\hfill & \hfill -0.37\hfill & \hfill -0.01\hfill & \hfill 1.000\hfill \end{array}\end{array}\right].\end{array}$$

   The other parameters in the simulation remain the same except that the sample size is increased to n = 400,000 (Table 4).

   Table 4. Multistrike Call Option with d Underlying Assets

   Table 4. Multistrike Call Option with d Underlying Assets

   	$d=2$		$d=4$		$d=6$		$d=8$		$d=10$		$d=12$
   	Dynamic	Universal	Dynamic	Universal	Dynamic	Universal	Dynamic	Universal	Dynamic	Universal	Dynamic	Universal
   Estimate	2.01E-13	1.96E-13	4.10E-6	3.73E-6	4.24E-6	4.60E-6	5.52E-5	5.21E-5	5.51E-5	5.42E-5	1.5E-3	1.6E-3
   R.E., %	0.31	2.21	0.24	3.44	0.35	6.38	0.29	7.70	0.28	9.95	0.28	8.07

   
   

   ---

   Note. R.E. denotes the relative error of the importance sampling estimate (standard error divided by the estimate).

   The variance ratios to the plain Monte Carlo estimates are not presented because it is often the case that the plain Monte Carlo cannot produce meaningful estimates/variances for $d=2,4,6,8,10$. For example, a typical plain Monte Carlo estimate for $d=2$ is zero, and a typical estimate for $d=8$ is 1.16E-5 with a $40\%$ relative error (in these cases, the true values of the expectation and variance are both vastly underestimated, a common occurrence in rare event simulation) (Glasserman and Wang 1997, Dupuis and Wang 2004). It is only for $d=12$ that the plain Monte Carlo produces meaning results, in which case the variance ratios are 2,300 and 4, respectively. In all of these simulations, one can see that both the dynamic IS and the universal IS outperform the plain Monte Carlo simulation. The dynamic IS, when feasible, is more efficient than the universal IS.

Example 4.

In this example, we show the advantage of universal IS by considering three different types of rainbow options. Let K, $\{c_i:i=1,\dots ,d\}$, and $\{K_i:i=1,\dots ,d\}$ be strictly positive constants:

1. basket option: $h(S)={[c_1{S}_1(T)+\cdots +c_d{S}_d(T)-K]}^{+}$;

2. pyramid option: $h(S)={[|S_1(T)-K_1|+\cdots +|S_d(T)-K_d|-K]}^{+}$; and

3. Madonna option: $h(S)={[\sqrt{|S_1(T)-K_1{|}^2+\cdots +|S_d(T)-K_d{|}^2}-K]}^{+}$.

For all of these options, the corresponding F cannot be expressed as the maximum of concave functions. Therefore, it is not clear how one can construct appropriate subsolutions and consequently, dynamic IS algorithms.

However, the universal IS algorithm is still easily applicable. For all simulations, we consider $d=3$ underlying assets with parameters $r=0.05$, $T=1$, and

$$\begin{array}{l}S(0)=\left[\begin{array}{c}40\\ 35\\ 30\end{array}\right], \left[\begin{array}{c}{\sigma }_1\\ {\sigma }_2\\ {\sigma }_3\end{array}\right]=\left[\begin{array}{c}0.2\\ 0.2\\ 0.1\end{array}\right], \mathrm{\Sigma }=\left[\begin{array}{ccc}1& 0.2& 0.3\\ 0.2& 1& 0.4\\ 0.3& 0.4& 1\end{array}\right],\\ \left[\begin{array}{c}{c}_1\\ c_2\\ c_3\end{array}\right]=\left[\begin{array}{c}0.3\\ 0.3\\ 0.4\end{array}\right], \left[\begin{array}{c}{K}_1\\ K_2\\ K_3\end{array}\right]=\left[\begin{array}{c}35\\ 35\\ 35\end{array}\right].\end{array}$$

We use the MATLAB function fmincon to solve the minimization problem in Lemma 2. There is no guarantee that the local minimum found by the function is the global minimum because of the lack of concavity (even though we have used random initial search conditions for fmincov to alleviate the problem to some degree). However, the simulation shows great improvement in performance, which is our ultimate goal. In Tables 5–7, “plain MC” stands for the plain Monte Carlo scheme.

Table 5. Basket Call Option with Three Underlying Assets

Table 6. Pyramid Call Option with Three Underlying Assets

Table 7. Madonna Call Option with Three Underlying Assets

6.1. Summary

The numerical results have shown that both dynamic IS and universal IS greatly improve the efficiency of the Monte Carlo simulation, especially when the options are deep out of the money. When applicable, the dynamic IS outperforms the universal IS, which is not surprising considering that the universal scheme only uses very limited knowledge of the model. On the other hand, universal IS (still provably asymptotically optimal in the alternative scaling) leads to great variance reduction compared with the plain Monte Carlo algorithm, and it is much more broadly applicable. Further, it also is much easier to construct for various rainbow options with complicated payoffs.

Remark 8.

The implementation of the dynamic importance sampling or the universal importance sampling requires the partition of the time interval. Even though we have picked a fairly large partition with $k=50$ here, the results remain nearly identical even if k is much smaller: for example, $k=10$. Note that in this paper, the plain Monte Carlo simulation does not require the partition of the time interval, whereas the state-dependent importance sampling does. Such computational overhead is insignificant compared with the reduction in variance, especially for the simulation of deep out-of-money options.

7. Further Remarks

This paper discusses two state-dependent importance sampling schemes: the “dynamic IS” based on the classical large deviations embedding and the “universal IS” based on an alternative large deviations embedding. Both of them greatly improve the efficiency of Monte Carlo simulation of rainbow options. Even though the theory is based on market models with lognormal stock prices and path-independent options, in principle one may be able to extended the algorithms to more general market models and path-dependent options. For the dynamic IS, the challenge lies in the overly complicated structure of the HJB equations, which is inevitable when one deviates from the simple lognormal stock price model and/or introduces path dependency, as well as the construction of appropriate subsolutions. On the other hand, the universal IS will reduce the complexity of the problem considerably by only dealing with probabilities. Even so, it still leads to a number of challenging yet interesting future research directions.

For example, some commonly used stock price models are jump diffusion models. The building blocks in those models, for instance, can take the form

where W is a Brownian motion, N is a Poisson process, and $\{Y_1,Y_2,\dots \}$ determine jump sizes. It would be interesting to construct an explicit and asymptotically efficient algorithm similar to (20) by analyzing the large deviation rate functions associated with those building blocks. Another example is path-dependent options, such as barrier options and Asian options. Regardless of whether such options are discretely or continuously monitored, one has to deal with a very high-dimensional or infinite-dimensional system. The construction and analysis of universal schemes within this context are not immediately clear. The last example that we would like to mention is with respect to American options. To the best of our knowledge, there has been no work on importance sampling for deep-out-of-the-money American options with provable asymptotic optimality. It would be interesting to see if any rigorous analysis can be done under the alternative large deviations embedding.

8. Heuristic Connection Between Two Embeddings

Finally, we would like to discuss the connection between the two large deviations embedding heuristically. It is difficult to give any rigorous justification, and hence, we should only proceed heuristically. The first large deviation scaling (4) leads to the rate function given by (5), whereas the second scaling (16) leads to (17). Intuitively speaking, the most important point for importance sampling under each large deviations scaling is the optimizing point in (5) and (17). Our goal is to heuristically justify that these two optimizing points are indeed close to each other in practice. In the special case where G is an indicator function (i.e., binary options), these two points are identical. More generally, let us assume that $G(z)={[f(z)-K]}^{+}$ for some continuous function f (i.e., various call options). We also let $T=1$ for notational convenience. In this case, (5) is just

where $D≔\{f(z)>K\}=\{G(z)>0\}$. Assume that f is continuously differentiable and $c_{1}f(z)\le \Vert \nabla f(z)\Vert \le c_{2}f(z)$ for some constants $c_1,c_2>0$ on the set $\overline{D}$. Given the form of the stock prices, this assumption is very natural. The maximizing $z^{*}\in \overline{D}$ should satisfy the equationand thus,

Because we are interested in the case where K is large and $\Vert z\Vert \to \infty$ as $K\to \infty$ if $z\in \overline{D}$ (because of the continuity of f), we can assume that $\Vert z^{*}\Vert >c_2$, and thus, the preceding inequality amounts to

Thus, the maximization problem for the first large deviations scaling or Equation (22) can be taken over all z such that

For all such z (note that ǁzǁis large when K is large), we have

where $\varphi (z)=o(\Vert z{\Vert }^2)$ because

Therefore, intuitively, the optimization problem is more or less equivalent to the following:

However, recalling the definition of $D^{\prime }$ in (23) and that for each $i=1,2$,

we conclude that $D^{\prime }$ is very close to $\partial D$. This justifies that the optimizer for (5) is very close to that of (17).