Recursive Marginal Quantization of Higher-Order Schemes

TL;DR

This paper extends recursive marginal quantization to higher-order schemes like Milstein, improving accuracy in stochastic process simulations and option pricing, with efficient implementation and boundary handling.

Contribution

It introduces recursive marginal quantization for Milstein and weak order 2.0 schemes, enhancing numerical accuracy and efficiency in stochastic differential equation applications.

Findings

Higher order schemes show improved weak order convergence.

Numerical results demonstrate better approximation of analytical distributions.

Enhanced option pricing accuracy for European, Bermudan, and Barrier options.

Abstract

Quantization techniques have been applied in many challenging finance applications, including pricing claims with path dependence and early exercise features, stochastic optimal control, filtering problems and efficient calibration of large derivative books. Recursive Marginal Quantization of the Euler scheme has recently been proposed as an efficient numerical method for evaluating functionals of solutions of stochastic differential equations. This method involves recursively quantizing the conditional marginals of the discrete-time Euler approximation of the underlying process. By generalizing this approach, we show that it is possible to perform recursive marginal quantization for two higher-order schemes: the Milstein scheme and a simplified weak order 2.0 scheme. As part of this generalization a simple matrix formulation is presented, allowing efficient implementation. We further extend the applicability of recursive marginal quantization by showing how absorption and reflection at the zero boundary may be incorporated, when this is necessary. To illustrate the improved accuracy of the higher order schemes, various computations are performed using geometric Brownian motion and its generalization, the constant elasticity of variance model. For both processes, we show numerical evidence of improved weak order convergence and we compare the marginal distributions implied by the three schemes to the known analytical distributions. By pricing European, Bermudan and Barrier options, further evidence of improved accuracy of the higher order schemes is demonstrated.

Equations

$$\mathbb{E}\left[H(X)\right]=\int_{\mathbb{R}}H(x)\,d\mathbb{P}(X\leq x)\approx\sum_{\gamma\in\Gamma}H(\gamma)\mathbb{P}\big{(}\widehat{X}=\gamma\big{)}.$$

$$\displaystyle\pi_{\Gamma}(X)=\big{\{}\gamma^{i}\in\Gamma\,\big{|}\,\,\|X-\gamma^{i}\|\leq\|X-\gamma^{j}\|\text{ for all $j=1,\ldots,N;$ }$$

$$\displaystyle\text{holds only for $i<j$}\big{\}}.$$

$$R_{i}(\Gamma)=\big{\{}x\in\mathbb{R}\,\big{|}\,\,\pi_{\Gamma}(x)=\gamma^{i}\big{\}}.$$

$$R_i\cap R_j=\varnothing \text{ for }i\ne j\text{and}{\cup }_{i=1}^N{R}_i=\mathbb{R}.$$

$$r^{i-}=\frac{{\gamma }^{i-1}+{\gamma }^i}{2}\text{and}{r}^{i+}=\frac{{\gamma }^i+{\gamma }^{i+1}}{2},$$

$$D(\Gamma )$$

$$={\int }_{\mathbb{R}}\Vert x-{\pi }_{\Gamma }(x){\Vert }^{2}d\mathbb{P}(X\le x)$$

$$=\sum _{i=1}^N{\int }_{R_i(\Gamma )}\Vert x-{\gamma }^i{\Vert }^{2}d\mathbb{P}(X\le x).$$

$$\frac{\partial D(\Gamma )}{\partial {\gamma }^i}=2\left({\int }_{R_i(\Gamma )}({\gamma }^i-x)d\mathbb{P}(X\le x)\right),$$

$${\gamma }_{n+1}^i=\frac{{\int }_{R_i({\Gamma }^{(n)})}xd\mathbb{P}(X\le x)}{{\int }_{R_i({\Gamma }^{(n)})}d\mathbb{P}(X\le x)},$$

$${\mathbf{\Gamma }}^{(n+1)}={\mathbf{\Gamma }}^{(n)}-{\left[{\nabla }^{2}D\left({\mathbf{\Gamma }}^{(n)}\right)\right]}^{-1}\nabla D\left({\mathbf{\Gamma }}^{(n)}\right)$$

$$M_X^p(x)=\mathbb{E}\left[X^p{\mathbb{I}}_{\left\{X<x\right\}}\right],$$

$$D(\Gamma )$$

$$\displaystyle=\sum_{i=1}^{N}\Big{[}M_{X}^{2}(r^{i+})-M_{X}^{2}(r^{i-})-2\gamma^{i}\left(M_{X}^{1}(r^{i+})-M_{X}^{1}(r^{i-})\right)$$

$$\displaystyle\qquad\quad\quad+(\gamma^{i})^{2}\left(F_{X}(r^{i+})-F_{X}(r^{i-})\right)\Big{]}.$$

$$\frac{\partial D(\Gamma )}{\partial {\gamma }^i}=2{\gamma }^i\left(F_X(r^{i+})-F_X(r^{i-})\right)-2\left(M_X^1(r^{i+})-M_X^1(r^{i-})\right),$$

$$\frac{{\partial }^{2}D(\Gamma )}{\partial ({\gamma }^i{)}^2}$$

$$\frac{{\partial }^{2}D(\Gamma )}{\partial {\gamma }^i\partial {\gamma }^{i+1}}$$

$$\frac{{\partial }^{2}D(\Gamma )}{\partial {\gamma }^i\partial {\gamma }^{i-1}}$$

$$[\mathbf{\Gamma }{]}_i$$

$$[\mathbf{f}{]}_i$$

$$[\mathbf{p}{]}_i$$

$$=\mathbb{P}(X\in R_i(\Gamma ))$$

$$=F_X(r^{i+})-F_X(r^{i-}),$$

$$\mathbb{E}\left[H(X)\right]=\sum_{i=1}^{N}H(\gamma^{i})\mathbb{P}\big{(}\widehat{X}=\gamma^{i}\big{)}=\mathbf{p}H(\mathbf{\Gamma}),$$

$$\nabla D(\mathbf{\Gamma })=2\mathbf{\Gamma }\circ {\mathbf{p}}^{\top }-2\mathbf{M},$$

$${\mathbf{h}}_{\mathsf{off}}=-\frac{1}{2}[\mathbf{f}\circ \Delta \mathbf{\Gamma }{]}^{\top },$$

$${\mathbf{h}}_{\mathsf{main}}=2\mathbf{p}+[{\mathbf{h}}_{\mathsf{off}}|0]+[0|{\mathbf{h}}_{\mathsf{off}}],$$

$$f_X(x)$$

$$F_X(x)$$

$$M_X^1(x)$$

$${\gamma }^n=\frac{5.5n}{N+1}-2.75,$$

$$f_X(x)$$

$$F_X(x)$$

$$M_X^1(x)$$

$$x^{\pm }$$

$${\gamma }^n=\{\begin{array}{ll}{\left(\frac{(3+\sqrt{\lambda })n}{N}\right)}^2& \text{for λ<2.5}\\ {\left(\frac{5n}{N+1}-2.5+\sqrt{\lambda }\right)}^2& \text{for λ≥2.5,}\end{array}$$

Full text

1 Introduction

Quantization techniques have been shown to be effective in a number of quantitative finance applications. In particular, this includes the pricing of contingent claims with path dependency and early exercise [Pàges and Wilbertz, 2009; Sagna, 2011; Bormetti et al., 2016], stochastic control problems [Pàges et al., 2004] and nonlinear filtering [Pàges and Pham, 2005]. To improve numerical efficiency, Pagès and Sagna [2015] introduced a technique known as recursive marginal quantization. This technique has been shown to be effective for fast calibration of large derivative books by Callegaro et al. [2014, 2015a], another challenging application in finance. This paper aims to significantly improve the accuracy of recursive marginal quantization.

Quantization is a lossy compression technique that produces a discrete representation of a signal using less information than the original. The technique originated in the field of signal compression, but has found application in fields as far-reaching as signal processing, pattern recognition, data mining, integration theory and, more recently, numerical probability. For general overviews of the mathematics and applications of quantization see Du et al. [1999] and Pagès [2014].

Vector quantization of probability distributions was formalized in the work of Graf and Luschgy [2000] and has been applied to the field of mathematical finance since its inception. It is a technique for optimally representing a continuous distribution by a discrete distribution, where a measure of the ‘distance’ between the two, called the distortion, is minimized. The distortion is most commonly measured using the squared Euclidean error.

The application of vector quantization to the solution of finance-related problems generally proceeds by discretising time and then quantizing the corresponding marginal distributions of the system of stochastic differential equations specific to the problem. The quantized grids and their associated weights are then used to compute the expectations required in pricing contingent claims (including claims with early exercise) or for performing the optimizations required in stochastic control problems [Pagès et al., 2004]. Due to the reliance on Lloyd’s Algorithm [Lloyd, 1982], or variants thereof, these approaches generally incur a heavy computational burden.

A more efficient approach for the single-factor case has recently been proposed by Pagès and Sagna [2015]. Known as recursive marginal quantization (RMQ), it makes use of a Newton-Raphson iteration to quantize the Euler-Maruyama [Maruyama, 1955] updates of the underlying SDE. This technique has been used to provide fast calibration of a local volatility model by Callegaro et al. [2014, 2015a], and extended for use with two factor SDEs and applied to stochastic volatility models [Callegaro et al., 2015b].

In the present work, the RMQ algorithm is generalized, allowing the implementation of schemes of higher order than the Euler-Maruyama scheme. We now provide an overview of the rest of the paper.

Section Two provides a review of vector quantization (VQ) as applied to probability distributions. In particular, we strive to simultaneously provide a precise, concise and intuitive description of the methodology. The resulting algorithm is presented using a matrix formulation, allowing for efficient implementation. The section concludes by showing examples of vector quantization applied to the Gaussian and noncentral chi-squared distributions.

The third section introduces recursive marginal quantization applied to stochastic differential equations. Our formulation of the problem is presented in more generality than the original formulation by Pagès and Sagna [2015]. A matrix formulation, which demonstrates the connection with Markov chains, is provided, allowing easy and efficient implementation.

In the section that follows we extend the RMQ algorithm to higher-order updates, specifically the Milstein [1975] scheme and a simplified weak order 2.0 Taylor scheme of Kloeden and Platen [1999]. Geometric Brownian motion (GBM) and the constant elasticity of variance (CEV) process serve as examples to illustrate the improved error in the quantized marginal distributions and the improved weak order convergence.

When performing a Monte Carlo simulation of a discrete-time approximation of a process, non-negativity of the solution is usually enforced by implementing absorption or reflection. Under certain circumstances this is also required when using RMQ. In Section Five we present the modifications of the RMQ algorithm necessary to ensure an absorbing or reflecting boundary at zero. These modifications allow the RMQ algorithm to be applied to the CEV process for parameter sets that would otherwise be problematic under the original formulation.

Section Six presents numerical results of the application of the three RMQ schemes to option pricing. European, barrier and Bermudan options are priced under the GBM and CEV models. Where possible, the results are compared to available closed-form solutions, otherwise they are compared to high-resolution finite difference or Monte Carlo implementations.

One of the goals of this paper has been to produce a more general and easily accessible introduction to the theory of VQ and RMQ. The primary contribution is, however, the more general formulation of RMQ that enables the systematic extension of the work of Pagès and Sagna [2015]. The paper concludes with a discussion of ongoing work.

2 Vector Quantization

Vector quantization is a lossy compression technique that provides a way to encode a vector space using a discrete subspace. While the technique is applicable more generally, we shall only consider the quantization of one dimensional distributions. The vector quantization problem we aim to address in this section may be specified intuitively as follows:

Find the discrete distribution that “best” represents the continuous distribution function associated with a random variable $X$.

This is depicted in the Figure 1, which shows a density function of a continuous random variable and its corresponding quantized version. Here, for ease of visualization, we have chosen to plot the probability density function (of the continuous random variable) and the probability mass function of the quantizer instead of the continuous and discrete distribution functions.

We now provide a more rigorous specification of the problem. Let $X$ be a continuous random variable, taking values in $\mathbb{R}$, and defined on a probability space $(\Omega,\mathscr{F},\mathbb{P})$. The above question may be reposed as:

How does one optimally approximate $X$, in a least-squares sense, by a discrete random variable $\widehat{X}:\Omega\rightarrow\Gamma$, where $\Gamma$ is a finite set of elements in $\mathbb{R}$?

The reason that quantization is useful is that it allows efficient approximations of expectations of functionals $H(X)$ of $X$, that is,

[TABLE]

Here, for example, $H$ may be the discounted payoff of a financial claim and $\mathbb{P}$ the risk-neutral probability measure.

Consider the approximation of $X$ given by $\hat{X}$, a discrete random vector defined as the nearest-neighbour projection of $X$ onto $\Gamma=\{\gamma^{1},\gamma^{2},\ldots,\gamma^{N}\}$, a set of distinct points in $\mathbb{R}$ with finite cardinality $N\in\mathbb{N^{+}}$. We shall refer to $\Gamma$ as the quantizer and its elements as codewords. The nearest neighbor projection operator $\pi_{\Gamma}:\mathbb{R}\rightarrow\Gamma$ is defined as

[TABLE]

Here, $\|\cdot\|$ is the Euclidean norm. Associated with the quantizer, the regions $R_{i}(\Gamma)\subset\mathbb{R}$ are the subset of values of $X$ that are mapped to each codeword $\gamma^{i}$:

[TABLE]

These are also known as Voronoi regions. For the sake of brevity we shall use $R_{i}$ to refer to $R_{i}(\Gamma)$ when it is clear that the quantizer we are referring to is $\Gamma$. The set of regions $\{R_{i}\}_{i=1}^{N}$ is called a tessellation of $\mathbb{R}$, and has the following properties:

[TABLE]

Since we are working in one dimension, the regions $R_{i}$ may be defined directly as $R_{i}=\{x\,|\,r^{i-}<x\leq r^{i+}\}$ with

[TABLE]

for $1\leq i\leq N$, where, by definition, $r^{1-}=-\infty$ and $r^{N+}=\infty$. If the distribution under consideration is not defined over the whole real line, then $r^{1-}$ and $r^{N+}$ are adjusted to reflect the interval of support. Figure 2 shows a simple graphical representation of these regions.

With these definitions in place, we are now in a position to precisely define the optimization problem. We wish to find the quantizer $\Gamma$ such that $\widehat{X}=\pi_{\Gamma}(X)$ best approximates $X$. The sense in which the quantizer $\Gamma$ is “best” is specified by the distortion function

[TABLE]

We require the $\Gamma$ that minimizes $D(\Gamma)$. The probability weights then follow directly as a result of the nearest neighbor projection operator, i.e., $\mathbb{P}(\widehat{X}=\gamma^{i})=\mathbb{P}(X\in R_{i})$. A necessary condition on the optimal $\Gamma$ is that the gradient of the distortion function is zero, that is, $\nabla D(\Gamma)=\bar{0}$, where the elements of $\nabla D(\Gamma)$ are given by

[TABLE]

for $1\leq i\leq N$. Intuitively, this means that the first moment conditioned on the outcomes of a region equals the respective codeword. Thus, one way to solve this system of equations is to set up a fixed-point iteration using the above expression for the gradient. This is the basis for Lloyd’s algorithm which starts with an initial guess for the quantizer, $\Gamma^{(0)}$, and generates successive updates, $\Gamma^{(n+1)}$, with the new codewords, $\gamma_{n+1}^{i}$, computed as the centroids of the regions associated with the previous updates $\Gamma^{(n)}$:

[TABLE]

for $1\leq i\leq N$ and $0\leq n<n_{\mathrm{max}}$.

Another approach is to represent the quantizer by the column vector $\mathbf{\Gamma}$, derive the Hessian, $\nabla^{2}D(\mathbf{\Gamma})$, and compute the updated estimates of the quantizer using an iterative Newton-Raphson method

[TABLE]

for $0\leq n<n_{\mathrm{max}}$. We now develop this approach further.

Suppose $f_{X}$ and $F_{X}$ are the PDF and CDF of $X$, respectively. We define the $p$-th lower partial expectation as

[TABLE]

where $M_{X}^{0}(X)=F_{X}(x)$ represents the distribution function of $X$. Then, direct integration of the distortion function (1) gives

[TABLE]

Consequently, the elements of the vector $\nabla D(\mathbf{\Gamma})$ are given by

[TABLE]

for $1\leq i\leq N$.

Similarly, the tridiagonal Hessian matrix, $\nabla^{2}D(\mathbf{\Gamma})$, may be computed. It has diagonal elements given by

[TABLE]

respectively. Note that the quantities required to compute a Newton-Raphson iteration (i.e., the gradient and Hessian) only require the PDF, CDF and first lower partial expectation to be known. The second lower partial expectation is required only if one wishes to compute a numerical estimate of the distortion.

2.1 Efficient Implementation

We now provide a matrix formulation of the above Newton iteration intended to aid efficient implementation. As stated previously, the quantizer is represented by a column vector $\mathbf{\Gamma}$. This vector and three other column vectors required are defined by

[TABLE]

Defining $\mathbf{p}$ as a row vector is convenient since the expectation of a functional $H$ applied to the quantizer is

[TABLE]

where $H$ is applied element-wise to $\mathbf{\Gamma}$. Moreover, this will be compatible with a Markov chain formulation of the recursive marginal quantization technique presented later.

Using these vectors the gradient of the distortion function is then

[TABLE]

where $\circ$ indicates the element-wise Hadamard product.

The super- and sub-diagonal (or off-diagonal) entries of the Hessian matrix $\nabla^{2}D(\mathbf{\Gamma})$ are given by the length-$(N-1)$ row vector

[TABLE]

with the main diagonal given by

[TABLE]

where the copies of the $\mathbf{h}_{\mathsf{off}}$ vector are appended and prepended with a zero. It is now straightforward to set up the Newton-Raphson iteration in terms of the quantities.

2.2 Examples

In this section we apply the above methodology to the Gaussian distribution and the noncentral chi-squared distribution with one degree of freedom. The latter is important for the higher-order recursive marginal quantization schemes that we explore later in the paper.

2.2.1 The Standard Normal Distribution

When $X$ is a standard normal random variable we have

[TABLE]

where $\phi(\cdot)$ and $\Phi(\cdot)$ are the standard normal PDF and CDF, respectively. Here, a good guess for the initial quantizer $\mathbf{\Gamma}^{(0)}$ is

[TABLE]

for $1\leq n\leq N$. Figure 3 shows the quantizer of cardinality $N=50$, using this initial guess, after $n_{\mathrm{max}}=20$ Newton-Raphson iterations. At first glance, it is tempting to think of the quantizer (represented by bars) as one might think of a histogram and suspect that it does not adequately capture the features of the density because it has the “incorrect shape”. This is, however, a misleading analogy since a histogram represents the probability of realising a random variable in equally sized intervals on the $x$-axis. The quantizer, however, accumulates the probability mass over the Voronoi regions associated with each of the codewords and, as these regions are larger in the tails than in the body of the distribution, proportionally more probability mass is accumulated for codewords in the tails than for codewords in the body.

2.2.2 The Noncentral Chi-squared Distribution

While, in general, the noncentral chi-squared distribution must be specified using Bessel functions, this is not the case when the degree of freedom equals one. In particular, consider the random variable $X=(Z+\mu)^{2}$, where $Z\sim\mathscr{N}(0,1)$. Then $X\sim\chi^{\prime 2}(1,\lambda)$, is a noncentral chi-squared distributed with one degree of freedom and noncentrality parameter $\lambda=\mu^{2}$. Moreover, on $x\in\mathbb{R}^{+}$, we have

[TABLE]

This means that we may express the noncentral chi-squared distribution with one degree of freedom using the standard normal PDF and CDF, thus allowing efficient computation of a quantization scheme. This will be important for computational efficiency when we implement higher-order RMQ schemes later in the paper.

When implementing VQ, care must be taken when evaluating these functions at left and right limits. To ensure convergence we set $f_{X}(0)=F_{X}(0)=M_{X}^{1}(0)=0$, $f_{X}(\infty)=0$, $F_{X}(\infty)=1$ and $M_{X}^{1}(\infty)=1+\lambda$. Of course, all three functions are zero when $x$ is negative. A good initial guess for $\mathbf{\Gamma}^{(0)}$ is given by

[TABLE]

for $1\leq n\leq N$.

Figure 4 shows three examples of quantizers of cardinality $N=50$ for the noncentral chi-squared distribution with one degree of freedom for a range of noncentrality parameters. Note that for certain values of the noncentrality parameter (e.g. the central panel) the distribution is bimodal while for larger values it resembles a Gaussian distribution.

3 Recursive Marginal Quantization

Consider the continuous-time diffusion specified by the stochastic differential equation

[TABLE]

defined on the filtered probability space $(\Omega,\mathscr{F},(\mathscr{F}_{t})_{t\in[0,T]},\mathbb{P})$ with $a$ and $b$ sufficiently smooth and regular functions to ensure the existence of a weak solution. The question of interest is:

How does one optimally approximate $X_{t_{k}}:\Omega\rightarrow\mathbb{R}$, for some time discretisation point $t_{k}\in[0,T]$, when the distribution of $X_{t_{k}}$ is unknown?

Usually this is achieved by performing a Monte Carlo experiment using a discrete-time approximation scheme for the SDE, the simplest scheme being the Euler-Maruyama [Maruyama, 1955] update

[TABLE]

for $0\leq k<K$, where $\Delta t=T/K$ and independent $Z_{k+1}\sim\mathscr{N}(0,1)$, with initial value $\widetilde{X}_{0}=x_{0}$. The innovation of Pagès and Sagna [2015] was to show that a recursive procedure based on quantizing these updates is possible.

Since $\widetilde{X}_{1}$ has a Gaussian distribution, it is possible to use vector quantization to obtain $\Gamma_{1}$, an optimal quantization grid for the first step of the above scheme. One must, however, find a way to quantize the successive (marginal) distributions of $\widetilde{X}_{k+1}$. Given knowledge of the distribution of $\widetilde{X}_{k}$, the distortion of the quantizer $\Gamma_{k+1}$ may be written as

[TABLE]

Unfortunately, we do not know the exact distribution of $\widetilde{X}_{k}$ for $k>1$. The main result of Pagès and Sagna [2015] shows that if one uses the previously quantized distribution of $\widehat{X}_{k}$, instead of the continuous distribution of $\widetilde{X}_{k}$, the resultant procedure converges. Furthermore, the error associated with this procedure is bounded by a constant, which is dependent on the parameters used. As a result, the integral in (4) may be rewritten as a sum over the codewords in quantizer $\Gamma_{k}$ and their associated probabilities.

Then an approximate value for the distortion may be computed as

[TABLE]

Here, $N_{k}$ is the cardinality of the quantizer $\Gamma_{k}$ at time step $k$, which is allowed to vary. With this definition of D\big{(}\Gamma_{k+1}\big{)} we may now specify a Newton-Raphson iteration in order to compute the quantizer at time step $k+1$, which minimizes the distortion.

Given the quantizer at time $t_{k}$, represented as a column vector $\mathbf{\Gamma}_{k}$, and the associated probabilities, \mathbb{P}\big{(}\widehat{X}_{k}=\gamma_{k}^{i}\big{)} for $1\leq i\leq N_{k}$, the Newton-Raphson iteration for the quantizer $\mathbf{\Gamma}_{k+1}$, at time $t_{k+1}$, is given by

[TABLE]

for $0\leq n<n_{\mathrm{max}}$.

Before developing the mathematics further, we pause to provide an intuitive explanation of how the RMQ algorithm proceeds. Figure 5 is a depiction of the process that occurs. The top panel shows the quantizer at time step $k$. Conditional on each codeword, a Gaussian Euler update is propagated (second panel). In panel three, these updates are weighted by the probability of the associated originating codeword and summed to produce the marginal density at time step $k+1$, as shown in the final panel. The distribution associated with this marginal density is the distribution that is quantized to produce the quantizer at time step $k+1$. This process is repeated until the quantizer at the final time is produced.

We now proceed to derive the quantities required for the Newton-Raphson iterations. To summarize notation, we write the update in affine form as

[TABLE]

where

[TABLE]

with $Z_{k+1}^{i}\sim\mathscr{N}(0,1)$ identically distributed to $Z_{k+1}$. Here, we have introduced a new index $i$ for the random variable $Z^{i}_{k+1}$ anticipating that it may depend on $\gamma_{k}^{i}$. This is redundant in the case of the Euler update because $Z_{k+1}$ is a standard normal random variate irrespective of starting point. This more general notation will become necessary when we analyse more general cases (see Section 4). We denote the corresponding density and distribution functions by $f_{Z_{k+1}^{i}}$ and $F_{Z_{k+1}^{i}}$ respectively.

With this notation in place, the approximate marginal distribution of $\widetilde{X}_{k+1}$ is

[TABLE]

where $H(\cdot)$ is the Heaviside step function and $\,\mathrm{sgn}(\cdot)$ is the signum function. The last step follows due to the fact that the left-hand probability on the penultimate line may be written as

[TABLE]

It should be noted that, for the Euler update, $m_{k}^{i}$ is a proxy for the volatility of the SDE and its positivity is usually guaranteed, in which case (7) may be simplified. However, we persist with this formulation because in the general case $m_{k}^{i}$ may not be guaranteed to be positive.

The elements of the gradient of the distortion $\nabla D(\Gamma_{k+1})$ may then be written as

[TABLE]

where $1\leq j\leq N_{k+1}$ is the index tracking the elements of the $t_{k+1}$ quantizer, and the index associated with the $t_{k}$ quantizer is $i$. The integration bounds in (8) must now be expressed in terms of the variable of integration.

Here, as in Section 2, $U_{k+1}^{i}\in R_{j}(\Gamma_{k+1})$ is equivalent to the inequality

[TABLE]

and $r^{1-}_{k+1}=-\infty$ and $r^{N_{k+1}+}_{k+1}=\infty$ by definition. Defining the conditionally normalized region boundaries,

[TABLE]

allows the inequality to be written in terms of the random variable $Z^{i}_{k+1}$ as

[TABLE]

which can now be used as the range over which the integration is taken.

Directly evaluating (8), each element of the gradient of the distortion at time $t_{k+1}$ is given by

[TABLE]

respectively.

Although these expressions may appear complex, they are simply summations over the density function, cumulative distribution function and first lower partial expectation of the random variable, $Z^{i}_{k+1}$, and are thus easy to compute when these functions are known.

All the detail required to implement the Newton iteration (5) has now been provided with the exception of the initial guess. In all applications considered in this paper, we have assumed that $N_{k}=N$ for $1\leq k\leq K$, and used the quantizer from the previous time step as the initial guess, i.e., $\mathbf{\Gamma}_{k+1}^{(0)}=\mathbf{\Gamma}_{k}$.

3.1 Efficient Implementation

As in Section 2.1, where we provided an efficient matrix formulation for the Newton-Raphson iteration required for VQ, RMQ is also amenable to a matrix specification. This aids simple and computationally efficient implementation.

Aside from a guess for $\mathbf{\Gamma}_{k+1}$, we require the following time-indexed column vectors

[TABLE]

The row vector of probabilities

[TABLE]

is retained and a row-vector of ones of length $d$ is denoted by $\mathbf{j}_{d}.$ With the exception of $\Delta\mathbf{\Gamma}_{k+1}$, which must be recomputed before each Newton-Raphson iteration, the other vectors are computed once per time step.

Before each Newton-Raphson iteration, three matrices must be computed in terms of the new estimate of $\mathbf{\Gamma}_{k+1}$: an $N_{k}\times N_{k+1}$ matrix of transition probabilities

[TABLE]

The gradient of the distortion function at time step $k+1$ may then be written in terms of these vectors and matrices as

[TABLE]

where $\circ$ is the Hadamard (or element-wise) product.

The super and sub-diagonal elements of the (tridiagonal) Hessian matrix, $\nabla^{2}D\left(\mathbf{\Gamma}_{k+1}\right)$, are given by the vector

[TABLE]

while the main diagonal is given by

[TABLE]

where ${\circ-1}$ in the exponent refers to the element-wise inverse.

Equations (11), (12) and (13) provide the necessary components required for implementation of the Newton-Raphson iteration in (5). After the requisite number of iterations, the probabilities associated with the final quantizer $\mathbf{\Gamma}_{k+1}$ are computed using

[TABLE]

where $\mathbf{P}_{k+1}$ must be recomputed in terms of the final $\mathbf{\Gamma}_{k+1}$. Thus, the matrix formulation presented here allows RMQ to be interpreted as the propagation of an inhomogeneous, discrete time Markov chain, where $\mathbf{\Gamma}_{k}$ represents the Markov states at time step $k$, the probability of being in those states is $\mathbf{p}_{k}$, and the associated transition probability matrix is $\mathbf{P}_{k+1}$. Sometimes in the literature the transition probability matrix between time step $k$ and $k+1$ is represented as $\mathbf{P}_{k,k+1}$, we have chosen to omit the first index.

3.2 Example

As a first example of the RMQ algorithm, Figure 6 shows the evolution of the quantizers through time for the canonical (risk-neutral) geometric Brownian motion process

[TABLE]

using parameters $S_{0}=100$, $r=5\%$ and $\sigma=30\%$. The RMQ parameters used were $T=1$, $K=12$, $\Delta t=T/K$ and $N_{k}=200$ for all $k$, with $n_{\mathrm{max}}=50$ for the VQ algorithm and $n_{\mathrm{max}}=5$ for the RMQ algorithm. Unless otherwise stated, these parameters are used wherever geometric Brownian motion is used for other calculations in the paper. In these plots, the colour of the quantizer indicates the associated time step, with lines in blue closer to initial time and lines in green closer to final time. This convention is kept throughout.

4 Higher-order RMQ Schemes

Given the general formulation of the RMQ algorithm in the previous section, we now explore two high-order extensions: the Milstein scheme and a simplified weak order 2.0 scheme. Any numerical scheme for an SDE that can be written in the affine form (6), may be used with the RMQ algorithm as long as the CDF, PDF and lower partial expectation of the random variable $Z^{i}_{k+1}$ can be computed for all $i$ and $k$.

4.1 The Milstein Scheme

The Milstein [1975] scheme for the SDE in (3) is given by

[TABLE]

for $0\leq k<K$, where $\Delta t=T/K$ and $Z_{k+1}\sim\mathscr{N}(0,1)$ with initial value $\widetilde{X}_{0}=x_{0}$. By completion of the square, this may be written as

[TABLE]

Thus, the update (6) may be written in the affine form required as

[TABLE]

The random variable $Z_{k+1}^{i}$ is now noncentral chi-squared distributed with one degree of freedom and noncentrality parameter

[TABLE]

It is important to note that, unlike the Euler-Maruyama case, the distribution of the random variable $Z^{i}_{k+1}\sim\chi^{\prime 2}(1,\lambda_{k+1}^{i})$ now depends on the codeword $\gamma_{k}^{i}$.

Although the Milstein scheme possesses a strong order of convergence of $1$, compared to the Euler scheme, which only has strong order of convergence of $\frac{1}{2}$, both schemes have a weak order of convergence of $1$. Thus, while the Milstein scheme is more accurate in a strong sense than the Euler scheme, we require a different update for higher weak order convergence, which is needed when approximating expectations of financial payoffs. We now explore such a scheme.

4.2 A Weak Order $2.0$ Taylor Scheme

While it is not possible to write a weak order 2.0 Taylor scheme in the affine form required, the simplified weak order 2.0 scheme of Kloeden and Platen [1999] is amenable. This scheme is given by

[TABLE]

for $0\leq k<K$, where $\Delta t=T/K$ and $Z_{k+1}\sim\mathscr{N}(0,1)$ with initial value $\widetilde{X}_{0}=x_{0}$. Again, completion of squares is used to write this update in the required affine form,

[TABLE]

Here, $Z_{k+1}^{i}$ is again noncentral chi-squared distributed with one degree of freedom, with noncentrality parameter given by

[TABLE]

or, more succinctly, $Z^{i}_{k+1}\sim\chi^{\prime 2}(1,\lambda_{k+1}^{i})$.

4.3 Examples

To illustrate the accuracy of the above schemes, the RMQ algorithm is applied to geometric Brownian motion, as previously described by (14), and the constant elasticity of variance (CEV) process. The SDE for the CEV process is

[TABLE]

and, in the examples that follow, the process-specific parameters chosen were $S_{0}=100$, $r=5\%$ and $\alpha=0.7$ and $\sigma=\sigma_{\text{LN}}S_{0}^{1-\alpha}$, with $\sigma$ given in terms of the instantaneous log-normal volatility $\sigma_{\text{LN}}=30\%$. In the case of GBM the parameters in Section 3.2 were used. For both GBM and CEV, the RMQ specific parameters mentioned in that section were also used, unless stated otherwise.

In Figure 7 we show numerical evidence for weak order convergence. The absolute value of the difference between the first moment of the resulting terminal quantizer and the true moment is plotted for a range of time step sizes. In each case, the horizontal axis provides the base-2 logarithm of the step size and the vertical axis the base-2 logarithm of the error in the first moment. Thus, the slope of each graph reflects the power of the step size and hence the order of weak convergence for the error. In the figure legends, the regressed gradients of the graphs, denoted by $\beta$, indicate the weak orders of convergence. As is theoretically expected, see Kloeden and Platen [1999], both the Euler and Milstein schemes have approximately weak order one convergence, whereas the simplified weak order 2.0 scheme reaches a weak order close to two. Therefore, the latter scheme is by magnitude more accurate and can be expected to produce results with substantially lower error than the Euler scheme when valuing contingent claims. For increased accuracy in these plots, the cardinality used was $N_{k}=1\,000$.

Since the true conditional distributions for GBM and CEV are known in closed form, we can compare the approximate marginal distributions at time step $k+1$, as implied by the quantizer at time step $k$ and computed using (7), with the exact distributions at time step $k+1$. In the case of the CEV process we used the analytical expression for the distribution given by Lindsay and Brecher [2012] adjusted for the drift component in the SDE.

For each of the three schemes, under the GBM and CEV test cases, the difference between the exact marginal distribution and the implied marginal distribution is plotted in Figure 8. Here we reverted to using a cardinality of $N_{k}=200$. Note the scale of the $y$-axes of the graphs in the figure — from top to bottom, the magnitude of the error decreases by an order of magnitude in each successive row. This gives an indication of the improvement that can be expected when these higher order schemes are used to price contingent claims.

5 The Zero Boundary

Sometimes discrete-time approximations of an SDE may exhibit behaviour that is inconsistent with the true solution. For example, an Euler-Maruyama approximation of geometric Brownian motion or the CEV process can, under certain circumstances, generate negative values, even though the SDE specification guarantees non-negativity in each case. As a result, discrete-time Monte Carlo simulations are often modified to generate reflecting or absorbing behaviour at zero, see for example Lord et al. [2010]. In this section, we describe how the RMQ algorithm may be modified in a similar manner.

5.1 Absorbing Boundary

To model an absorbing boundary, the domain of the approximate marginal distribution of $\widetilde{X}$, see (7), must be left-truncated at zero. Implementing the RMQ algorithm with a left limit of zero results in a quantizer that has probabilities at each time step that do not sum to unity. The probability that is not accounted for as a result of the domain truncation is the mass accumulated at the absorbing zero boundary. To compensate for this, the quantizer at each time step can be augmented with an extra codeword, which has a value of zero and a probability equal to one minus the sum of the probabilities associated with all the other codewords at that time step. The transition probability matrix may also be augmented in a consistent manner by realising that once the process attains the zero state it must remain in that state indefinitely, i.e., the conditional probability of moving from the absorbing state to any other state is zero, and, correspondingly, the conditional probability of remaining in the absorbing state is one.

Modifying the algorithm is straightforward and incurs no additional computational burden. Given that the elements of the previous quantizer $\mathbf{\Gamma}_{k}$ are all positive, the affine form of the update (6), will be negative when

[TABLE]

This implies that the domain of each $Z_{k+1}^{i}$ must be left-truncated at $-\tfrac{c_{k}^{i}}{m_{k}^{i}}$ to ensure only positive codewords at time-step $k+1$. This is achieved by setting

[TABLE]

for $1\leq i\leq N_{k}$, in the implementation described in Section 3.1. This is equivalent to assuming $r^{1-}_{k+1}=0$ in (9). The rest of the algorithm proceeds without modification.

Of course, this all depends on the fact that the quantizer at the first time step, $\mathbf{\Gamma}_{1}$, also has positive elements. This is achieved using an analogous truncation in the vector quantization algorithm. The initial guess for the Newton iteration must also ensure positivity.

5.2 Reflecting Boundary

Figure 9 shows $f(x)$, a density function — in this case a standard Gaussian density. The red line represents a reflecting boundary at $\bar{x}=-1.5$. The $f$ values to the left of the boundary are reflected and depicted by the dashed yellow line in the figure. These values are given by $f(2\bar{x}-x)$ for $x>\bar{x}$.

Thus, restricting the domain to $[\bar{x},\infty)$ the reflected density, denoted $\bar{f}(x)$, is given as the sum

[TABLE]

Direct integration of this expression over the integration limits from $\bar{x}$ to $x\in[\bar{x},\infty)$ gives the reflected distribution function

[TABLE]

and the first lower partial expectation function

[TABLE]

where $F(x)$ and $M^{1}(x)$ are the un-reflected distribution and first lower expectation functions associated with $f$. This may be applied, not only to the Gaussian case, but also to the noncentral chi-squared cases required for the higher order updates.

Modifying the RMQ algorithm to allow for a reflecting boundary at zero requires two changes to the implementation described in Section 3.1. Firstly, the lower bound for the integration, i.e., the domain of the $Z_{k+1}^{i}$ random variable in each affine update, must be left-truncated by replacing the furthest left region boundary as in (15) above. Secondly, the density, distribution and first lower partial expectation functions associated with each random variable, must be replaced by their reflected counterparts

[TABLE]

for $x\in[\bar{x}^{i}_{k},\infty)$, where $\bar{x}^{i}_{k}=-\tfrac{c_{k}^{i}}{m_{k}^{i}}$. The remainder of the algorithm proceeds as normal. The astute reader will have noticed that there are two terms missing in (17) when compared with (16). The reason for this omission is that these terms are constants for each $i$ and that the RMQ algorithm always only requires differences of partial moment terms, as seen in (10). There is, therefore, a cancelation of the constant terms when this difference is taken and thus we may use a definition that excludes them.

As in the case of the absorbing boundary, the analogous reflection must be applied in the vector quantization algorithm to ensure that $\mathbf{\Gamma}_{1}$ is consistent.

5.3 Examples

It is well known that when $0<\alpha<0.5$, the CEV process may reach zero and that this state may be either absorbing or reflecting. Lindsay and Brecher [2012] give the corresponding marginal distributions for both these cases (it is easy to adjust their formulations to account for the drift term in the SDE for the CEV process). In Section 4.3 we considered the CEV process with $\alpha>0.5$, which only allows absorption at zero. Now consider the case where $S_{0}=0.5$, $\alpha=0.35$ and $\sigma_{\text{LN}}=50\%$, with the rest of the parameters as before. Figure 10 shows the difference between the exact marginal distribution and the marginal distribution implied by RMQ for the three schemes as modified to account for an absorbing boundary (left) and a reflecting boundary (right). As before, the scale of the graphs changes from top to bottom, indicating the improvement as a result of the choice of the schemes.

Note that, under this choice of parameters for the CEV model, the standard RMQ formulation of Section 3 fails. Without implementing the modifications for either absorption or reflection proposed in this section, some codewords become negative at a certain point in the execution of the RMQ algorithm, leading to discrete-time updates with imaginary values.

6 Pricing

In this section, contingent claims are priced using the RMQ algorithm and the three update schemes are compared for accuracy. The claims priced include European, Bermudan and discretely-monitored barrier options under the dynamics of both geometric Brownian motion (GBM) and its generalization, the constant elasticity of variance (CEV) model.

The GBM model and its parameters are described in Section 3.2, whereas the specification for the CEV model can be found in Section 4.3. These parameters are used for pricing throughout this section. As is implied by these specifications, the continuously compounded interest rate is assumed constant with a value $r=5\%$. All option maturities are one year and the RMQ algorithm is executed using $K=12$, i.e., using monthly steps, with constant cardinality of $N_{k}=200$ for all $k$.

6.1 European Option Pricing

Once a terminal quantizer has been obtained using the RMQ algorithm, a European option with payoff function $H(S,X)$ at maturity $T=t_{K}$, where $S$ represents the asset process and $X$ the strike, may be priced directly by using the expectation defined in (2). The price is given by

[TABLE]

where $H_{0}$ is the value of the claim at initial time $t_{0}=0$ and $H(\mathbf{\Gamma}_{K},X)$ is the function $H$ applied element-wise to $\mathbf{\Gamma}_{K}$, which, as specified previously, is a column vector of length $N_{K}$.

Figure 11 shows the accuracy of put option prices for the GBM and CEV models for a wide range of strikes. The GBM option prices are compared against the Black-Scholes option pricing formula, whereas the CEV prices are compared against the analytical solution originally due to Schroder [1989] and reformulated in terms of the noncentral chi-squared distribution by Hsu et al. [2008]. In the graphs, the $x$-axis represents fixed-spot inverse moneyness, which is determined as the variable strike value over the initial asset price, $S_{0}$.

Even though the Euler scheme is reasonably accurate to start with, the increased accuracy of the Milstein and the simplified weak order 2.0 schemes is evident. For certain strikes the error is reduced by an order of magnitude.

6.2 Bermudan Option Pricing

Bermudan option prices are computed using the standard Backward Dynamic Programming Principle (BDPP), an important result from discrete-time optimal stopping theory. Pagès [2014] reviews the use of the BDPP as applied to grids that result from a quantization.

Once quantization grids and corresponding transition probability matrices have been computed using the RMQ algorithm, the high-level algorithm for Bermudan option pricing may be specified as follows:

Initialize $\mathbf{h}_{K}=H(\mathbf{\Gamma}_{K},X)$ 2. 2.

For $k=K-1,\dots,1$

      Set $\mathbf{h}_{k}=\max(H(\mathbf{\Gamma}_{k},X),e^{-r\Delta t}\mathbf{P}_{k+1}\mathbf{h}_{k+1})$ 3. 3.

Set $H_{0}=e^{-r\Delta t}\mathbf{p}_{1}\mathbf{h}_{1}$

Here the $\max$ function is applied element-wise with its second argument being the continuation value, which is easily computed as a conditional expectation due to availability of the transition probability matrix at each time step. The initial value of the Bermudan claim is given by $H_{0}$.

In Figure 12 the accuracy of a Bermudan put option with monthly exercise opportunities is shown for the GBM and CEV models. The reference price is computed using a high resolution Crank-Nicholson finite difference scheme using 600 time steps and 800 stock increments, equally spaced between zero and $4\times S_{0}$.

All three RMQ algorithms result in low absolute errors, with the simplified weak order 2.0 scheme again producing errors that are an order of magnitude smaller.

6.3 Barrier Option Pricing

The pricing of barrier options has previously been explored in the context of quantization by Sagna [2011]. This work showed that the barrier-crossing approach described in Section 6.4 of Glasserman [2003] may be applied to marginal quantization using a so-called transition kernel formulation. Using our notation, we now combine this approach with our proposed more accurate schemes and apply it to discretely monitored barrier options.

Consider expression (18) for pricing European options, which may be re-written as

[TABLE]

To price a knock-out barrier option the transition probability matrix at each time step in this expression must to be modified to take into account the possibility that the underlying process breaches the barrier. Thus, we rescale the transition probabilities by multiplying them by the probability of not having crossed the barrier.

Let $g(x,y)$ be the probability of transitioning between states $x$ and $y$ without crossing the barrier. If we form an $N_{k}$ by $N_{k+1}$ matrix of values

[TABLE]

then $\mathbf{P}_{k+1}\circ\mathbf{G}_{k+1}$ defines the transition kernel. The barrier option may then be priced using

[TABLE]

where $\mathbf{g}_{1}=[g(S_{0},\gamma_{1}^{1}),\dots,g(S_{0},\gamma_{1}^{N_{1}})]$ is a row vector.

In the case of discretely-monitored up-and-out barrier options with barrier level $L$, the function $g$ is given simply as the indicator function

[TABLE]

We refer to Glasserman [2003] and Sagna [2011] for the continuous monitoring case using the Rayleigh distribution.

In Figure 13 the accuracy of discretely-monitored up-and-out put option prices generated using RMQ is compared to a Monte Carlo implementation under the GBM and CEV models. The barrier levels ($x$-axis) are expressed as multiples of the at-the-money strike. Since we have chosen $K=12$ the barrier is monitored monthly.

The reference prices are provided by a one million path Monte Carlo experiment. The Monte Carlo paths are generated using Euler-Maruyama updates with $1\,200$ time steps, while ensuring that the barriers are only monitored at monthly intervals. We used the exact transition density to generate Monte Carlo samples for GBM to confirm that results were consistent and generating the correct standard deviations.

The results show a similar pattern to those in the previous sections, with an important caveat: the simplified weak order 2.0 scheme produces prices that, for the majority of the barrier values considered, lie within the three standard deviation bound of the million-path Monte Carlo experiment. The other two RMQ schemes are producing results that are statistically significantly incorrect when compared to the Monte Carlo simulation.

7 Conclusion

In this paper, the recursive marginal quantization methodology of Pagès and Sagna [2015] has been extended from the standard Euler-Maruyama scheme to higher-order numerical schemes, specifically the Milstein scheme and the simplified weak order 2.0 scheme of Kloeden and Platen [1999]. This entailed introducing noncentral chi-squared updates and generalising the formulation of Recursive Marginal Quantization (RMQ).

Is is also shown how to augment the RMQ algorithm in order to implement absorption or reflection at the zero boundary, thus ensuring non-negativity of solutions. This allows RMQ to be applied in important cases, where the algorithm may previously have failed.

Improved approximation of the marginal distributions by the proposed new methods has been demonstrated, using geometric Browninan motion and constant elasticity of variance dynamics. All the schemes were used successfully to price European, Bermudan and discrete barrier options. The pricing results show that once RMQ is implemented with second weak order updates, it provides extremely accurate contingent claim pricing, well suited for the extremely fast calibration of entire derivative books.