ML2R estimator simulation

Multilevel Richardson Romberg estimator

We propose and analyze a new estimator called Multilevel Richardson-Romberg (ML2R), introduced by Lemaire & Pagès (see [ML2R]), which combines the higher order bias cancellation of the Multistep Richardson-Romberg method and the variance control resulting from the stratification in the Multilevel Monte Carlo (MLMC) method. The Multilevel Monte Carlo method has been extensively applied to various fields of numerical probability from last years. We refer to the webpage for a comprehensive list of references in this field.

Goal

We aim at computing a very accurate approximation of \(I_0 = \mathbb{E}[Y_0]\) by a Monte Carlo type estimator where the (non-degenerate) random variable \(Y_0 \in L^2(\mathbb{P})\) cannot be simulated (exactly) at a reasonable cost. Let \((\Omega, \mathcal{A}, \mathbb{P})\) be a probability space and \((Y_h)_{h \in \mathcal{H}}\) be a family of real-valued random variables in \(L^2(\mathbb{P})\) associated to \(Y_0\) where \(\mathcal{H} \subset (0, {\mathbf h}]\) is an admissible subset of parameters having \(0\) as a limiting value and such that \(\frac{\mathcal{H}}{n}\subset \mathcal{H}\) for every integer \(n\ge 1\). Typically the random variable \(Y_h\) results from a time discretization scheme of parameter \(h\) so that we will speak of \(h\) as the bias parameter in what follows.

Framework

Weak and strong error

It is classical background that the family \((Y_h)_{h \in \mathcal{H}}\) satisfies the weak and the strong error assumptions:

\[ {(\mathbf{Weak \, Err})_{\alpha}}\quad \equiv\quad \exists \,\alpha > 0, c_k, \bar R \ge 1, \quad \mathbb{E}[Y_h] - \mathbb{E}[Y_0] = \sum_{k=1}^{\bar R} c_k h^{\alpha k} + o(h^{ \alpha (\bar R+1)}), \]

\[ {(\mathbf{Strong \, Err})_{\beta}}\quad \equiv\quad \exists\, \beta > 0, V_1 \ge 0, \quad \|Y_h - Y_0\|_2^2 = \mathbb{E}[\left| Y_h - Y_0\right|^2] \le V_1 h^\beta \]

Computational cost

We assume that the computational cost of \(Y_h\) is inverse linear with respect to \(h\), \(i.e.\)

\[\mathrm{Cost}(Y_h) = \kappa h^{-1},\]

where \(\kappa\) is a constant depending on each specific machine.

Crude Monte Carlo estimator

The well known crude Monte Carlo estimator is

\[I^{N}_h = \frac{1}{N} \sum_{k=1}^N Y^{(k)}_h\]

where \(\left(Y^{(k)}_h\right)_{k=1, \ldots, N} \) are i.i.d..

Multilevel Monte Carlo estimator

Taking advantage from the telescopic summation

\[\mathbb{E}\left[Y_{\frac h {n_R}}\right] = \mathbb{E}\left[Y_h\right] + \sum_{j=2}^R \mathbb{E}\left[Y_{\frac h {n_j}} - Y_{\frac h {n_{j-1}}}\right]\]

a Multilevel Monte Carlo estimator writes

\[ I^N_{h,R,M,q} = \frac{1}{N_1} \sum_{k=1}^{N_1} Y_{h}^{(1),k} + \sum_{j=2}^R \frac{1}{N_j} \sum_{k=1}^{N_j} \left(Y_{h_j}^{(j),k} - Y_{h_{j-1}}^{(j),k}\right) \]

where for every \(j\), \((Y^{(j),k})_{k\ge 1}\) is a sequence of independent copies of \(Y^{(j)}\), \(h_j = \frac h {n_j}\), \(n_j = M^{j-1}\) and \(N_j = \lceil N q_j \rceil\) where \(q=(q_1,\dots,q_R) \in \mathcal{S}_+(R) = \left\{q \in (0,1)^R, \sum_{j=1}^R q_j = 1\right\}\).

Multilevel Richardson Romberg estimator

The Multilevel Richardson Romberg estimator is a weighted version of the Multilevel Monte Carlo estimator, which writes

\[ I^N_{h,R,M,q} = \frac{1}{N_1} \sum_{k=1}^{N_1} Y_{h}^{(1),k} + \sum_{j=2}^R \frac{\mathbf{W}_j}{N_j} \sum_{k=1}^{N_j} \left(Y_{h_j}^{(j),k} - Y_{h_{j-1}}^{(j),k}\right) \]

with \(\mathbf{W}_j = \sum_{i=j}^R \mathbf{w}_i\) where \(\mathbf{w} = (\mathbf{w}_1, \ldots, \mathbf{w}_R)\) is the unique solution to the Vandermonde system \(V \mathbf{w} = e_1\) where

\[ V = V(1, n_2^{-\alpha}, \ldots, n_R^{-\alpha}) = \left( \begin{array}{cccc} 1 & 1 & \ldots & 1\\ 1 & n_2^{-\alpha} & \ldots & n_R^{-\alpha} \\ \vdots & \vdots & \ldots & \vdots \\ 1 & n_2^{-\alpha(R-1)} & \ldots & n_R^{-\alpha(R-1)} \end{array}\right) \]

Diffusion : Black-Scholes

We consider a geometric Brownian motion \((X_t)_{t\in [0,T]}\) such that

\[X_t = x_0 e^{\left( r - \frac {\sigma^2} 2 \right) t + \sigma W_t}\]

where \(r\) denotes the riskless interest rate, \(\sigma\) denotes the volatility and \((W_t)_{t \in [0,T]}\) is a standard Brownian motion. \((X_t)_{t \in [0,T]}\) is the solution of the SDE

\[dX_t = X_t(rdt + \sigma dW_t), \; X_0 = x_0\]

If the dimension of the problem is \(d>1\), we consider \((X_t)_{t\in [0,T]} = (X_t^1, \ldots, X_t^d)\) such that

\[\frac {dX_t^i} {X_t^i} = rdt + \sigma^i dW_t^i, \; X_0^i = x_0^i\]

where \((W_t^i)_{t \in [0,T]}\) are \(d\) Brownian motions with correlation matrix \(\rho\), i.e. \(d\langle W^i, W^j \rangle _t = \rho_{ij}dt\).

Let \(f\) be a real-valued \(\beta\) Hölder continuous function, \(\beta \in (0, 1]\). Our aim is to compute

\[\mathbb{E}(Y_0) = \mathbb{E}[f(X_T)]\]

(see for example vanilla options). In the case of path dependent options, we will consider

\[\mathbb{E}(Y_0) = \mathbb{E}[f(X_T, \sup_{t\in [0,T]} X_t)] \quad \text{or} \quad \mathbb{E}(Y_0) = \mathbb{E}[f(X_T, \inf_{t\in [0,T]} X_t)]\]

We consider the continuous Euler scheme \(\bar X^{h} = (\bar X^h_t)_{t \in [0,T]}\) with bias parameter \(h = \frac{T}{n}\) defined by \[ \bar X^h_t = X_0 + \int_0^t b (\bar X^h_{\underline s})\, ds + \int_0^t \sigma (\bar X^h_{\underline s})\, d W_s, \quad \text{where} \quad \underline s = kh \text{ on } [kh, (k+1) h) \]

We make the approximations \(Y_h = f(\bar X^h_T)\) for vanilla options and \(Y_h = f(\bar X^h_T, \sup_{t\in [0,T]} \bar X_t)\) for path dependent options (similarly for inf dependence).

Nested

The aim of a Nested Monte Carlo is to compute \[ \mathbb{E} (Y_0) = \mathbb{E}\left[f\left(\mathbb{E} \left[ X \middle\vert Y \right]\right)\right] \]

where \((X,Y)\) is a couple of \(\mathbb{R} \times \mathbb{R}^q\)-valued random variables defined on a probability space \((\Omega, \mathcal{A}, \mathbb{P})\) with \(X \in \mathbf{L}^2(\mathbb{P})\) and \(f:\mathbb{R} \rightarrow \mathbb{R}\) is a Lipschitz continuous function.

We make the assumption that there exists a Borel function \(F: \mathbb{R}^{q_Z} \times \mathbb{R}^{q_Y} \rightarrow \mathbb{R}\) and a random variable \(Z:(\Omega, \mathcal{A}) \rightarrow \mathbb{R}^{q_Z}\) independent of \(Y\) such that

\[X = F(Z,Y)\]

Then, if \(X\in \mathbf{L}^2\), we can write

\[\mathbb{E} \left( X \middle\vert Y\right)(\omega) = \left( \mathbb{E}[F(Z,y)] \right)_{ | y=Y(\omega)})\]

Then we can set

\[ \mathcal{H} = \left\{ 1/K, K\geq 1 \right\}, \; Y_0 = f\left( \mathbb{E}[X \middle\vert Y] \right), \; Y_{\frac 1 K} = f\left( \frac 1 K \sum_{k=1}^K F(Z_k,Y) \right) \]

The price of a vanilla option with payoff \(\varphi\) is given by \(e^{-rT} \mathbb{E} \left[\varphi \left(S_T\right)\right]\) for a monodimensional problem and by \(e^{-rT} \mathbb{E} \left[\varphi \left(\sum_{i=1}^d \frac 1 d S_T^i\right)\right]\) for dimension \(d>1\).

The price of a path dependent option with functional payoff \(\varphi\) is given by \(e^{-rT} \mathbb{E} \left[\varphi \left(\left(S_t\right)_{t\in [0,T]}\right)\right]\) for a monodimensional problem and by \(e^{-rT} \mathbb{E} \left[\varphi \left(\sum_{i=1}^d \frac 1 d \left(S_t^i\right)_{t\in [0,T]}\right)\right]\) for dimension \(d>1\).

Here a list of payoff :

Vanilla

Call option on stock price \((S_t)_{t\in [0,T]}\), with strike \(K\) :

\[\varphi(S_T) = (S_T - K)^+\]

Put option on stock price \((S_t)_{t\in [0,T]}\), with strike \(K\) :

\[\varphi(S_T) = (K-S_T)^+\]

Barrier Call

Barrier up-and-out call option on \((S_t)_{t\in [0,T]}\) with strike \(K\), barrier \(B>K\) and current price \(S_0<B\) :

\[\varphi(x) = \left( x(T) - K \right)^+ \mathbb{1} _{\left\{ \max_{t\in [0,T]} x(t)\leq B\right\}}\]

Barrier up-and-in call option on \((S_t)_{t\in [0,T]}\) with strike \(K\), barrier \(B\) and current price \(S_0<B\) :

\[\varphi(x) = \left( x(T) - K \right)^+ \mathbb{1} _{\left\{ \max_{t\in [0,T]} x(t)\geq B\right\}}\]

Barrier down-and-out call option on \((S_t)_{t\in [0,T]}\) with strike \(K\), barrier \(B\) and current price \(S_0>B\) :

\[\varphi(x) = \left( x(T) - K \right)^+ \mathbb{1} _{\left\{ \min_{t\in [0,T]} x(t)\geq B\right\}}\]

Barrier down-and-in call option on \((S_t)_{t\in [0,T]}\) with strike \(K\), barrier \(B>K\) and current price \(S_0>B\) :

\[\varphi(x) = \left( x(T) - K \right)^+ \mathbb{1} _{\left\{ \min_{t\in [0,T]} x(t)\leq B\right\}}\]

Barrier Put

Barrier up-and-out put option on \((S_t)_{t\in [0,T]}\) with strike \(K\), barrier \(B\) and current price \(S_0<B\) :

\[\varphi(x) = \left( K - x(T) \right)^+ \mathbb{1} _{\left\{ \max_{t\in [0,T]} x(t)\leq B\right\}}\]

Barrier up-and-in put option on \((S_t)_{t\in [0,T]}\) with strike \(K\), barrier \(B<K\) and current price \(S_0<B\) :

\[\varphi(x) = \left( K -x(T) \right)^+ \mathbb{1} _{\left\{ \max_{t\in [0,T]} x(t)\geq B \right\}}\]

Barrier down-and-out put option on \((S_t)_{t\in [0,T]}\) with strike \(K\), barrier \(B<K\) and current price \(S_0>B\) :

\[\varphi(x) = \left( K -x(T) \right)^+ \mathbb{1} _{\left\{ \min_{t\in [0,T]} x(t)\geq B\right\}}\]

Barrier down-and-in put option on \((S_t)_{t\in [0,T]}\) with strike \(K\), barrier \(B\) and current price \(S_0>B\) :

\[\varphi(x) = \left( K -x(T) \right)^+ \mathbb{1} _{\left\{ \min_{t\in [0,T]} x(t)\leq B\right\}}\]

Lookback

Lookback Call with floating strike option :

\[\varphi(x) = \left( x(T) - \lambda \min_{t \in [0,T]} x(t)\right)^+ \]

Lookback Put with floating strike option :

\[\varphi(x) = \left( \lambda \max_{t \in [0,T]} x(t) - x(T) \right)^+ \]

Nested Put-on-Call

A Put-on-Call option with expiration dates \(T_1 < T_2\) and strikes \(K_1\) and \(K_2\) on a stock price \(S_t\) is an option where the holder has at time \(T_1\) the right to sell, using the strike \(K_1\), a new Call with strike \(K_2\) and maturity \(T_2\). The payoff of such an option writes

\[\left( K_1 -\mathbb{E} \left[\left(S_{T_2} - K_2 \right) ^+ \middle\vert S_{T_1} \right] \right)^+\]