Multilevel Richardson-Romberg and Importance Sampling in Derivative Pricing

In [1], Giles explored the Richardson extrapolation in the context of both the MLMC and the standard Monte Carlo. The MLMC on its own was significantly better than the Richardson extrapolation. However, taken together, they worked even better. Lemaire and Pages [10] took this approach and undertook further comprehensive error analysis. They combined the method developed in [1] and Multistep Richardson extrapolation in order to minimize the cost of simulation.

Suppose, we are interested in estimating the expected payoff value i.e., $\mathbb{E}\left(P(X_{T})\right)$ in an option pricing problem, with finite time horizon $T>0$, where $(X_{t})_{0\leq t\leq T}$ is a process with values in $\mathbb{R}^{d}$, governed by the following SDE,

where, $W:=\begin{pmatrix}W_{1}&W_{2}&\dots&W_{q}\end{pmatrix}$ is a $q$-dimensional Brownian motion on a filtered probability space
$\left(\Omega,(\mathcal{F}_{t}\right)_{0\leq t\leq T},\mathbb{P})$, with $b:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ and $\sigma_{j}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ being the functions satisfying the following condition:

The assumption (A.1) ensures the existence and uniqueness [11] of the strong solution of (1). Now, in order to estimate $\mathbb{E}\left(P(X_{T})\right)$ one should be able to simulate $(X_{t})_{0\leq t\leq T}$. However, except for a handful of cases, where one can devise an analytical or semi-analytical solution, we must resort to discretization schemes to perform the simulation. As discussed above, we use the Euler-Maruyama discretization scheme in order to simulate equation (1).

Let, $M\in\mathbb{N}$ be the refinement factor and let $n_{l}=M^{l-1}$, where $l=1,\dots,L$ is the level refinement. Let $\displaystyle{h_{l}=\frac{T}{n_{l}}}$ be the time step size on level $l$. For example, the Euler scheme on level $l$ is given by,

where, $\eta_{n_{l}(t)}:=\lfloor t/h_{l}\rfloor h_{l}$. We approximate $\mathbb{E}\left[P\left(X_{T}\right)\right]$ by $\mathbb{E}\left[P\left(X_{T}^{n_{L}}\right)\right]$, where $L$ is the finest level. In [10] the authors proposed a methodology combining the order bias cancellation of Multilevel Richardson Romberg (ML2R), with variance control of MLMC, to solve the problem of improving the computational complexity, along with determining the optimal parameters in order to achieve the desired root-mean-squared error, with minimum computational effort. The study was based on weak error and strong error assumptions, on the sequence $\left(P(X_{T}^{h})\right)_{h\in\mathcal{H}}$, where $\mathcal{H}=\left\{\frac{\mathbf{h}}{n},n>1\right\}$. The assumptions [10, 11, 12] are stated as follows:

1. (A)

   Weak Error:

   $$\exists\alpha>0,\bar{L}\geq 1,(c_{l})_{1\leq l\leq\bar{L}},\quad\mathbb{E}[P(X^{h}_{T})]-\mathbb{E}[P(X^{0}_{T})]=\sum\limits_{k=1}^{\bar{L}}c_{k}h^{\alpha k}+h^{\alpha\bar{L}}\eta_{\bar{L}}(h).$$

   (2)

2. (B)

   Strong Error:

   $$\exists\beta>0,V_{1}\geq 1,\|P(X^{h}_{T})-P(X^{0}_{T})\|^{2}_{2}=\mathbb{E}\left[|P(X^{h}_{T})-P(X^{0}_{T})|^{2}\right]\leq V_{1}h^{\beta}.$$

   (3)

With the consideration of the above assumptions, the ML2R estimator is defined as,

where $\pi=(h,\mu,L)$ are the optimal parameters obtained as the solution to,

where the cost function is given by [12],

Further, in the above equation, $N_{l}$ is the number of sample paths on level $l$, and $\widetilde{W}_{l}$ are the weights given by,

where $\textbf{w}=\left(w_{l}\right)_{1\leq l\leq L}$ is the solution to the Vandermonde system $V\textbf{w}=e_{1}$, with the Vandermonde matrix being defined by,

Here, $\alpha$ is the weak error rate as defined above. The interested reader may refer to [10] for the construction of the optimal parameters, as the closed solution to equation (5). Here we tabulate the explicit values of these parameters required to achieve the root-mean-squared error $\epsilon$, with the following constants being used:

1. (A)

   $\displaystyle{\lambda=\sqrt{\frac{V_{1}}{\text{Var}\left(P(X^{0}_{T})\right)}}\leavevmode\nobreak\ \text{and}\leavevmode\nobreak\ \widetilde{c}_{\infty}=\lim_{L\rightarrow\infty}|c_{L}|^{1/\alpha}\in(0,\infty)}$.

2. (B)

   $\displaystyle{\underline{C}_{M,\beta}=\frac{1+M^{\beta/2}}{\sqrt{1+M^{-1}}}\leavevmode\nobreak\ \text{and}\leavevmode\nobreak\ \overline{C}_{M,\beta}=\left(1+M^{\beta/2}\right)\sqrt{1+M^{-1}}}$.

The above estimator was highly effective whenever there is a strong order of convergence i.e., $\beta\leq 1$, as it achieves $O\left(\epsilon^{-2}\log(1/\epsilon)\right)$ for $\beta=1$ and $O\left(\epsilon^{-2}e^{\frac{1-\beta}{\sqrt{\alpha}}\sqrt{2\log(1/\epsilon)\log(M)}}\right)$ for $\beta<1$, contrary to $O\left(\epsilon^{-2}\log(1/\epsilon)^{2}\right)$ and $O\left(\epsilon^{-2-\frac{1-\beta}{\alpha}}\right)$, respectively, which is achieved by the standard MLMC.

Now following the idea of [13], we consider a parametric family of stochastic process $\left(X_{t}(\theta)\right)_{0\leq t\leq T}$, with $\theta\in\mathbb{R}^{d}$, governed by the following SDE,

By the Girsanov’s Theorem, we know that there exists a probability measure $\mathbb{P}_{\theta}$ equivalent to $\mathbb{P}$ such that,

under which, the process $\left(\theta t+W_{t}\right)_{0\leq t\leq T}$ is a Brownian motion. Therefore,

The efficiency of the Monte Carlo method is considerably dependent on how small the variance is in the process of estimation. Among the many variance reduction techniques, available in literature (for improvement of efficiency of the Monte Carlo method), importance sampling is one such variance reduction technique, which is known for its efficiency [14]. In general, the idea of parametric importance sampling is to introduce the parametric transformation such that $\forall\theta\in\mathbb{R}^{d}$,

More specifically, we consider the general problem of estimating $\mathbb{E}[P(X)]$, where $X$ is a $d$-dimensional random variable. Also if $f(x)$ is the multivariate density, then,

where $\displaystyle{h(\theta,x)=\frac{P(x+\theta)f(x+\theta)}{f(x)}}$. Now the idea of importance sampling Monte Carlo method is to estimate $\mathbb{E}\left(P(X_{T})\right)$, where $\theta$ is given by,

In order to solve the above problem one can use the so-called Robbins-Monro algorithm that deals with a sequence of random variable $\left(\theta_{i}\right)_{i\in\mathbb{N}}$, which approximates $\theta^{*}$ accurately. However, the convergence of this algorithm requires restrictive conditions known as the non explosion condition given by [15],

To overcome these restrictive conditions, a truncation based procedure was introduced by Chen in [16, 17] which was furthered in [18, 19]. This was then proposed in the context of importance sampling [13]. An unconstrained procedure to approximate $\theta^{*}$, by using the regularity of the involved density extensively, was introduced in [20] along with the proof of the convergence of the algorithm.

The authors in [15] studied importance sampling in the context of statistical Romberg integration, which is a one-step generalization of the standard Monte Carlo integration. The paper used both constrained, as well as unconstrained optimization routines to approximate $\theta^{*}$. Further, they proved the almost sure convergence of both constrained and unconstrained versions of the optimization algorithm, towards the optimal parameter $\theta^{*}$, under discretized diffusion setting. The idea was then extended by [21] and [22] in the context of MLMC. The study carried out in [21] is the extension of [15] in the paradigm of MLMC. However [22] used sample average approximation to approximate $\theta^{*}$. Also, their algorithm performed an approximation procedure on each level of resolution, thus further enhancing the overall performance of the adaptive algorithm.

In this paper, we develop a hybrid algorithm incorporating ML2R and adaptive importance sampling procedure developed in [13] in the context of standard Monte Carlo in order to reduce the overall computational time, while achieving the desired accuracy. In [15], the authors studied the discretized version of the asymptotic variance obtained in the central limit theorem associated with the statistical Romberg method. The study was performed to determine the optimal parameter $\theta^{*}$, to minimize the variance of the Monte Carlo estimator. In [22], the authors extended the algorithm developed in [15] in the context of MLMC, building upon the Central Limit Theorem for the Euler MLMC studied in [23]. As discussed in [23], the optimal value of $\theta$ is obtained by using a constrained and unconstrained version of the stochastic Robbins-Monro algorithm, whereas in [22] it is determined using deterministic Newton-Raphson’s method. However, the algorithm was developed in the context of Euler discretization, and only European option pricing was studied, with no reference to path-dependent functionals. In this paper, we develop upon the Central Limit Theorem, for the general multilevel approach studied in [12], that incorporates the results for path-dependent functionals as well. We look at the stochastic optimization approach to approximate $\theta^{*}$. Also, we use the Milstein scheme, as well as Euler in order to discretize the underlying SDE. In order to highlight the novelty and contribution of this work, we present a comparison of the work accomplished, with the existing literature, in Table 2

In this subsection, we deal with the problem of finding the optimal value of $\theta_{l}$, to minimize the variance on level $l$. For the sake of brevity, we let,

Further let,

We start our discussion by recalling the Central Limit Theorems for ML2R from [12],

Based of the two theorems stated above we develop the optimization algorithm for $\beta>1$ and $\beta\in(0,1]$.

In this section we outline the stochastic approximation algorithm to estimate $\theta_{l}^{*}$. Based on the discussion in the previous section, we perform the discretization, such as Euler or Milstein, of the underlying SDE to approximate $Y_{l}$’s. Whereas for stochastic approximation we resort to the Robbins-Monro algorithm, to approximate the value of $\theta_{l}^{*},\leavevmode\nobreak\ l=1,\dots,L$. The aim is to construct a sequence of $(\theta_{l}^{n})_{n\in\mathbb{N}}$, such that, $\displaystyle{\lim_{n\rightarrow\infty}\theta_{l}^{n}=\theta_{l}^{*}}$.

Let $\Theta\subset\mathbb{R}^{q}$ be a compact convex subset such that $\theta\in\text{int}(\Theta)$. Then the sequence is constructed recursively as follows,

where $\displaystyle{\text{Proj}_{\Theta}(\theta)=\min_{\theta\in\Theta}\lvert\theta-\theta_{0}\rvert}$. Also $(\gamma_{n})_{n\geq 1}$ is a decreasing sequence of positive real numbers satisfying,

where for $\beta>1$,

and for $\beta\in(0,1]$,

It may be noted that in this paper, we use the constrained version of the Robbins-Monro algorithm in order to estimate the approximate value of $\theta_{l}^{*}$ for various level of resolutions.

Before, establishing the existence and uniqueness of $\theta_{l}^{*}$, we recall results related to the weights $(\widetilde{W}_{l})_{l=1,\dots,L}$, from [12] necessary for our study. Accordingly, let us define,

The following result was proved in [12],

With above results at our disposal, we present a series of results establishing the existence and uniqueness of the optimal parameter $\theta_{l}^{*}$ on various level of resolution. For the most part the proof follows line of argument similar to that presented in [22, 15, 13]. We start our discussion bu proving the following lemma, necessary for the existence and uniqueness results.