Remarks on stochastic automatic adjoint differentiation and financial models calibration

TL;DR

This paper examines the use of Automatic Adjoint Differentiation (AAD) for calibrating stochastic models, highlighting its SIMD parallelization capabilities and validating theoretical findings with numerical experiments.

Contribution

It demonstrates the SIMD parallelization of AAD in stochastic model calibration and confirms the theoretical results through numerical experiments.

Findings

AAD enables perfect SIMD parallelization in calibration tasks.

Theoretical analysis aligns with numerical experiment results.

Provides computational cost assessment for AAD in this context.

Abstract

In this work, we discuss the Automatic Adjoint Differentiation (AAD) for functions of the form $G=\frac{1}{2}\sum_1^m (Ey_i-C_i)^2$, which often appear in the calibration of stochastic models. { We demonstrate that it allows a perfect SIMD\footnote{Single Input Multiple Data} parallelization and provide its relative computational cost. In addition we demonstrate that this theoretical result is in concordance with numeric experiments.}

Equations

$$\left\{\lambda \cdot \frac{\partial f(x)}{\partial x_i},i=1\dots n\right\}$$

$$G=\frac{1}{2}\sum _1^m(E{y}_i-C_i{)}^2,$$

$$F:{\mathbb{R}}^{M+N}\to {\mathbb{R}}_y^m,$$

$$c=8\times Number\_of\_Cores.$$

$$\mathrm{Cost}({\mathrm{F}}_{\mathrm{v}})=\frac{{\mathrm{K}}_{\mathrm{F}}}{\mathrm{c}}\mathrm{Cost}(\mathrm{F}),\mathrm{Cost}({\mathrm{R}}_{\mathrm{v}})=\frac{{\mathrm{K}}_{\mathrm{R}}}{\mathrm{c}}\mathrm{Cost}(\mathrm{F})$$

$$Ey\sim \frac{1}{number\_of\_Paths}\sum _{k=1}^{number\_of\_Paths}y({\omega }_k),$$

$$\mathrm{Cost}(\{\mathrm{Ey}\})=\frac{{\mathrm{K}}_{\mathrm{F}}}{\mathrm{c}}\times number\_of\_Paths\times \mathrm{Cost}(\mathrm{F}).$$

$$\frac{\partial G}{\partial x_j}=E\left(\sum _{i=1}^m(E{y}_i-C_i)\frac{\partial y_i}{\partial x_j}\right),j=1,\dots ,M.$$

$$\left\{\sum _{i=1}^m(E{y}_i-C_i)\frac{\partial y_i(\omega )}{\partial x_j},j=1\dots M\right\}$$

$$\mathrm{Cost}(\mathrm{G})=\frac{2{\mathrm{K}}_{\mathrm{F}}+{\mathrm{K}}_{\mathrm{R}}}{\mathrm{c}}\times number\_of\_Paths\times \mathrm{Cost}(\mathrm{F}).$$

$$\mathrm{Cost}(\mathrm{G})=\frac{2{\mathrm{K}}_{\mathrm{F}}+{\mathrm{K}}_{\mathrm{R}}}{8\times number\_of\_Cores}\times number\_of\_Paths\times \mathrm{Cost}(\mathrm{F}).$$

$$\left\{\begin{array}[]{l}dS_{t}=\mu S_{t}dt+\sqrt{V_{t}}L(t,S_{t})S_{t}dW_{t}^{S},\\ dV_{t}=\kappa(\theta-V_{t})dt+\xi\sqrt{V_{t}}dW_{t}^{V},\\ dW_{t}^{S}dW_{t}^{v}=\rho dt.\end{array}\right.$$

$$K_F/c$$

$$y:=x_N$$

$$x_m:=f_m(x_{{\tau }_{m}1},\dots ,x_{{\tau }_m{i}_m}),1\le m<N.$$

$$\overline{D}_{j}=\left\{\begin{array}[]{ll}\overline{D}_{j}+\overline{D}_{m}\frac{\partial f_{m}}{\partial x_{j}},&m\neq k,\\ \overline{D}_{j}+E\overline{D}_{m},&m=k.\end{array}\right.$$

$$E\frac{\partial y}{\partial x_i}=E{\overline{D}}_i.$$

Full text

Key words: Automatic Adjoint Differentiation, automatic vectorization, Single instruction multiple data, AAD-Compiler

The Automatic Adjoint Differentiation (AAD) is a rapidly growing field with a wide range of applications including image restoration [5], computer vision [9] and maching learning in general see e.g. [2], [10]. The long list of AAD applications can be found in the site of Autodiff community, see [1].

For the additional reading one can mention monograph by Antoine Savine [12], articles by L.Capriotti [4] and by other authors e.g. [13], [14].

The AAD has become a widespread tool in applications due to the following property: If one has an algorithm for a function $f:\mathbb{R}^{n}_{x}\to\mathbb{R}^{m}_{y}$, then AAD maps each $\lambda\in\mathbb{R}^{m}$ to the linear combinations of $\partial_{i}f$

[TABLE]

the computation cost does not exceed that of $f$ which is usually between 2 and 10 (depending on a specific AAD tool used).

The utilization of parallel computations appears naturally in calculation of expectations since it assumes the numerous independent calculation of an integrand. Suppose that one needs to calculate the $\frac{d}{dx}Ey$, the algorithm for the $\frac{d}{dx}y$ is first determined and then the average using parallel computations. From that point of view, the AD should avoid differentiate expectations. However, in case of the Adjoint AD, this approach becomes mandatory. The point is that the adjoint differentiation algorithm is not local (e.g. [4]) i.e. it requires analysis of the whole algorithm for $f$, thus the differentiation of expectation requires much more memory than taking expectation of the derivative.

Some important problems involve computations that include an expectation as an intermediate operation. It is not a trivial question how one can use AAD and avoid differentiation of expectations. In a recent article [7], Fries suggest a general recipe for this problem (see Appendix). However, Fries did not provide the accurate analysis for the computational cost of the proposed algorithm. This work conducts a step-by-step computation of the following simple but important example.

Consider a functional of the type

[TABLE]

where $y_{i}$ are random variables on a filtered probability space $(\Omega,\mathbb{Q},\{\mathcal{F}_{t}\})$ and $C_{i}$ are some given target conditions, $i=1,\ldots n$. It is assumed that we are provided with a *forward * algorithm

[TABLE]

which calculates $y_{i}$ for a given set of $M$ parameters and $N$ independent random variables in frames of a time-discretization of the stochastic process.

We assume that the AAD tool provides the algorithm $R$ (so called ”reverse” algorithm, e.g. [4]) which calculates the full set of derivatives $(\ref{0711A})$222it can be applied only after the $F$ was launched. We also assume that the AAD tool produces parallelized versions $F_{v}$ and $R_{v}$. This means that $F_{v}$ (or $R_{v}$) provides the execution of $F$ (or $R$) for $c$ independent sets of input data. For example, the natural value of $c$ for the case of an AAD-tool tuned to the Intel AVX512 architecture is

[TABLE]

Is noted that the first factor in (3) equals to $1$ in all known to us C++ AAD tools.

We consider the operations $F\to F_{v}$ and $F\to R_{v}$ as a simple way to parallelize calculations and we are interested to check if there is a way to use this tool effectively.

Correction coefficients $K_{F}$ and $K_{R}$ are defined via the following expression

[TABLE]

and they reflect the quality of the $AAD$ tool. Ideally, $K_{F}=1$, but in practice, the different software optimization and different hardware specifics can never make it perfect. One reason is that the executable version of $F$ is produced with all compiler’s optimization abilities turned On, including the ability SIMD vectorization.

Some remarks on the Monte-Carlo simulations. We approximate expectations by

[TABLE]

where the set $\{\omega_{i}\}$ contains Monte-Carlo simulations of a given stochastic process and each $\omega_{i}$ is a simulated sample path of sequences of random variables. It is assumed that drawings are $\mathbb{Q}$-uniform.

In particular,the cost of the independent evaluation for the set $\{Ey_{i},i=1,\ldots,m\}$ is

[TABLE]

Calculation of the gradient of $G$ (introduced in (2)). For any parameter $x_{j}$ we get that

[TABLE]

It leads to the following algorithm :

Calculate $Ey_{i}$ using $F_{v}$. The costs of this calculation is given by (5). 2. 2.

Fix a path $\omega$. Using formula (1) for the vector $\lambda$ with components $\lambda_{i}=Ey_{i}-C_{i}$ we get the set

[TABLE]

For $c$ paths, it costs $(K_{F}+K_{R})\rm{Cost}(F)$ since for each path $\omega$ the reverse algorithm can be executed only after the forward algorithm has been launched, unless execution results of the first step can be stored in the memory for each $\omega$. In case memory usage is not constrained, the cost is $K_{R}\rm{Cost}(F)$ 3. 3.

Summation over paths $\omega$. To calculate the cost of this operation we take into account that the integrand should be calculated $\rm{number\_of\_Paths}$ times.

Summarizing, we get the total cost

[TABLE]

For C++, AAD-compiler produced by MathLogic LTD can be rewritten as

[TABLE]

where $8$ reflects that the AVX512 architecture allows $8$ doubles per vector register.

Test on the Heston Stochastic Local Volatility model calibration. For financial institutions, the model calibration is a routine day-to-day procedure. Although the Partial Differential Equation (PDE)-based techniques became quite popular e.g. [6], [11], practitioners more often use direct numeric simulations due to the simplicity of the technique.

Consider the Heston SLV model given by a process

[TABLE]

Our implementation utilizes the standard Euler discretization (e.g. 3.4 of [8]) and European options in the quality of the $y_{i},i=,\ldots,m$. The set of optimization parameters consists values of the piecewise constant Leverage function $\{L(t_{i},S_{j})\}$ (where the set of interpolation nodes $\{t_{i},S_{j}\}$ is fixed) and five standard Heston model parameters. We applied the AADC333C++ AAD-tool produced by matlogica.com by MathLogic LTD and observed the following values for coefficients $K_{F}$ and $K_{R}$ defined in (4).

[TABLE]

Tests were executed on one core. The values of the aforementioned coefficients became almost constant when number of time intervals and number of optimization instruments $m$ grow. As mentioned, all optimization parameters (like vectorization) were turned on while the evaluation of the $\rm{Cost(F)}$.

Conclusion. We consider a situation when a straighforward application of the AAD software is not possible i.e. for functions of the form $G=\frac{1}{2}\sum_{1}^{m}(Ey_{i}-C_{i})^{2}$. We propose a two-step algorithm which leads to the estimate (6) and therefore, provides a perfect SIMD parallelization. A numeric implementation of the algorithm is realized for the Heston model using the AAD-Compiler by matlogica.com. Presented results are in good concordace with (6).

Appendix. Expected Backward Automatic Differentiation Algorithm following [7]. Consider a scalar function $y$ given by a sequence of operations:

[TABLE]

where the number of variables $i_{m}$ is either 0 and it means that the $x_{m}$ is an independent variable or $1\leq i_{m}<m$. Evidently for function $\tau_{m}$ we have $1\leq\tau_{m}<m$. We assume that $k-$th operator is an expectation operator and others $f_{m},m\neq k$ given by a closed-form444i.e. it can be evaluated in a finite number of the ”well-known” operations, see e.g. wikipedia https://en.wikipedia.org/wiki/Closed-form_expression for the detailed definition..

Sequentially,

• •

Initialise $\overline{D}_{N}=1$ and $\overline{D}_{m}=0,~{}m<N$.

• •

For all $m=N,N-1,...,1$ (iterating backward through the operator list)

• –

for all $j=\tau_{m}1,...,\tau_{m}i_{m}$ (iterating through the argument list)

[TABLE]

Then, for all $1\leq i\leq N$

[TABLE]

Acknowledgments. E.L. was partially supported by Portuguese funds through the CIDMA - Center for Research and Development in Mathematics and Applications and the Portuguese Foundation for Science and Technology (“FCT–Fundção para a Ciência e a Tecnologia”), within project UID/MAT/ 0416/2019.