A framework for adaptive Monte-Carlo procedures

The sequence $(\theta_{n})_{n}$ can be any sequence adapted to $(X_{n})_{n\geq 1}$ convergent or not. For instance, $(\theta_{n})_{n}$ can be an ergodic Markov chain distributed around the minimizer $\theta^{\star}$ such as Monte Carlo Markov Chain algorithms.

To derive a central limit theorem for the adaptive estimator $\xi_{n}$, we need a central limit theorem for locally square integrable martingales, whose convergence rate has been extensively studied. We refer to the works of Rebolledo [21], Jacod and Shiryaev [7], Hall and Heyde [6] and Whitt [25] to find different statements of central limit theorems for locally square integrable càdlàg martingales in continuous time, from which theorems can easily be deduced for discrete time locally square integrable martingales.

The assumptions of Theorem 2.3 are fairly easy to check in practice since they are formulated independently of the sequence $(\theta_{n})_{n}$. When $\theta_{\infty}=\theta^{\star}$, which is nonetheless not required, the limiting variance is optimal in the sense that a crude Monte Carlo computation with the optimal parameter $\theta^{\star}$ would have lead to the same limiting variance. These assumptions are satisfied in the frameworks introduced in Section 4.

Let ${(X_{n})}_{n\geq 1}$ be an i.i.d sequence of random variables following the law of $X$ and ${(\gamma_{n})}_{n\geq 1}$ be a decreasing sequence of a positive real numbers satisfying

$$\sum_{n}\gamma_{n}=\infty\quad\mbox{and}\quad\sum_{n}\gamma_{n}^{2}<\infty.$$

(5)

The sequence $(\gamma_{n})_{n}$ is often called the gain sequence or the step sequence. We define the $\sigma-$field ${\mathcal{F}}_{n}=\sigma(X_{k},\;k\leq n)$. We introduce an increasing sequence of compact sets $(K_{j})_{j}$ of ${\mathbb{R}}^{d}$

$$\bigcup_{n=0}^{\infty}K_{n}\>=\>{\mathbb{R}}^{d}\quad\mbox{and}\quad K_{n}\subsetneq\mathring{K}_{n+1}$$

(6)

Now, we can present the randomly truncated stochastic algorithm introduced in [3], which essentially consists in a truncation of the Robbins Monro algorithm on an increasing sequence of compact sets. For $\theta_{0}\in K_{0}$ and $\alpha_{0}=0$, we define the sequences of random variables ${(\theta_{n})}_{n}$ and ${(\alpha_{n})}_{n}$ by

$$\begin{cases}&\theta_{n+\frac{1}{2}}=\theta_{n}-\gamma_{n+1}U(\theta_{n},X_{n+1}),\\ \text{if $\theta_{n+\frac{1}{2}}\in\mathcal{K}_{\alpha_{n}}$}&\theta_{n+1}=\theta_{n+\frac{1}{2}}\quad\mbox{ and }\quad\alpha_{n+1}=\alpha_{n},\\ \text{if $\theta_{n+\frac{1}{2}}\notin\mathcal{K}_{\alpha_{n}}$}&\theta_{n+1}=\theta_{0}\quad\mbox{ and }\quad\alpha_{n+1}=\alpha_{n}+1.\end{cases}$$

(7)

$\theta_{n+\frac{1}{2}}$ is the new sample we draw, either we accept it and set $\theta_{n+1}=\theta_{n+\frac{1}{2}}$ or we reject it and reset the algorithm to $\theta_{0}$ when it tries to jump too far ahead in a single step. Note that $\theta_{n+\frac{1}{2}}$ is actually drawn along the dynamics of the Robbins Monro algorithm and either we accept it as the new iterate or we reject it when the algorithm tries to jump to far ahead and in this case we reset the new iterate to $\theta_{0}$. In the following, we write Equation (7) in a more condensed form

$$\theta_{n+1}={\cal T}_{K_{\alpha_{n}}}\left(\theta_{n}-\gamma_{n+1}U(\theta_{n},X_{n+1})\right),$$

(8)

where ${\cal T}_{K_{\alpha_{n}}}$ denotes the truncation on the compact sets $K_{\alpha_{n}}$.

The use of truncations enables to relax the hypotheses required to ensure the convergence. From the recent results of [16], we can state the following convergence result

Note that the assumptions required to ensure the convergence are very weak and are formulated independently of the algorithm trajectories, which makes them easy to check. Since the variance reduction technique we settle here aims at being automatic in the sense that it does not require any fiddling with the gain sequence depending on the function $U$, it is quite natural to average the procedure defined by Equation (7).

This section is based on the remark that Cesaro type averages tend to smooth the behaviour of convergent estimators at least from a theoretical point of view. Such averaging techniques have already been studied and proved to provide asymptotically efficient estimators (see for instance [20], [14] or [19]).

At the same time, it is well known that true Cesaro averages are not so efficient from a practical point of view because the rate at which the impact of the first iterates vanishes in the average is too slow and it induces some kind of a numerical bias which in turn dramatically slows down the convergence. Combining these two facts has led us to consider a moving window average of Algorithm (7).
In this section, we restrict to gain sequences of the form $\gamma_{n}=\frac{\gamma}{(n+1)^{a}}$ with $\frac{1}{2}<a<1$. Let $\tau>0$ be the length of the window used for averaging. For $n\geq 1$, we introduce

$$\hat{\theta}_{n}(\tau)=\frac{\gamma_{p}}{\tau}\sum_{i=p}^{p+\lfloor\tau/\gamma_{p}\rfloor}\theta_{i}\quad\text{ with }p=\sup\{k\geq 1:k+\tau/\gamma_{k}\leq n\}\wedge n.$$

(9)

We use the convention $\sup\emptyset=+\infty$. The almost sure convergence of $(\hat{\theta}_{n}(\tau))_{n}$ can easily be deduced from Theorem 3.1. The asymptotic normality of the sequence $(\hat{\theta}_{n}(\tau))_{n}$ has been studied in [15]. The definition of $\hat{\theta}_{n}$ is a little different from the one used in [15] because we want to ensure that the sequence $(\hat{\theta}_{n})_{n}$ is adapted to the filtration ${\mathcal{F}}_{n}$ in view of the use of $\hat{\theta}_{n}$ as an estimator of $\theta^{\star}$ in Algorithm 1.2.

Let $G=(G_{1},\ldots,G_{d})$ be a $d-$dimensional standard normal random vector. For any measurable function $h:{\mathbb{R}}^{d}\longrightarrow{\mathbb{R}}$ such that ${\mathbb{E}}(|h(G)|)<\infty$, one has for all $\theta\in{\mathbb{R}}^{d}$

$${\mathbb{E}}\left(h(G)\right)={\mathbb{E}}\left(e^{-\theta\cdot G-\frac{|\theta|^{2}}{2}}h(G+\theta)\right).$$

(10)

Assume we want to compute ${\mathbb{E}}(f(G))$ for a measurable function $f:{\mathbb{R}}^{d}\longrightarrow{\mathbb{R}}$ such that $f(G)$ is integrable. By applying equality (10) to $h=f$ and $h(x)=f^{2}(x)\mathop{\mathrm{e}^{-\theta\cdot x+\frac{|\theta|^{2}}{2}}}$, one obtains that the expectation and the variance of the random variable $f(G+\theta)\mathop{\mathrm{e}^{-\theta\cdot G\nobreakspace-\frac{|\theta|^{2}}{2}}}$ are respectively equal to ${\mathbb{E}}(f(G))$ and $v(\theta)-{\mathbb{E}}^{2}(f(G))$ where

$$v(\theta)={\mathbb{E}}\left(f^{2}(G)\mathop{\mathrm{e}^{-\theta\cdot G+\frac{|\theta|^{2}}{2}}}\right).$$

The strict convexity of the function $v$ is already known from [23] for instance. For the sake of completeness, we prove here a slightly improved version of this result.

Proposition 4.1 implies that $v$ has a unique minimiser $\theta^{\star}$ characterised by $\nabla v(\theta^{\star})=0$, i.e. ${\mathbb{E}}\left((\theta^{\star}-G)e^{-\theta^{\star}\cdot G+\frac{|\theta^{\star}|^{2}}{2}}f^{2}(G)\right)=0$.

Equality (10) can actually be extended to the Brownian motion framework using Girsanov’s theorem. Let $(W_{t},0\leq t\leq T)$ be a $d-$dimensional Brownian motion and ${\mathcal{F}}$ its natural filtration. For any measurable and ${\mathcal{F}}-$predictable process $(\theta_{t},0\leq t\leq T)$ such that ${\mathbb{E}}\left(\mathop{\mathrm{e}^{\frac{1}{2}\int_{0}^{T}|\theta_{t}|^{2}dt}}\right)<\infty$, one has

$${\mathbb{E}}\left(f(W_{t},0\leq t\leq T)\right)={\mathbb{E}}\left(e^{-\int_{0}^{T}\theta_{t}\cdot dW_{t}-\frac{1}{2}\int_{0}^{T}|\theta_{t}|^{2}dt}f\left(W_{t}+\int_{0}^{t}\theta_{s}ds,0\leq t\leq T\right)\right).$$

Assume $\theta_{t}=\theta\in{\mathbb{R}}^{d}$, for all $t\in[0,T]$. The variance of $\mathop{\mathrm{e}^{-\theta\cdot W_{T}-\frac{\theta^{2}T}{2}}}f(W_{t},0\leq t\leq T)$ writes down $v(\theta)-{\mathbb{E}}\left(f^{2}(W_{t},0\leq t\leq T)\right)$ with

$$v(\theta)={\mathbb{E}}\left(e^{-\theta\cdot W_{T}+\frac{|\theta|^{2}}{2}T}f^{2}\left(W_{t}+\theta t,\,0\leq t\leq T\right)\right).$$

A similar result to Proposition 4.1 holds; in particular $v$ is infinitely continuously differentiable, strictly convex and goes to infinity at infinity.

For more general processes $(\theta_{t},0\leq t\leq T)$, we refer the reader to [17].

The idea of tilting some probability measure to find the ones that minimises the variance is a very common idea which can be also be applied to a wide range of distribution, see for instance the recent results of Kawai [11, 10] in which he applied an exponential change of measure to Lévy processes, also known as the Esscher transform.

Consider a random variable $X$ with values in ${\mathbb{R}}^{d}$ and cumulative generating function $\psi(\theta)=\log{\mathbb{E}}\big{(}\mathop{\mathrm{e}^{\theta\cdot X}}\big{)}$. We assume that $\psi(\theta)<\infty$ for all $\theta\in{\mathbb{R}}^{d}$. Let $p$ denote the density of $X$. We define the density $p_{\theta}$ by

$$p_{\theta}(x)=p(x)\mathop{\mathrm{e}^{\theta\cdot x-\psi(\theta)}},\quad x\in{\mathbb{R}}^{d}.$$

Let $X^{(\theta)}$ have $p_{\theta}$ as a density, then

$${\mathbb{E}}(f(X))={\mathbb{E}}\left[f(X^{(\theta)})\frac{p(X^{(\theta)})}{p_{\theta}(X^{(\theta)})}\right].$$

The variance of $f(X^{(\theta)})\frac{p(X^{(\theta)})}{p_{\theta}(X^{(\theta)})}$ writes $v(\theta)-{\mathbb{E}}\Big{(}f(X^{(\theta)})^{2}\frac{p(X^{(\theta)})^{2}}{p_{\theta}(X^{(\theta)})^{2}}\Big{)}$ with

$$v(\theta)={\mathbb{E}}\left(\left|f(X^{(\theta)})\frac{p(X^{(\theta)})}{p_{\theta}(X^{(\theta)})}\right|^{2}\right)={\mathbb{E}}\left(f(X)^{2}\mathop{\mathrm{e}^{-\theta\cdot X+\psi(\theta)}}\right).$$

Obviously, this change of measure is only valuable as a variance reduction technique if $X^{(\theta)}$ can be simulated at approximately the same cost as $X$.

We consider a $D-$multidimensional local volatility model in which each asset is supposed to be driven by the following dynamics under the risk neutral measure.

$$dS^{i}_{t}=S^{i}_{t}(rdt+\sigma(t,S^{i}_{t})\cdot dW^{i}_{t}),\ S^{i}_{0}=s^{i}.$$

$W=(W^{1},\dots,W^{D})^{*}$ is a vector of correlated standard Brownian motions. The covariance structure of the Brownian motions is given by $\langle W,W\rangle_{t}=\Gamma t$ where $\Gamma$ is a definite positive matrix with a diagonal filled with ones. In our numerical examples, we take $\Gamma_{ij}={\bf 1}_{\{i=j\}}+\rho{\bf 1}_{\{i\neq j\}}$ with $\rho\in(\frac{-1}{D-1},1)$ to ensure that the matrix $\Gamma$ is positive definite. The function $\sigma$ is the local volatility function, $r$ is the instantaneous interest rate and the vector $(s^{1},\dots,s^{D})$ is the vector of the spot values. In this model, we want to price path-dependent options whose payoffs can be written as a function of $(S_{t},t\leq T)$. Hence, the price is given by the discounted expectation $\mathop{\mathrm{e}^{-rT}}{\mathbb{E}}(\psi(S_{t},t\leq T))$. Most of the time, this expectation must be computed by Monte Carlo methods and one has to consider an approximation of $\psi(S_{t},t\leq T)$ on a time grid $0=t_{0}<t_{1}<\dots<t_{N}=T$. Then, the quantity of interest becomes

$$\mathop{\mathrm{e}^{-rT}}{\mathbb{E}}(\hat{\psi}(S_{t_{0}},S_{t_{1}},\ldots,S_{t_{N}})).$$

The discretisation of the asset $S$ can for instance be obtained using an Euler scheme, which means that the function $\hat{\psi}$ can be expressed in terms of the Brownian increments or equivalently using a random normal vector. These remarks finally turn the original pricing problem into the computation of an expectation of the form ${\mathbb{E}}(\phi(G))$ where $G$ is a standard normal random vector in ${\mathbb{R}}^{ND}$ and $\phi:{\mathbb{R}}^{N\times D}\longmapsto{\mathbb{R}}$ is a measurable and integrable function. Using Equation (10), we have for all $\theta\in{\mathbb{R}}^{d}$,

$${\mathbb{E}}(\phi(G))={\mathbb{E}}\left(\phi(G+A\theta)e^{-A\theta\cdot G-\frac{|A\theta|^{2}}{2}}\right),$$

(15)

where $A$ is $d\times ND$ matrix. The particular choice $d=ND$ and $A=I_{d}$ corresponds to Equation (10). When $D=1$, the choice $d=1$ and $A=(\sqrt{t_{1}},\sqrt{t_{2}-t_{1}},\dots,\sqrt{t_{N}-t_{N-1}})^{*}$ corresponds to adding a linear drift to the one dimensional standard Brownian motion $W$ and we recover the Cameron-Martin formula.

Transformation (15) actually relies on an importance sampling change of measure. Other strategies may be applicable such as stratification for instance as it is explained by Glasserman et al. in [5].

It ensues from Proposition 4.1, that the second moment

$$v(\theta)={\mathbb{E}}\left(\psi(G+A\theta)^{2}e^{-2A\theta\cdot G-|A\theta|^{2}}\right)={\mathbb{E}}\left(\phi(G)^{2}e^{-A\theta\cdot G+\frac{|A\theta|^{2}}{2}}\right)$$

(16)

is strongly convex, infinitely differentiable and

$$\nabla v(\theta)={\mathbb{E}}\left(A^{*}(A\theta-G)\phi(G)^{2}e^{-A\theta\cdot G+\frac{|A\theta|^{2}}{2}}\right).$$

(17)

If we apply Equation (10) again, we obtain an other expression for

$$\nabla v(\theta)={\mathbb{E}}\left(-A^{*}G\phi(G+A\theta)^{2}e^{-2A\theta\cdot G+|A\theta|^{2}}\right).$$

(18)

Let us introduce the following two functions

	$\displaystyle U^{1}(\theta,G)$	$\displaystyle=A^{*}(A\theta-G)\phi(G)^{2}e^{-A\theta\cdot G+\frac{|A\theta|^{2}}{2}},$		(19)
	$\displaystyle U^{2}(\theta,G)$	$\displaystyle=-A^{*}G\phi(G+A\theta)^{2}e^{-2A\theta\cdot G+|A\theta|^{2}}.$		(20)

Using either Equation (17) or (18), we can write $\nabla v(\theta)={\mathbb{E}}(U^{2}(\theta,G))={\mathbb{E}}(U^{1}(\theta,G))$ and these two functions $U^{1}$ and $U^{2}$ fit in the framework of Section 3 and enable to construct two estimators of $\theta^{\star}$ $(\theta^{1}_{n})_{n}$ and $(\theta^{2}_{n})_{n}$ following Equation (7)

	$\displaystyle\theta^{1}_{n+1}$	$\displaystyle={\cal T}_{K_{\alpha_{n}}}\left(\theta^{1}_{n}-\gamma_{n+1}U^{1}(\theta^{1}_{n},G_{n+1})\right),$		(21)
	$\displaystyle\theta^{2}_{n+1}$	$\displaystyle={\cal T}_{K_{\alpha_{n}}}\left(\theta^{1}_{n}-\gamma_{n+1}U^{2}(\theta^{2}_{n},G_{n+1})\right),$		(22)

where $G_{n}$ is an i.i.d sequence of random variables following the law of $G$. We also introduce their corresponding averaging versions $(\widehat{\theta^{1}}_{n})_{n}$ and $(\widehat{\theta^{2}}_{n})_{n}$ following Equation (9). Based on Equation (15), we define

$$H(\theta,G)=\phi(G+A\theta)e^{-A\theta\cdot G-\frac{|A\theta|^{2}}{2}}.$$

Corresponding to the different estimators of $\theta^{\star}$ listed above, we can define as many approximations of ${\mathbb{E}}(\phi(G))$ following Equation (3)

	$\displaystyle\xi^{1}_{n}=\frac{1}{n}\sum_{i=1}^{n}H(\theta^{1}_{i-1},G_{i}),\quad$	$\displaystyle\xi^{2}_{n}=\frac{1}{n}\sum_{i=1}^{n}H(\theta^{2}_{i-1},G_{i}),$
	$\displaystyle\hat{\xi}^{1}_{n}=\frac{1}{n}\sum_{i=1}^{n}H(\hat{\theta}^{1}_{i-1},G_{i}),\quad$	$\displaystyle\hat{\xi}^{2}_{n}=\frac{1}{n}\sum_{i=1}^{n}H(\hat{\theta}^{2}_{i-1},G_{i}),$

where the sequence $G_{i}$ has already been used to build the $(\theta_{n})_{n}$ estimators. From Proposition 4.1 and Theorems 3.1, 2.1 and 2.3, we can deduce the following result.