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)^+\]

Set here below the constants of the problem as defined in the option and the model description.

\(\varepsilon\) is the precision of the problem. For these products we set \(0<\varepsilon < 1\). We strongly recommend to choose \(\varepsilon<0.1\). Nevertheless, very low values of \(\varepsilon\) can lead to very high computational costs and consequently to long execution times. The heaviest computation is the one of the Crude Monte Carlo estimator, make sure you do not have a too big computational cost for this estimator before launching the computation.

The seed used to initialize the pseudorandom number generator changes at each call of the algorithm if we set the variable \(seed = 0\). If we choose an integer different from zero, this seed will be taken to initialize the pseudorandom generator. This can be useful when we want to reproduce twice the same simulation.