# Boosting convergence rate of numerical simulations for Monte-Carlo and PDE (Partial Differential Equations) based methods

## Introduction

In this post, we briefly present some numerical results that we found interesting for numerical simulations using Monte-Carlo, but also for PDE (Partial Differential Equations) methods.

These results are based on a new numerical method which allows us to explicitly compute Monte Carlo samples for a very large class of stochastic processes: we can produce these samples for any stochastic processes defined by an SDE (Equation to stochastic derivatives ), and also in very high dimensions (we tested up to 64 dimensions).

In the sequel, we illustrate the qualitative properties of these samples. We illustrate also the performance of this sampling method, when used together with a PDE solver, studying an application that comes from financial mathematics: the computations of prices, Greeks (deltas, gammas), as well as some regulatory risk measures, that are VAR (Value At Risk), or EE (Expected Exposure), for a large portfolio of AutoCall.

## Qualitative Property of sampling.

The computed samples have the particularity of being quite regularly distributed, recalling the distribution that would be obtained by optimal quantification methods. In particular, for the computation of integral (expectation), we measure convergence rates on certain classes of functions that are better than the Monte Carlo methods classically used in industry, from a scale factor.To illustrate the regularity of this distribution, we present in the following figures three samples of N = 256 points using three sampling methods. The stochastic process studied is a bivariate log-normal martingale process, negatively correlated (-0.5), each of the two marginal distributions having a volatility of 10%, the process is studied at time t = 1.

The first figure at left-hand side (labeled MT19937) represents the generation of 256 samples computed by a pseudo-random generator Mersenne Twister, widely used in industry. This generator has a well known convergence rate of statistical nature:

where the constant appearing on the right side of the previous equation is the variance of P, considered as a random variable at time 1.

The second figure (labeled Sobol) at center represents the generation of 256 samples computed from a Sobol generator, which is a quasi-random generator used industrially. This generator, motivated by a better distribution of sampling points, has a rate of convergence which is still conjectured to our knowledge:

where the constant appearing on the right hand side is the bounded variation of the sampled function P, N the number of sampling points, and D the dimension. These convergence rates are better than a pseudo-random generator when the sample number N is large. However, the quality of the sample deteriorates rapidly as the dimension of problem D increases.

The third figure to the right (labeled CoDeFi) represents the generation of 256 samples computed from this new method. Numerically, these samples present some interesting properties, besides the striking regularity of their distribution. First, their convergence rate is always at least as good as the two other sequences (MT19937 and SOBOL) on all tests we performed. On the other hand, these samples exhibit convergence rates of kind

This rate of convergence is measured for a set of functions P that are quite widely used, particularly in mathematical Finance. The constant at right-side measures a form of regularity of P and which, numerically, seems independent of the dimension D. This is a surprising result: four orders of magnitude above a Monte Carlo method for certain functions means, for instance, that a sample of 100 points compute expectations as accurately as 100 million sampling using a MT19937 generator, regardless of the number of dimensions. This could therefore be an appreciable gain in terms of algorithmic complexity. Of course, the counterpart is that these sampling sequences are more difficult and takes more time to compute than pseudo or quasi-random type sequences.

## Application to PDE methods :a regulatory risk type computation.

The result of the test presented here was carried out with CoDeFi, a framework allowing the valuation of complex derivative products, or portfolios, as well as the computations of their risk indicators, or evaluation of complex strategies as can be risk limits.

CoDeFi is based on a hybrid PDE / Monte-Carlo technique. The main idea of this hybrid PDE solver is to use Monte-Carlo trajectories to produce computing grids used by the PDE solver. In particular, we can use as grid the trajectories produced by a generator MT19937, or Sobol, but CoDeFi also proposes to calculate its own trajectories. For instance, the previous figure corresponds to three grids usable by CoDeFi at time 1, in dimension 2, for PDE computations. In this test, we use the grid produced by CoDeFi, that is corresponding to the third figure for the two-dimensional case (CoDeFi label). The purpose of this test is to qualify these samples for the PDE engine.

This test is specified as follows: first define a univariate (one-dimensional, D=1) stochastic process, and a function defined over the process filtration:- The process is a log-normal univariate martingale process, 10% volatility.- The function under consideration is a widely sold financial instrument, called an AutoCall. Its characteristics are: maturity five years, date of observation of annual maturity (1,2,3,4,5 years). The AutoCall barrier is 100%, the coupon barrier is 80%. Once this setting done, the test presented depends only on two parameters: D the dimension, N the number of points. For a given dimension D, a D-dimensional process is first defined in the following way:

– Each marginal process is defined as a copy of the above univariate process.

– A constant correlation is imposed between each marginal.

Then we define an AutoCall portfolio as follows:

– First define all univariate Autocall, as a copy of the Autocall defined above, by changing the underlying. There is therefore D univariate Autocalls in this portfolio.

– Second, we define all bi-variate worst-off Autocalls, as a copy of the Autocall defined above, by selecting all combination of two underlyings. There is therefore D (D-1) / 2 bi-variate Autocalls in this portfolio.

– Etc..until the definition of the last Autocall, which is a worst-of AutoCall written over D underlyings.

This portfolio contains exactly 2 power D Autocalls. However, in the context of this test, only D Autocalls are actually really different: each d-varied AutoCall, with d between 1 and D. Hence we can numerically approximate prices and derivatives (deltas) of each of these AutoCall by a Monte method-Carlo using a large number of trajectories to benchmark our results.

Finally, we pick up a number N, generate the computation grid and solve the retrograde Kolmogorov equation (also called the Black and Scholes equation) using a PDE grid of N samples. For instance, the following figure illustrates the case D = 10, N = 512: we plot the results using the ten prices computed by Monte-Carlo on the abscissa, since ordinates is used for the 1024 prices calculated by our PDE solver.

For this test, we measured two things:

– The overall error on the entire portfolio of 1024 AutoCall is 0.2% (relative error on price), which corresponds to a convergence factor of rate 1 / N with N = 512.

– The time consumed by the algorithm to compute the prices of these 1024 AutoCalls, as well as their 1024 * 10 first-order sensitivities (delta), as well as the 1024 * 100 second-order sensitivities (gammas), as well as a simple market VAR (Value At Risk) of this portfolio, as well as the EE (expected exposure) of this portfolio, using a half-core of an Intel I7 processor, is 20 seconds. It is also noted that the execution times of this test, for dimensions 11 (2048 AutoCalls), 12 (4096 Autocalls), and 13 (8192 AutoCalls), are below 30 seconds, with stable convergence rates.

## Conclusions.

We chose to illustrate the algorithmic performance gains obtained by these sampled sequences in a case that is disadvantageous[i]. However, even in this case, the proposed method shows a quite appreciable gain in terms of performance for regulatory-type calculations[ii].

[i] Indeed, the function corresponding to an AutoCall is a bounded variation function, a function class for which we know that the convergence rates of a sampling method can not exceed 1 / N, not 1/N^2

[ii] In particular, we do not know of any other numerical method capable of performing this test, with this level of precision (0.2% on prices, on VAR and EPE), of course for the computational power at our disposal.

[i] Indeed, the function corresponding to an AutoCall is a bounded variation function, a function class for which we know that the convergence rates of a sampling method can not exceed 1 / N. We could provide such tests for more regular functions showing 1/N^2 convergence rate.

[ii] In particular, we do not know of any other numerical method capable of performing this test, with this level of precision (0.2% on prices, on VAR and EPE), of course for the computational power at our disposal