[C] left=1.3in, right=1.3in, top=1in,bottom=1in, includefoot,
headheight=13.6pt

#  On Randomization of Affine Diffusion Processes
with Application to Pricing of Options on VIX and S&P 500

LECH A. GRZELAK<EMAIL_ADDRESS>Mathematical Institute, Utrecht University,
Utrecht, the Netherlands  Financial Engineering, Rabobank, Utrecht, the
Netherlands

###### Abstract

The class of Affine (Jump) Diffusion [8] (AD) has, due to its closed form
characteristic function (ChF), gained tremendous popularity among
practitioners and researchers. However, there is clear evidence that a
linearity constraint is insufficient for precise and consistent option
pricing. Any non-affine model must pass the strict requirement of quick
calibration- which is often challenging. We focus here on Randomized AD (RAnD)
models, i.e., we allow for exogenous stochasticity of the model parameters.
Randomization of a pricing model occurs outside the affine model and,
therefore, forms a generalization that relaxes the affinity constraints. The
method is generic and can apply to any model parameter. It relies on the
existence of moments of the so-called randomizer- a random variable for the
stochastic parameter. The RAnD model allows flexibility while benefiting from
fast calibration and well-established, large-step Monte Carlo simulation,
often available for AD processes. The article will discuss theoretical and
practical aspects of the RAnD method, like derivations of the corresponding
ChF, simulation, and computations of sensitivities. We will also illustrate
the advantages of the randomized stochastic volatility models in the
consistent pricing of options on the S&P 500 and VIX.

###### keywords:

Randomized Affine-Diffusion, RAnD Method, Stochastic Parameters, Heston,
Bates, VIX, S&P 500, Short Maturity Options, Stochastic Collocation,
Quantization.

††footnotetext: The views expressed in this paper are the personal views of
the authors and do not necessarily reflect the views or policies of their
current or past employers. The authors have no competing interests.

##  1 Introduction

The class of Affine Jump Diffusion, which we abbreviate by AD, asset price
models constitutes a family of processes that, under certain linearity
conditions, allow for a semi-closed form characteristic function (ChF) [8].
Once the ChF is derived, it can be used, via Fourier inversion, for pricing
certain financial derivatives, thus facilitating fast asset price model
calibration. Over many years, various affine models have been studied, leading
to significant contributions like efficient large-step Monte Carlo
simulations, PDE discretizations, quick calibrations or even analytical,
perturbation-based approximations for implied volatilities. These
accomplishments popularized the class of AD processes.

Although the high speed of “regular” option pricing is a necessary condition
that allows for calibration, it is insufficient for a model to be used for
exotic derivative pricing. The main limitation of classical dynamic stochastic
volatility models driven solely by possibly correlated Brownian motion is the
limited flexibility for short maturities, i.e., the slope of the skew is
insufficient to fit the market data. By adding jumps, this effect can, to some
extent, be mitigated; however, it is often unsatisfactory, especially in a
volatile market. The class of affine models is inherently inflexible, with
limited possibilities for model extensions. In particular, with the rise in
popularity of products such as options on volatility, it is evident that the
standard affine, stochastic volatility models are insufficient for accurate
and consistent option pricing. This observation has initiated a gold rush to
discover new classes of models that would address this pricing problem.

One possible route to extend the existing affine models was suggested by P.
Carr and L. Wu in [7], where the problem of insufficient skew was reported,
and the remedy of the form of parameter randomization was suggested: “Starting
from the jump-diffusion stochastic volatility model of Bates [2], it would be
tempting to try to capture stochastic skewness by randomizing the mean jump
size parameter or the correlation parameter (…) However, randomizing either
parameter is not amenable to analytic solution techniques that greatly aid
econometric estimation.” Although the authors considered randomization as a
technique to incorporate more flexibility into a stochastic model, thus a way
to improve the calibration quality, the concept of randomizing is more
fundamental, i.e., it represents the incorporation of the uncertainty of
potentially hidden states that are not adequately captured by deterministic
parameters. In this article, we focus on the randomization of the class of
affine models and allow for additional flexibility by letting the model
parameters be stochastic. In particular, we will randomize the parameters of
the Bates model and examine their impact on the shapes of implied
volatilities.

Attempts to randomize individual model parameters are already known in the
literature. Starting from the simplest model, the Black-Scholes model with the
randomized volatility coefficient- based on a discrete distribution- was
presented in [4]. An extension to continuous distributions was analyzed in
[19]. In their work, the authors discussed various stochastic variants for the
volatility and their level of explosiveness\- this led to several valuable
expansions for short maturity options. In [20, 25], for example, it was shown
that randomization of the initial variance of the Heston model significantly
increased the steepness of the implied volatility, primarily when the
randomizer is defined in an unbounded domain. This idea was further studied in
[24] where the term stationary, to describe random initial parameter, was
used. In [27] the concept of Recursive Marginal Quantization for one-
dimensional diffusion discretization schemes was discussed, with the extension
to higher order schemes in [29] and in [6] a method to discretize a stochastic
process via quantization optimally was given. Another randomization attempt
was presented in [10], where the volatility-of-volatility parameter in the
Heston model was assumed to be stochastic, and approximations, based on
perturbation, were derived. In the case of a randomizer with a finite number
of states, one may think of a regime-switching model [13]. Such a setup was
discussed [28], where the hidden regime continuous-time Markov chain was
considered.

Customarily, an extension with a stochastic parameter would require that the
model satisfies the linearity conditions- this, however, does not need to be
the case under the framework proposed here. Instead, we only require the
affinity to hold for a particular realization of a stochastic parameter. We
build an external layer over the class of AD models and allow the parameter to
be stochastic in that layer. In a standard approach, this would imply an
infinite number of affine models, making this method impractical. The solution
to this problem is at the heart of this article, in which we develop a method
for constructing an appropriate ChF based on only a few, critical parameter
realizations. Selecting these points is based on a quadrature rule based on
the moments of the stochastic parameter. The procedure of building a
randomized ChF, we call the RAnD method, which corresponds to “Randomized
Affine-Diffusion”.

The presented structure is generic, and we believe it may be applied in
various asset classes. To highlight its flexibility in modern finance, we
apply this technique to the challenging problem of consistent pricing options
on S&P 500 and VIX. This pricing problem is not new, and many innovative
solutions have already been proposed.

The major challenge regarding this pricing task is to find a stochastic model
that can sufficiently explain both the negative skew for equity options on the
S&P 500 index and the positive skew for options on the VIX. In addition,
different regime settings for implied volatility cause a calibration dilemma
where the volatility of the volatility parameter under the Heston model is
either too large for VIX or too small for S&P 500 (or both). The first
attempts to address this pricing problem aimed to introduce additional degrees
of freedom in multidimensional models. In [30], for example, the variance
process in the Bates model was enriched with an additional jump process. In
[1], the so-called 3/2 process and inclusion of jumps for the stock process
could capture the observed upward skew for VIX, although the model was unable
to calibrate to both derivatives simultaneously. More recently, in [23],
another attempt to include jumps in both returns and volatility added to
stochastic volatility was discussed. The authors concluded that more
flexibility, especially with a relaxed Feller condition, is necessary to
provide a sufficient fit. A promising extension through the Hawkes process for
the jumps was discussed [22]\- however, additional degrees of freedom are
insufficient for calibration purposes. A breakthrough in volatility modelling
came with the so-called “rough volatility” models. In [21], the asymptotic
expansions of the rough Bergomi [3] model to jointly describe the VIX and the
S&P 500 were presented. It was further extended in [18] with the application
of modulated Volterra processes. Recently, in [11], Gatheral et al. have
introduced and successfully calibrated the so-called quadratic rough Heston.
Another angle to address this intriguing pricing problem was proposed in [16],
where a discrete-time dispersion-constrained martingale transport problem was
applied.

The references enumerated above unquestionably address the problem of
consistent options pricing on an index and its corresponding volatility.
However, in our view, this is insufficient to claim victory- many of the
proposed solutions only work for short maturities or lack effective numerical
methods to be perceived as practical. We will show that by randomizing AD
models, one can benefit from well-established numerical techniques and the
ability to price long-term options while consistently pricing volatility
products.

The article is organized as follows: in Section 2, we start with the
definition of the class of AD and give an integral form for the randomized
ChF. In Section 2.1, the RAnD method is introduced, and extensions for
bivariate and multivariate parameters cases are presented. Section 2.3
presents an illustrative example of the randomization of the Black-Scholes
model. We analyze here the computational aspects of the method and the
evolution of implied volatilities. Section 3 focuses on the randomized Bates
model; here, we start with the analysis of randomization and its impact on the
implied volatilities; also, in this section, we address what was mentioned
earlier in the introduction- the randomization of jump sizes. Further, in
Section 3.1, we provide the pricing details for VIX options, and in Section
3.2, we use derived formulae to perform simultaneous model calibration to S&P
and VIX options. Discussion about hedging under the RAnD framework is covered
in Section 3.3 and in Sections 3.4 and 3.5, the details regarding Monte Carlo
simulation with possible improvements are given. Section 4 concludes.

##  2 Affine (Jump) Diffusion processes and randomization

The class of stochastic AD processes for the asset dynamics refers to a fixed
probability space $\left(\Omega,\mathcal{F}(t),\mathbb{Q}\right)$ and a
Markovian $n$-dimensional affine process
${\bf{X}}(t)=[X_{1}(t),\dots,X_{n}(t)]^{\rm T}$ in some space
$\mathbb{R}\subset\mathbb{R}^{n}$.

The stochastic model of interest can be expressed by the following stochastic
differential form:

${\rm d}{\bf{X}}(t)=\bar{\boldsymbol{\mu}}(t,{\bf{X}}(t)){\rm
d}t+\bar{\boldsymbol{\sigma}}(t,{\bf{X}}(t)){\rm
d}{\widetilde{\bf{W}}}(t)+{\bf J}(t)^{\rm T}{\rm d}{\bf X}_{\mathcal{P}}(t),$
(2.1)

where ${\widetilde{\bf{W}}}(t)$ is an $\mathcal{F}(t)$-standard column vector
of independent Brownian motions in $\mathbb{R}^{n}$,
$\bar{\boldsymbol{\mu}}{\left(t,{\bf{X}}(t)\right)}:\mathbb{R}\rightarrow\mathbb{R}^{n}$,
$\bar{\boldsymbol{\sigma}}\left(t,{\bf{X}}(t)\right):\mathbb{R}\rightarrow\mathbb{R}^{n\times
n}$, and ${\bf X}_{\mathcal{P}}(t)\in\mathbb{R}^{n}$ is a vector of orthogonal
Poisson processes, characterized by an intensity vector
$\bar{\boldsymbol{\xi}}(t,{\bf X}(t))\in\mathbb{R}^{n}$.

${\bf J}(t)\in\mathbb{R}^{n}$ is a vector governing the amplitudes of the
jumps and is assumed to be a matrix of correlated random variables, that are
independent of the state vector ${\bf X}(t)$ and of ${\bf
X}_{\mathcal{P}}$(t).

Additionally, we consider an orthogonal vector
$\Theta=[{\vartheta}_{1},\dots,{\vartheta}_{n}]^{\rm T}$, $n\in\mathbb{N}$,
where each ${\vartheta}_{i}$ is an independent, time-invariant, random
variable 111We consider here $n\in\mathbb{N}$ stochastic parameters, this is
however not a necessary constraint.. A realization of ${\vartheta_{i}}$ we
indicate by $\theta_{i}$, $\vartheta(\omega)=\theta$, per se the realization
for ${\Theta}$ is indicated by $\theta=[\theta_{1},\dots,\theta_{n}]^{\rm T}$.

For processes in the AD class, it is required that for a particular parameter
realization, $\theta$, drift $\bar{\boldsymbol{\mu}}(t,{\bf{X}}(t))$, interest
rate component $\bar{\boldsymbol{r}}(t,{{\bf{X}}(t)})$, and the covariance
matrix $\bar{\boldsymbol{\sigma}}(t,{\bf
X}(t))\bar{\boldsymbol{\sigma}}(t,{\bf X}(t))^{\rm T}$ are of the affine form,
i.e.

$\displaystyle\bar{\boldsymbol{\mu}}(t,{\bf{X}}(t))$ $\displaystyle=$
$\displaystyle a_{0}(\theta)+a_{1}(\theta){\bf{X}}(t),\text{ for any
}(a_{0}(\theta),a_{1}(\theta))\in\mathbb{R}^{n}\times\mathbb{R}^{n\times n},$
(2.2) $\displaystyle\bar{\boldsymbol{r}}(t,{\bf{X}}(t))$ $\displaystyle=$
$\displaystyle r_{0}(\theta)+r_{1}(\theta)^{\rm T}{\bf{X}}(t),\text{ for
}(r_{0}(\theta),r_{1}(\theta))\in\mathbb{R}\times\mathbb{R}^{n},$ (2.3)
$\displaystyle(\bar{\boldsymbol{\sigma}}(t,{\bf{X}}(t))\bar{\boldsymbol{\sigma}}(t,{\bf{X}}(t))^{\rm
T})_{i,j}$ $\displaystyle=$
$\displaystyle(c_{0}(\theta))_{i,j}+(c_{1}(\theta))^{\rm
T}_{i,j}{\bf{X}}_{j}(t),$ (2.4) $\displaystyle\bar{\xi}(t,{\bf X}(t))$
$\displaystyle=$ $\displaystyle l_{0}(\theta)+l_{1}(\theta){\bf X}(t),\text{
with }(l_{0}(\theta),l_{1}(\theta))\in\mathbb{R}^{n}\times\mathbb{R}^{n},$
(2.5)

with $(c_{0}(\theta),c_{1}(\theta))\in\mathbb{R}^{n\times
n}\times\mathbb{R}^{n\times n\times n}$, and where in (2.4) it is meant that
each element in the matrix
$\bar{\boldsymbol{\sigma}}(t,{\bf{X}}(t))\bar{\boldsymbol{\sigma}}(t,{\bf{X}}(t))^{\rm
T}$ should be affine, as well as each vector element in the drift and interest
rate vectors.

###  2.1 The RAnD method: details on construction of the characteristic
function

For a given realization of $\Theta$, $\theta$, we consider ${\bf
X}_{\theta}(t):={\bf X}(t)|{\bf\Theta}=\theta$, ${\bf J}_{\theta}(t):={\bf
J}(t)|{\bf\Theta}=\theta$. It is shown [8] that in this class, for a state
vector ${\bf{X}}_{\theta}(t)$, the discounted characteristic function is also
of the following form:

$\phi_{{\bf X}_{\theta}}({\bf{u}};t,T)=\mathbb{E}_{t}\left[{\rm
e}^{-\int_{t}^{T}r(s){\rm d}s+i{\bf u}^{\rm
T}{\bf{X}}_{\theta}(T)}\right]={\rm e}^{\bar{A}({\bf
u};\tau,\theta)+{\bf{\bar{B}^{\rm T}}}({\bf
u};\tau,\theta){\bf{X}}_{\theta}(t)},$ (2.6)

with the expectation under risk-neutral measure $\mathbb{Q}$ for $\tau=T-t.$

The coefficients $\bar{A}:=\bar{A}({\bf u};\tau,\theta)$ and $\bar{\bf
B}:={\bf{\bar{B}}}^{\rm T}({\bf u};\tau,\theta)$ have to satisfy the following
complex-valued Riccati ODEs, see the work by Duffie-Pan-Singleton [8]:

$\displaystyle\frac{{\rm d}\bar{A}}{{\rm d}\tau}$
$\displaystyle=-r_{0}(\theta)+{\bf{\bar{B}}}^{\rm
T}a_{0}(\theta)+\frac{1}{2}{\bf{\bar{B}}}^{\rm
T}c_{0}(\theta){\bf{\bar{B}}}+l_{0}^{\rm T}\mathbb{E}\left[{\rm e}^{{\bf
J}_{\theta}(\tau){\bf\bar{B}}}-1\right],$ (2.7) $\displaystyle\frac{{\rm
d}\bf{\bar{B}}}{{\rm d}\tau}$ $\displaystyle=-r_{1}(\theta)+a_{1}(\theta)^{\rm
T}{\bf{\bar{B}}}+\frac{1}{2}{\bf{\bar{B}}}^{\rm
T}c_{1}(\theta){\bf{\bar{B}}}+l_{1}(\theta)^{\rm T}\mathbb{E}\left[{\rm
e}^{{\bf J}_{\theta}(\tau){\bf\bar{B}}}-1\right],$

where the expectation, $\mathbb{E}[\cdot]$ in (2.7), is taken with respect to
the jump amplitude ${\bf J}_{\theta}(t)$.

The dimension of the (complex-valued) ODEs for ${\bf{\bar{B}}}({\bf
u};\tau,\theta)$ corresponds to the dimension of the state vector
${\bf{X}}(t)$. Then, for stochastic parameter ${\bf\vartheta}$, the ChF is
given by:

$\displaystyle\phi_{{\bf X}}({\bf u};t,T):=\mathbb{E}_{t}[{\rm
e}^{-\int_{t}^{T}r(s){\rm d}s+i{\bf u}^{\rm T}{\bf
X}(T)}]=\mathbb{E}_{t}\left[\mathbb{E}_{t}\big{[}{\rm
e}^{-\int_{t}^{T}r(s){\rm d}s+i{\bf u}^{\rm T}{\bf
X}_{\theta}(T)}\big{|}{\bf\Theta}={\theta}\big{]}\right].$

The inner expectation can be recognized as the conditional ChF in (2.6); thus,
by definition of the ChF and integration over all the parameter space, we
find,

$\displaystyle\phi_{{\bf X}}({\bf u};t,T)=\mathbb{E}_{t}\left[\phi_{{\bf
X}|{\Theta}}({\bf u};t,T)\right]=\int_{\mathbb{R}^{n}}\phi_{{\bf
X}|{\Theta}=\theta}({\bf u};t,T)f_{\Theta}(\theta){\rm d}{\theta}.$ (2.8)

We aim to provide numerically efficient methods for computation of ChF in
(2.8).

######  Remark 1 (Simplified notation).

The representation above is rather generic in terms of multidimensional ChF
and possibly a multidimensional set of stochastic parameters. However, such a
generic case is rather exceptional, as, in a typical pricing application, a
payoff function will only depend on a single, marginal distribution (even so,
possibly, driven by the multidimensional system of the SDEs); therefore, for
the sake of simplicity, from now on, we focus on the derivations of the Chf
for ${\bf u}=[u,0,\dots]$, and we denote it by $u$. Moreover, although the
presented framework allows for simultaneous randomization of multiple
parameters, we will consider a single stochastic parameter $\vartheta$ with
$\theta$ indicating its realization. In the following sections, we consider a
finite number of realizations of $\vartheta$; therefore, we will denote them
as $\theta_{1},\dots,\theta_{N}$, for some $N\in\mathbb{N}$. These “specific”
realizations we will interchangeably call either “collocation” [15] or
“quadrature” points. Finally, by $\vartheta(\hat{a},\hat{b},\cdot)$ we denote
that $\theta$ is a random variable driven by parameters, $\hat{a}$, $\hat{b}$.

To determine the ChF of an affine model with a randomized parameter, one needs
to integrate the parameter’s probability density function over the whole
domain (2.8). This can be avoided, i.e., the complicated integrand can be
factored into a set of pairs
$\\{\\{\omega_{1},\theta_{1}\\},\dots,\\{\omega_{N},\theta_{N}\\}\\}$,
$N\in\mathbb{N}$, with a nonnegative “weights” function, $\omega_{i}\geq 0$,
such that $\sum_{n=1}^{N}\omega_{n}=1$ and specific, collocation, points
$\theta_{n}$. Once the number of evaluations, $N$, is low, we can
significantly reduce the computational cost. The key element here, however, is
that the pairs, $\\{\omega_{n},\theta_{n}\\}_{n=1}^{N}$, cannot be chosen
arbitrarily but need to be computed based on the parameter’s distribution,
$\vartheta$. We consider a random parameter $\vartheta$ with its PDF,
$f_{\vartheta}(\cdot)$, such that for a fixed number $N\in\mathbb{N}$, moments
are finite, i.e., $\mathbb{E}[\vartheta^{2N}]<\infty$. Then we are able to
determine a mapping, often called quadrature rule, function
$\zeta(\vartheta):\mathbb{R}\rightarrow\\{\omega_{n},\theta_{n}\\}_{n=1}^{N}$.
We follow the approach presented by G.H. Golub and J.H. Welsch in [12] where
$\omega_{n}$ are the quadrature weights determined based on the moments of the
random parameter, $\vartheta$. As explained in [12], one can calculate pairs
$\\{\omega_{n},\theta_{n}\\}_{n=1}^{N}$ based on a three-term recurrence
relation which is well-known for orthogonal polynomials generated by
$f_{\vartheta}(x)$, or equivalently, by the moments of $\vartheta$. Although
the technique to compute these pairs of points is well-founded, for the
readers’ convenience, the theoretical aspects and the algorithmic details
222The accompanying Python and MATLAB codes are available at
https://github.com/LechGrzelak/Randomization are included in A.1. Throughout
the article, the method to compute these points will be denoted as
$\zeta(\vartheta):\mathbb{R}\rightarrow\\{\omega_{n},\theta_{n}\\}_{n=1}^{N}$.

The RAnD method for determining the ChF with a randomized model parameter is
presented in the theorem below.

######  Theorem 2.1 (Characteristic Function for Randomized Affine Jump
Diffusion Processes).

Consider a random variable $\vartheta$ defined on some finite domain
$D_{\vartheta}:=[a,b]$, with its PDF, $f_{\vartheta}(x)$, CDF,
$F_{\vartheta}(x)$ and a realization $\theta$, $\vartheta(\omega)=\theta$ such
that for some $N\in\mathbb{N}$ the moments are finite,
$\mathbb{E}[\vartheta^{2N}]<\infty$. Let ${\bf X}(t)$ represent an affine
state vector with some constant parameter $\theta$. Assuming that the
corresponding ChF, $\phi_{{\bf X}|\vartheta=\theta}(\cdot)$, is well defined
and $2N$ times differentiable w.r.t $\theta$, the unconditional ChF for the
randomized ${\bf X}$, exists and is given by:

$\displaystyle\phi_{{\bf X}}({u};t,T)=\sum_{n=1}^{N}\omega_{n}\phi_{{\bf
X}|\vartheta=\theta_{n}}({u};t,T)+\epsilon_{N}=\sum_{n=1}^{N}\omega_{n}{\rm
e}^{\bar{A}({u};\tau,\theta_{n})+{\bf{\bar{B}^{\rm
T}}}({u};\tau,\theta_{n}){\bf{X}}(t)}+\epsilon_{N},$ (2.9)

with $\bar{A}(u;\tau,\theta_{n})$, $\bar{{\bf B}}^{\rm T}(u;\tau,\theta_{n})$
are defined in (2.7) and where

$\epsilon_{N}=\frac{1}{(2N)!}\frac{\partial^{2N}}{\partial\xi^{2N}}\phi_{{\bf
X}|\vartheta=\xi}({u};t,T),$ (2.10)

for $a<\xi<b$ and where the pairs $\\{\omega_{n},\theta_{n}\\}_{n=1}^{N}$ are
the Gauss-quadrature weights and the nodes based on the parameter
distribution, $f_{\vartheta}(\cdot),$ determined by
$\zeta(\vartheta):\mathbb{R}\rightarrow\\{\omega_{n},\theta_{n}\\}_{n=1}^{N}$
defined in A.1.

###### Proof.

Starting with the definition of the ChF and conditioning on a parameter
$\vartheta,$ we find:

$\displaystyle\phi_{{\bf X}}({\bf u};t,T):=\mathbb{E}_{t}[{\rm
e}^{-\int_{t}^{T}r(s){\rm d}s+i{\bf u}^{\rm T}{\bf
X}(T)}]=\mathbb{E}_{t}\left[\mathbb{E}_{t}\big{[}{\rm
e}^{-\int_{t}^{T}r(s){\rm d}s+i{\bf u}^{\rm T}{\bf
X}_{\theta}(T)}\big{|}{\vartheta}={\theta}\big{]}\right].$

Then, by definition of the ChF,

$\displaystyle\phi_{{\bf X}}({\bf u};t,T)$ $\displaystyle=$
$\displaystyle\mathbb{E}_{t}\left[\phi_{{\bf X}|\vartheta=\theta}({\bf
u};t,T)\right]=\int_{\mathbb{R}}\phi_{{\bf X}|\vartheta=\theta}({\bf
u};t,T)f_{\vartheta}(x){\rm d}x$ $\displaystyle=$
$\displaystyle\sum_{n=1}^{N}\omega_{n}\phi_{{\bf X}|\vartheta=\theta_{n}}({\bf
u};t,T)+\epsilon_{N}.$

First, we prove that a convex linear combination of characteristic functions,
$\sum_{n=1}^{N}\omega_{n}\phi_{{\bf X}|\vartheta=\theta_{n}}({\bf u};t,T)$,
and $\sum_{n=1}^{N}\omega_{n}=1$ with $\omega_{n}>0$ of a finite or a
countable number of characteristic function is the characteristic function.

Since $\omega_{n}$ can be associated with probabilities, we can define a
discrete random variable $\mathbb{P}[\vartheta=\theta_{n}]=\omega_{n}$, so we
have:

$\displaystyle\sum_{n=1}^{N}\omega_{n}\phi_{{\bf X}|\vartheta=\theta_{n}}({\bf
u};t,T)=\sum_{n=1}^{N}\mathbb{E}[1_{\vartheta=\theta_{n}}]\mathbb{E}_{t}\big{[}{\rm
e}^{i{\bf u}^{\rm T}{\bf
X}_{\theta_{n}}(T)}\big{]}=\sum_{n=1}^{N}\mathbb{E}\big{[}1_{\vartheta=\theta_{n}}{\rm
e}^{i{\bf u}^{\rm T}{\bf X}_{\theta_{n}}(T)}\big{]},$

since $\vartheta$ is independent of ${\bf X}(t)$. For Fubini’s theorem to
hold, the boundness condition needs to be satisfied,

$\displaystyle\sum_{n=1}^{N}\mathbb{E}\big{[}\big{|}1_{\vartheta=\theta_{n}}{\rm
e}^{i{\bf u}^{\rm T}{\bf
X}_{\theta_{n}}(T)}\big{|}\big{]}=\sum_{n=1}^{N}\mathbb{E}\big{[}1_{\vartheta=\theta_{n}}\big{|}{\rm
e}^{i{\bf u}^{\rm T}{\bf
X}_{\theta_{n}}(T)}\big{|}\big{]}\leq\sum_{n=1}^{N}\mathbb{P}\big{[}\vartheta=\theta_{n}\big{]}=1<\infty,$

therefore:

$\displaystyle\sum_{n=1}^{N}\omega_{n}\phi_{{\bf X}|\vartheta=\theta_{n}}({\bf
u};t,T)=\mathbb{E}\big{[}\sum_{n=1}^{N}1_{\vartheta=\theta_{n}}{\rm e}^{i{\bf
u}^{\rm T}{\bf X}_{\theta_{n}}(T)}\big{]}=\mathbb{E}\big{[}\exp\big{(}{i{\bf
u}^{\rm T}\sum_{n=1}^{N}1_{\vartheta=\theta_{n}}{\bf
X}_{\theta_{n}}(T)}\big{)}=:\phi_{Y}({\bf u}),$

with $\phi_{Y}(\cdot)$ being a ChF of
$\sum_{n=1}^{N}1_{\vartheta=\theta_{n}}{\bf X}_{\theta_{n}}(T)$.

Using the fact that the collocation points $\theta_{n}$ for $n=1,\dots,N$ are
optimal and correspond to the zeros of orthogonal polynomials, we relate the
collocation method to Gaussian quadrature where the integral of a complex-
valued function, $\phi_{{\bf X}|\vartheta=\theta}({\bf u};t,T)$, can be
approximated by a polynomial:

$\int_{\mathbb{R}}\phi_{{\bf X}|\vartheta=x}({\bf u};t,T)f_{\vartheta}(x){\rm
d}x-\sum_{n=1}^{N}\phi_{{\bf X}|\vartheta=\theta_{n}}({\bf
u};t,T)\omega_{i}=\frac{1}{(2N)!}\frac{\partial^{2N}}{\partial\xi^{2N}}\phi_{{\bf
X}|p=\xi}({\bf u};t,T)=:\epsilon_{N},$

with $f_{\vartheta}(x)$ the weight function, and $\omega_{n}$ the quadrature
weights. The error is given by $\epsilon_{N}$, for $a<\xi<b$, and the pair
$\\{\omega_{n},\theta_{n}\\}_{n=1}^{N}$ are the Gauss-quadrature weights and
the nodes based on the parameter distribution as proven in [5] (p.180, Theorem
3.6.24). ∎

Theorem 2.1 illustrates that the ChF of the randomized AD model is a weighted
sum of a set of conditional ChFs evaluated at certain realizations,
$\theta_{n}$, of the underlying stochastic parameter $\vartheta$. The theorem
also shows, in Equation (2.10), the exponential decay of the error in terms of
$N$ associated with the quadrature approximation- suggesting that only a few
terms will be needed to reach high-precision.

The key feature of the RAnD method is its ability to facilitate fast
evaluation of the randomized ChF (2.9). The construction of the ChF is based
on the randomized variable, $\vartheta$, via its $2N$ moments, i.e., the
availability of the moments allows us to determine $\zeta(\vartheta)$ which
provides the quadrature points $\\{\omega_{n},\theta_{n}\\}_{n=1}^{N}$. It is,
therefore, essential that the randomizing variables enable the closed-form
moment calculation. Moreover, when specifying a randomizing variable, its
support must correspond to the physical range of the randomized parameter; for
example,– some of the model parameters may be only valid for specific ranges,
such as correlation or volatility coefficients. In Table 1, a few possible
distributions with their corresponding moments are provided. When analytical
moments are available, the computation of the corresponding points only
requires the computation of a Cholesky decomposition and certain eigenvalues
(see A.1); it is therefore computationally cheap.

Table 1: Selected distributions for the stochastic parameters. For the normal random variable with some mean, $\mu$, and variance, $\sigma^{2}$, it is sufficient to consider standard normal distribution, $\mathcal{N}(0,1)$, and properly scale the $\theta_{n}$ points, obtained from Algorithm 1. name | raw moment | domain
---|---|---
$\vartheta\sim\mathcal{U}([\hat{a},\hat{b}])$ | $\mathbb{E}[\vartheta^{n}]=\frac{\hat{b}^{n+1}-\hat{a}^{n+1}}{(n+1)(\hat{b}-\hat{a})}$ | $[\hat{a},\hat{b}]$
$\vartheta\sim\exp(\hat{a})$ | $\mathbb{E}[\vartheta^{n}]=\frac{n!}{\hat{a}^{n}}$ | $\mathbb{R}^{+}$
$\vartheta\sim\mathcal{N}(0,1)$ |  $\mathbb{E}[\vartheta^{n}]=(n-1)!!$ if $n$ even; $0$ otherwise | $\mathbb{R}$
$\vartheta\sim\Gamma(\hat{a},\hat{b})$ | $\mathbb{E}[\vartheta^{n}]=\hat{b}^{n}\Gamma(n+\hat{a})/\Gamma(\hat{b})$ | $\mathbb{R}^{+}$
$\vartheta\sim\chi^{2}(\hat{a},\hat{b})$ | $\mathbb{E}[\vartheta^{n}]=2^{n-1}(n-1)!(\hat{a}+n\hat{b})+\sum_{j=1}^{n-1}\frac{(n-1)!2^{j-1}}{(n-j)!}(\hat{a}+j\hat{b})\mathbb{E}[\vartheta^{n-j}]$ | $\mathbb{R}^{+}\cup\\{0\\}$

###  2.2 Pricing with the RAnD method

We will focus on pricing two different plain vanilla options, namely European
options on an asset, $S(t)$, and options on the associated volatility, VIX.
The pricing will rely on a Fourier inversion method, namely the COS method
[26], which, by a change of the pricing coefficients, $H_{k}$ and the
corresponding ChF allows for pricing of both contracts. The generic pricing
equation is given by:

$V(t_{0})={\rm
e}^{-r(T-t_{0})}\sideset{}{{}^{\prime}}{\sum}_{k=0}^{N_{c}-1}\Re\left[\phi_{\bf
X}\left(\frac{k\pi}{b-a};t_{0},T\right)\exp\left(-ik\pi\frac{a}{b-a}\right)\right]\cdot
H_{k}+\epsilon_{c_{1}},$ (2.11)

where $N_{c}$ represents the number of expansion terms, $\Re(\cdot)$ is the
real part and where $H_{k}$ for $k\geq 0$ are known in closed-form
coefficients corresponding to the payoff function. We will use the $H_{k}$
coefficients derived for European-style call/put options (see [26] for
details) and for options on VIX (presented later, in Section 3, Lemma 3.2).
Parameters, $a$ and $b$ are the tuning parameters used to determine the
integration range; while the error $\epsilon_{c_{1}}$, is exponentially
decaying in $N_{c}$.

######  Remark 2 (Inversion of the randomized ChF).

The application of Fourier inversion to the randomized ChF in (2.9) yields,

$\displaystyle f_{{\bf X}}(x)=\frac{1}{2\pi}\int_{\mathbb{R}}{\rm
e}^{-iux}\sum_{n=1}^{N}\omega_{n}\phi_{{\bf X}|\vartheta=\theta}({u};t,T){\rm
d}u=\sum_{n=1}^{N}\omega_{n}f_{{\bf X}|\vartheta=\theta}(x).$ (2.12)

Since $\omega_{1}+\dots+\omega_{N}=1$, $\omega_{n}\geq 0$, for $n=1,\dots,N$,
which implies that the density of the affine, randomized, system of SDEs,
${\bf X}(t)$ can be expressed as a, possibly multi-modal, mixture
distribution. Although the exact representation for the density $f_{\bf
X}(\cdot)$ may be challenging to derive analytically, one can use Fourier
inversion to recover the corresponding density and utilize it for swift option
pricing.

A natural extension of the RAnD method is to consider two parameters that are
stochastic and follow a bivariate distribution. Such an extension would
require that the conditional moments are known explicitly- which is somewhat
troublesome. To this day, only a few distributions allow for exact moments.
However, if we stay, for example, within the Gaussian world, such an extension
to a 2D case is trivial. Corallary 2.1 gives us the general details.

######  Corollary 2.1 (Random parameters with bivariate distribution).

Under a bivariate distribution $\Theta=[\vartheta_{1},\vartheta_{2}]$ with
$\zeta(\vartheta_{1})=\\{\omega_{1,n},\theta_{1,n}\\}_{n=1}^{N}$ and
conditioned on
$\zeta(\vartheta_{2}|\vartheta_{1})=\\{\omega_{2,n},\theta_{2,n}\\}_{n=1}^{N}$,
the randomized ChF is given by:

$\displaystyle\phi_{{\bf
X}}({u};t,T)=\sum_{n_{1}=1}^{N}\omega_{n_{1}}\sum_{n_{2}=1}^{N}\omega_{n_{2}}\phi_{{\bf
X}|\vartheta_{1}=\theta_{n_{1}},\vartheta_{2}=\theta_{n_{2}}}({u};t,T)+\tilde{\epsilon}_{N},$
(2.13)

where $N$ indicates the number of expansion terms,
$\vartheta_{2}|\vartheta_{1}$ indicates a conditional random variable, and the
remaining specification follows Theorem 2.1.

###  2.3 Illustrative case: the RAnD Black-Scholes model

As the first example for the RAnD method, let us consider the randomized
Black-Scholes (RAnD BS) model. Given the probability space
$(\Omega,\mathcal{F}(t),\mathbb{Q})$, we define a stock
$S(t):=\\{S(t,\omega):t\in T\\}$ and the randomizing volatility
$\sigma:=\\{\sigma(\omega):\Omega\rightarrow\mathbb{R}^{+}\\}$, with some
constant interest rate, $r$. We consider three different, continuous,
randomizers for $\sigma$, i.e., we assume the volatility parameter, $\sigma$,
to either follow a uniform distribution,
$\sigma\sim\mathcal{U}([\cdot,\cdot])$, gamma distribution,
$\sigma\sim\Gamma(\cdot,\cdot)$, or the non-central chi-squared distribution,
$\sigma\sim\chi^{2}\left(\cdot,\cdot\right)$; thus, also mimicking a similar
structure as in the Heston model where the CIR dynamics drive the variance
process (see Equations (3.3) and (3.4) for the exact relation).

Formally, the RAnD BS model is given by the following dynamics:

$\displaystyle{\rm d}S(t)$ $\displaystyle=$ $\displaystyle rS(t){\rm
d}t+\sigma S(t){\rm d}W(t),\;\;S(t_{0})>0,$ $\displaystyle\sigma$
$\displaystyle\sim$
$\displaystyle\mathcal{U}([\hat{a},\hat{b}]),\;\;\;\;{\text{or}}\;\;\;\;\sigma\sim\Gamma(\hat{a},\hat{b}),\;\;\;\;\text{or}\;\;\;\;\sigma\sim\hat{a}\chi^{2}(\hat{b},\hat{c}),$

for $\hat{a},\hat{b},\hat{c}\in\mathbb{R}^{+}.$ As it will be presented later,
the randomization of the BS model is elegant as it allows for a generation of
the European option prices in closed form (see [4] where discrete randomizers
for the BS model were covered). However, since the randomized volatility is
time-invariant, or stationary, such a model cannot compete with, e.g. the
model of Heston, where the volatility has a term structure, and it is driven
by a correlated stochastic process. Contrary to the Heston model, however, the
RAnD BS model does not require Fourier transformation to price options.

We have two techniques to compute option prices of the randomized Black-
Scholes model. First, a straightforward approach is based on the closed-form
solution of the “unrandomized” model. The second method is based on
randomization at the ChF level. By application of the RAnD method and
conditioning on the volatility realizations, the price of a European call
option, for $t_{0}=0$, reads:

$\displaystyle V(t_{0})$ $\displaystyle=$ $\displaystyle{\rm
e}^{-rT}\mathbb{E}_{t_{0}}\left[\max\left(S(T)-K,0\right)\right]$
$\displaystyle=$ $\displaystyle{\rm
e}^{-rT}\mathbb{E}_{t_{0}}\left[\mathbb{E}_{t_{0}}\left[\max\left(S(T)-K,0\right)\big{|}\sigma=\sigma_{n}\right]\right]={\rm
e}^{-rT}\mathbb{E}_{t_{0}}\left[V_{c}(t_{0},T,K,S(t_{0}),r,\sigma_{n})\right],$

with $\mathbb{E}_{t_{0}}[\cdot]$ being the expectation operator with
filtration $\mathcal{F}(t_{0})$, where
$V_{c}(t_{0},T,K,S(t_{0}),r,\sigma_{n})$ represents a European call option
price with strike $K$, expiry $T$ and the corresponding interest rate $r$, and
it is defined in (2.15) below. The price $V(t_{0})$ is now given by:

$\displaystyle V(t_{0})$ $\displaystyle=$ $\displaystyle{\rm
e}^{-rT}\int_{\mathbb{R}^{+}}V_{c}(t_{0},T,K,S(t_{0}),r,x)f_{\sigma}(x){\rm
d}x$ (2.14) $\displaystyle=$ $\displaystyle{\rm
e}^{-rT}\sum_{n=1}^{N}\omega_{n}V_{c}(t_{0},T,K,S(t_{0}),r,\sigma_{n})+\hat{\epsilon}_{N},$

where pairs $\\{\omega_{n},\sigma_{n}\\}$ are the weights and the quadrature
points, determined based on the distribution of $\sigma(\omega)$ and are
computed using Algorithm 1 in A.1, $\hat{\epsilon}_{N}$ is the quadrature
error defined in (2.9).

At each realization $\sigma_{n},$ the pricing is known in closed form and is
given the Black-Scholes formula:

$\displaystyle
V_{c}(t_{0},T,K,S(t_{0}),r,\sigma_{n})=S(t)F_{\mathcal{N}(0,1)}(d_{1})-KF_{\mathcal{N}(0,1)}(d_{2}){\rm
e}^{-r(T-t_{0})},$ (2.15)

with $d_{1}=\frac{1}{\sqrt{\sigma_{n}^{2}(T-t_{0})}}\left[\log
S(t_{0})/K+\left(r+1/2\sigma_{n}^{2}\right)(T-t_{0})\right],\;\;\;d_{2}=d_{1}-\sqrt{\sigma_{n}^{2}(T-t_{0})}.$

######  Remark 3 (Continuous randomizers for the BS model- previous attempts).

Early attempts to randomize the Black-Scholes model using continuous random
variables are already known in the literature. In [19], the authors have
randomized the variance of the Black-Scholes model by a CEV-generated
distribution. The primary aim of such randomization was to achieve the so-
called moderately explosive distribution allowing for enhanced skew. The
presented framework, however, although via the asymptotic expansion allowed
for short-maturity proxies for the implied volatilities, is very much problem
and randomizer specific. Moreover, due to the nature of the expansion, it
performs well only for short-expiry derivatives, contrary to the RAnD method
discussed in this article.

####  2.3.1 Randomization of the ChF under the Black-Scholes model

The presented randomization technique is not limited to the Black-Scholes
model, where a closed-form solution exists for each volatility realization,
$\sigma_{n}$. In order to achieve more flexibility and be able to price
challenging derivatives, the stochastic volatility models, like the Heston and
Bates models, will be considered. A numerical routine, like Fourier inversion,
needs to be used. Although such inversions are swift due to exponentially
convergent methods, like the COS method [26], it is important to keep the
number of inversions as low as possible. This requirement suggests that a
preferred way to compute a randomized option price is via randomization at the
ChF level, introduced in (2.9), which for the RAnD BS model is given by the
following expression,

$\displaystyle\phi_{{X}}({u};t_{0},T)=\sum_{n=1}^{N}\omega_{n}\exp\big{(}{iuX(t_{0})+(r-\frac{1}{2}\sigma_{n}^{2})iu(T-t_{0})-\frac{1}{2}\sigma_{n}^{2}u^{2}(T-t_{0})}\big{)}+\epsilon_{N},$
(2.16)

with $X(t_{0})=\log S(t_{0})$, the error $\epsilon_{N}$ defined in (2.9).

The pricing of European-style options can be performed using the expansion
formula in (2.11).

At this point, one can still employ the pricing equation (2.11) to (2.14);
however, it would require $N$ evaluations of the pricing equation in (2.11).
Since the strength of the COS method is its ability to compute option prices
for a range of strikes $K$, it is important to keep the number of vector-
matrix multiplications in (2.11) to a minimum; it is, therefore, beneficial to
evaluate the pricing only once- thus the randomization should take place at
the ChF level (as a ChF does not depend on strikes). An additional element
that needs to be considered is the generated error, i.e., depending on the
smoothness of the option value (2.14) or the corresponding ChF (2.16), error
$\hat{\epsilon}_{N}$ may differ from the $\epsilon_{N}$ in (2.14).

####  2.3.2 RAnD BS model: convergence and implied volatilities

In order to assess the convergence of the RAnD method and the impact of the
randomization on implied volatilities, we consider a numerical experiment,
where two pricing approaches, presented in the previous section, will be
compared. In particular, we are interested in the implied volatility behaviour
for short maturity options ($T=0.1$), as this is often considered a
challenging case. We start by specifying the parameters for the randomizing
random variables (see Table 2).

Table 2: Parameters used in the RAnD BS model simulation. For the non-central chi-square distributions, the parameters were chosen based on the Heston model (see Dynamics (3.3)), with: $\gamma=1.9$, $\kappa=0.5$, $\bar{v}=0.3$, $v_{0}=0.3$, $T=0.1.$ | $\hat{a}$ | $\hat{b}$ | $\hat{c}$
---|---|---|---
$\sigma\sim\mathcal{U}([\hat{a},\hat{b}])$ | 0.1 | 0.45 | -
$\sigma\sim\Gamma(\hat{a},\hat{b})$ | 2.55 | 0.1 | -
$\sigma\sim\hat{a}\chi^{2}(\hat{b},\hat{c})$ | 0.088 | 0.1662 | 3.2417

With the distributions given, we are able to determine the corresponding pairs
$\\{\omega_{n},\sigma_{n}\\},$ for $n=1,\dots,N$. These points are based on
the moments of the randomizing random variable (see Table 1) and utilization
of Algorithm 1. Figure 1 (left) illustrates the method’s convergence regarding
the number of collocation points $N$. Surprisingly, for the randomized Black-
Scholes model, it does not matter whether the randomization takes place at the
price level or the corresponding ChF. However, the convergence rate is
affected by the associated randomizer. We report the highest convergence speed
for the uniform distribution, lower for gamma, and the slowest convergence
rate for the non-central chi-square randomizer. The phenomenon can be
explained by the amount of randomness introduced to the pricing equation.

Figure 1: Left: The RAnD method convergence for the number of collocation
points $N$. Red markers correspond to the results obtained based on analytical
expression (2.14), and blue lines correspond to Fourier expansion (2.16). The
“error” is the difference between option prices against the reference, fetched
with 5 million paths. The summation is over the entire range of option
strikes. Right: Implied volatility surface for the RAnD BS model for
${\sigma\sim\mathcal{U}([\hat{a},\hat{b}])}$.

In Figure 1 (right), the impact of $\sigma\sim\mathcal{U}([\hat{a},\hat{b}]),$
on the implied volatilities is illustrated (the remaining cases are
illustrated in A.2 in Figure 10). As expected, the RAnD BS model can generate
an implied volatility smile. However, because of the independence between the
stock’s Brownian motion and the randomizer, the skew control is limited,
irrespective of the randomizing distribution. It is important to note that by
introducing stochastic but time-invariant volatility, we can generate an
implied volatility surface that develops in time from a very peaked, almost
explosive, to a flattening shape. The implied volatility surface pattern is
very similar in all three cases. The RAnD BS model could be further extended,
with, for example, a displacement parameter that would introduce skewness in
the model, thus further extending the applicability of the Black-Scholes
model’s standard.

In the final experiment for the RAnD BS model, we report the convergence rate
in terms of the implied volatilities in Table 3. We see that the convergence
is not affected by the maturity time $T$. We expect that, with stochastic
volatility models, the convergence will depend on the time-to-expiry. The
results are excellent; for $N=6$ points, the maximum absolute error is about
$0.05-0.07\%$, which is a convincing result and implies that the method can be
used for accurate model calibration and pricing. As expected, the reported
errors depend on the randomizing variable (compare Table 3 against Table 5 in
A.5).

Table 3: Maximum implied volatility, $\text{IV}(\cdot)$, error computed against Monte Carlo results with 10 million samples, with $K=S_{0}{\rm e}^{0.1\sqrt{T}\delta}$, $\delta=[-3,-2,-1,-0.5,0,0.5,1,2,3]$. $\text{error}=\max_{i}{\big{|}\text{IV}(K_{i})-\text{IV}_{\text{ref}}(K_{i})\big{|}}$ |
---|---
$\sigma\sim\Gamma(a,b)$ | $N=2$ | $N=3$ | $N=4$ | $N=5$ | $N=6$ | $N=7$ | $N=8$ | $N=9$
$T=1d$ | $0.56\;\%$ | $0.26\;\%$ | $0.15\;\%$ | $0.10\;\%$ | $0.07\;\%$ | $0.06\;\%$ | $0.05\;\%$ | $0.04\;\%$
$T=1w$ | $0.56\;\%$ | $0.26\;\%$ | $0.15\;\%$ | $0.10\;\%$ | $0.07\;\%$ | $0.06\;\%$ | $0.04\;\%$ | $0.04\;\%$
$T=2w$ | $0.55\;\%$ | $0.25\;\%$ | $0.14\;\%$ | $0.09\;\%$ | $0.06\;\%$ | $0.04\;\%$ | $0.03\;\%$ | $0.03\;\%$
$T=1m$ | $0.54\;\%$ | $0.24\;\%$ | $0.13\;\%$ | $0.08\;\%$ | $0.06\;\%$ | $0.04\;\%$ | $0.03\;\%$ | $0.03\;\%$
$T=3m$ | $0.54\;\%$ | $0.24\;\%$ | $0.13\;\%$ | $0.08\;\%$ | $0.05\;\%$ | $0.04\;\%$ | $0.03\;\%$ | $0.02\;\%$
$T=6m$ | $0.55\;\%$ | $0.25\;\%$ | $0.14\;\%$ | $0.09\;\%$ | $0.07\;\%$ | $0.05\;\%$ | $0.04\;\%$ | $0.03\;\%$
$T=12m$ | $0.56\;\%$ | $0.26\;\%$ | $0.15\;\%$ | $0.10\;\%$ | $0.07\;\%$ | $0.05\;\%$ | $0.04\;\%$ | $0.04\;\%$

Further analysis regarding randomization and implied volatility shapes will be
discussed later in the context of the Heston and the Bates models.

##  3 The RAnD framework and pricing of plain vanilla and VIX options

This section focuses on applying the RAnD method to stochastic volatility and
jump models. In particular, we will consider the Bates [2] model, which forms
an elegant extension of the model of Heston [17]. From the pricing
perspective, the characteristic function of both models is known analytically.
Although some may argue that managing jumps in the pricing framework is of
limited value, especially for hedging, it is well known that pure stochastic
volatility, affine, models cannot generate sufficient skew for short maturity
options. Adding jumps to calibrate short-maturity options is often a preferred
choice. However, it is often insufficient. Our goal is to use the RAnD method
and apply it to the challenging task of simultaneous pricing of plain-vanilla
and VIX options. In particular, we illustrate the ability to fit both option
markets simultaneously by randomizing some Bates model parameters.

Randomization of individual Heston model parameters is known in the
literature. Initial random volatility was discussed in [24, 20, 25], and in
[10], the perturbation-based approximations for stochastic vol-of-vol
parameter were derived.

We will focus on the Bates model: the corresponding ChF, pricing of plain-
vanilla and VIX options and calibration, under the RAnD method, to market
data.

The Bates model, under the $\mathbb{Q}$ measure, is described by the following
system of SDEs:

$\displaystyle\begin{array}[]{l}{\rm
d}S(t)/S(t)=\left(r-\lambda\mathbb{E}[{\rm e}^{J}-1]\right){\rm
d}t+\sqrt{v(t)}{\rm d}W_{x}(t)+\left({\rm e}^{J}-1\right){\rm
d}X_{\mathcal{P}}(t),\\\\[6.45831pt] {\rm
d}v(t)=\kappa\left(\bar{v}-v(t)\right){\rm d}t+\gamma\sqrt{v(t)}{\rm
d}W_{v}(t),\end{array}$ (3.3)

with Poisson process $X_{\mathcal{P}}(t)$, intensity $\lambda$, and normally
distributed jump sizes, $J\sim\mathcal{N}(\mu_{j},\sigma^{2}_{j})$, with
$\mathbb{E}[{\rm e}^{J}]={\rm e}^{\mu_{J}+\frac{1}{2}\sigma_{J}^{2}}.$
$X_{\mathcal{P}}(t)$ is assumed to be independent of the Brownian motions and
the jump sizes. There is a correlation $\rho$ between the governing Brownian
motions, $\rho{\rm d}t={\rm d}W_{x}(t){\rm d}W_{v}(t)$. Under this model the
variance process follows the non-central chi-square distribution,
$\chi^{2}(\delta,\bar{\kappa}(\cdot,\cdot))$ with $\delta$ degrees of freedom
and non centrality parameter $\bar{\kappa}(t_{0},t)$,

$\displaystyle
v(t)|v(t_{0})\sim\bar{c}(t_{0},t)\chi^{2}(\delta,\bar{\kappa}(t_{0},t)),$
(3.4)

where

$\displaystyle\bar{c}(t_{0},t)=\frac{\gamma^{2}}{4\kappa}(1-{\rm
e}^{-\kappa(t-t_{0})}),\;\;\;\delta=\frac{4\kappa\bar{v}}{\gamma^{2}},\;\;\;\bar{\kappa}(t_{0},t)=\frac{4\kappa{\rm
e}^{-\kappa(t-t_{0})}v(t_{0})}{\gamma^{2}(1-{\rm e}^{-\kappa(t-t_{0})})}.$
(3.5)

Under a log transformation, the dynamics of the stock, $X(t)=\log S(t),$
belongs to the class of affine jump-diffusion processes, and it reads:

$\displaystyle\begin{array}[]{l}{\rm
d}X(t)=\left(r-\frac{1}{2}v(t)-\lambda\mathbb{E}[{\rm e}^{J}-1]\right){\rm
d}t+\sqrt{v(t)}{\rm d}W_{x}(t)+J{\rm d}X_{\mathcal{P}}(t),\end{array}$

so that, given Theorem 2.1, we can derive the corresponding, randomized, ChF,
which is given by the following expression:

$\displaystyle\phi_{X}({u};t,T)=\sum_{n=1}^{N}\omega_{n}{\rm
e}^{iuX(t)+\bar{C}(u,\tau;\theta_{n})v(t)+\bar{A}_{B}({u},\tau;\theta_{n})}+\epsilon_{N},$
(3.7)

for $\tau=T-t$, with complex valued functions $\bar{A}_{B}(u;\tau,\theta_{n})$
and $\bar{C}(u;\tau,\theta_{m})$, given in A.3. The third argument
$\theta_{n}$ represents a particular realization of model parameter,
$\vartheta\in\\{v_{0},\lambda,\gamma,\rho,\mu_{J},\sigma_{J},\kappa,\bar{v}\\}$,
that we wish to randomize.

Once the RAnD ChF is derived, it is worth checking the convergence concerning
the number of summation terms, $N$, in (3.7). As an example, in Figure 2 the
convergence for randomization of the vol-vol parameter, $\gamma$, with two
distributions, uniform, $\mathcal{U}$, and gamma distribution, $\Gamma$, is
depicted. The convergence rate depends on the distribution itself but also on
its parameters. In both cases, the convergence rate is excellent; however, as
expected, it is much faster for a uniform distribution than for the gamma
distribution. In both cases, satisfactory results will be accomplished for
$N=5$. Randomization of other model parameters has shown equivalent
convergence patterns.

Figure 2: Convergence results for randomization of vol-of-vol parameter,
$\gamma$ under the Bates model with: $T=31/365$, $r=0$, $\mu_{J}=-0.25$,
$\sigma_{J}=0.05$, $\lambda=0.1$, $\kappa=0.5$, $\gamma=0.72$, $\bar{v}=0.1$,
$\rho=-0.85$ and $v_{0}=0.0625$. Left:
$\gamma\sim\mathcal{U}(\hat{a},\hat{b})$ for varying $\hat{a}$ and $\hat{b}$;
Right: $\gamma\sim\Gamma(\hat{a},\hat{b})$ for varying $\hat{a}$ and
$\hat{b}$.

The impact of the randomization for the Bates model can be analyzed with the
ChF in (3.7). In particular, the time evolution of the implied volatility
surface for randomized vol-vol, $\gamma$, and the jump amplitude parameter,
$\mu_{J}$, is presented in Figure 3. The results are insightful, showing that
the RAnD method can magnify the skewness of short-maturity options.
Intriguingly, the evolution of the skew differs for both parameters, i.e.,
randomized $\gamma$ allows for slow skew decay towards the benchmark- the
basic Bates model, while the randomization of $\mu_{J}$ generates much more
rapid decay. Depending on the market data, both shapes may be desired.

Figure 3: Implied volatility surface for RAnD Bates model. Left panel:
randomized vol-vol, $\gamma~{}\sim\Gamma(1,0.5)$. Right panel: randomized
jump’s mean, $\mu_{J}~{}\sim\mathcal{N}(-0.1,0.2)$. Other model parameters are
$r=0$, $\mu_{J}=-0.1$, $\sigma_{J}=0.06$, $\lambda=0.08$, $\kappa=0.5$,
$\gamma=0.5$, $\bar{v}=0.13$, $\rho=-0.7$, $T=1/12$, and $v_{0}=0.13$.

Because of the flexibility of the RAnD method, the randomization of all the
Bates model parameters can be analyzed. As discussed earlier, one can specify
a variety of distributions (see Table 1) for the model parameters. In Figure 4
and Figure 5, the impact of randomization with different distributions and
parameters is compared against the benchmark. Randomization enables flexible
control of the implied volatility shapes. In this experiment we choose the
reference model parameters ($r=0.0$, $\mu_{J}=-0.1$, $\sigma_{J}=0.06$,
$\lambda=0.08$, $\kappa=0.5$, $\gamma=0.5$, $\bar{v}=0.13$, $\rho=-0.7$,
$v_{0}=0.13$ and $T=1/12$), and perform one-by-one randomization.

Starting with $\gamma$, in the left panel of Figure 4, we report a significant
impact on the curvature and the skew. Randomization of $v_{0}$ (right panel of
Figure 4) mainly affects the implied volatility level; however, some curvature
effect is also present. The randomization of the correlation (left panel of
Figure 5) generates the rotation of the implied volatilities around the ATM
level. Another striking effect is reported for $\mu_{J}$ where only the
volatilities far from the ATM level are affected, i.e., one can control either
the left or the right tail while keeping the ATM level fixed. The remaining
parameters are presented in Figure 11 and Figure 12 in A.6.

Figure 4: Impact of randomized parameters on implied volatilities. Left:
randomized vol-vol, $\gamma$. Right: randomized initial vol, $v_{0}$.

Figure 5: Impact of randomized parameters on implied volatilities. Left:
randomized correlation, $\rho$. Right: randomized jump’s mean, $\mu_{J}$.

###  3.1 Vix option pricing under the randomized Heston model with jumps

For a given fixed time-horizon $[t,T],$ the volatility index of an asset
$S(t)$, denoted as $\mathrm{vix}(t,T)$ is defined as:

$\displaystyle{\overline{\mathrm{vix}}}^{2}(t,T)=100^{2}\times\frac{-2}{T-t}\mathbb{E}_{t}\left[\log\frac{S(T)}{S(t)}\right],$
(3.8)

where $\mathbb{E}_{t}[\cdot]$ indicates the expectation taken under under the
risk-neutral measure $\mathbb{Q}$ and the natural filtration $\mathcal{F}(t)$.
The VIX formulation in (3.8) represents the annualized square root of the
price of a contract, where the most liquid options are for $T-t=\text{30
days}.$ The scaling with $100$ relates to a percentage representation of the
volatility.

Depending on the dynamics chosen for stock $S(t)$, the expression in (3.8) may
be derived in closed-form. In particular, for the case of the Bates model: the
closed-form expression for the expectation and the probability density
function, $f_{\mathrm{vix}}(x;t,T)$, are known explicitly and are presented in
Lemma 3.1 below.

######  Lemma 3.1.

Under the Bates model (3.3), the VIX defined in (3.8) is expressed by:

$\displaystyle\overline{\mathrm{vix}}^{2}(t,T)$ $\displaystyle=$
$\displaystyle 100^{2}\times\mathrm{vix}^{2}(t,T),$ (3.9)
$\displaystyle\mathrm{vix}^{2}(t,T)$ $\displaystyle=$ $\displaystyle
a(t,T)v(t)+b(t,T)+c,$ (3.10)

where

$a(t,T)=\frac{1-{\rm
e}^{-\kappa(T-t)}}{\kappa(T-t)},\;\;\;b(t,T)=\bar{v}\left(1-a(t,T)\right),\;\;c=2\lambda\left({\rm
e}^{\mu_{J}+\frac{1}{2}\sigma_{J}^{2}}-\mu_{J}-1\right),$

and $v(t)$ follows a non-central chi-square distribution, defined by the CIR
process in (3.3). $\lambda$ represents the Poisson intensity, and $\mu_{J}$
and variance $\sigma^{2}_{J}$ are the parameters for the jumps magnitude,
$J\sim\mathcal{N}(\mu_{J},\sigma_{J}^{2})$. The PDF of $\mathrm{vix}$ in
(3.10), is given by:

$\displaystyle
f_{\mathrm{vix}}(x;t,T)=2\alpha_{1}xf_{\chi^{2}(\delta,\bar{\kappa}(t_{0},t))}\left(\alpha_{1}(x^{2}-\alpha_{2})\right),$
(3.11)

where $\alpha_{1}=\frac{1}{a(t,T)\bar{c}(t_{0},t)}$ and $\alpha_{2}=b(t,T)+c$
with $\delta$, $\bar{\kappa}(\cdot,\cdot)$ and $\bar{c}(\cdot,\cdot)$ defined
in (3.5). Finally, the ChF for $\mathrm{vix}^{2}(t,T)$ is given by:

$\displaystyle\phi_{\mathrm{vix}^{2}(t,T)}(u)={\rm
e}^{iu(b(t,T)+c)}\left({\frac{\alpha_{1}}{\alpha_{1}-2iu}}\right)^{\frac{1}{2}\delta}\exp\left(\frac{iua(t,T)\bar{c}(t_{0},t)\bar{\kappa}(t_{0},t)}{1-2iua(t,T)\bar{c}(t_{0},t)}\right).$
(3.12)

###### Proof.

The proof can be found in A.4. ∎

Given the results from Lemma 3.1 and the availability of the ChF (3.12) for
VIX defined in (3.8), we consider a European option with strike price $K$ and
expiry $T$, whose terminal payoff is described as
$\max(\mathrm{vix}(T,T+\delta T)-K,0)$, with $\delta T$ is equal to $30$ days.
According to the risk-neutral valuation theory, the option price at time $t,$
denoted by $V_{\mathrm{vix}}(t),$ can be expressed as the discounted
conditional expectation of the terminal payoff under the risk-neutral measure,
$\mathbb{Q}$:

$\displaystyle V_{\mathrm{vix}}(t)$ $\displaystyle=$ $\displaystyle{\rm
e}^{-r(T-t)}\mathbb{E}_{t}\big{[}\max(\overline{\mathrm{vix}}(T,T+\delta
T)-\overline{K},0)\big{]}$ (3.13) $\displaystyle=$ $\displaystyle
100\times{\rm
e}^{-r(T-t)}\int_{\mathbb{R}^{+}}\max(\sqrt{v}-{K},0)f_{\mathrm{vix}^{2}}(v;T,T+\delta
T){\rm d}v,$

for $\overline{K}=100\times K$ and $\overline{\mathrm{vix}}(\cdot,\cdot)$
defined in (3.9). Since the ChF for $\mathrm{vix}^{2}(\cdot;\cdot,\cdot)$ is
explicitly known (see Equation (3.12)), we utilize the Fourier expansion from
(2.11). The missing ingredients are the payoff-dependent coefficients $H_{k}$.
These coefficients are provided by Lemma 3.2.

######  Lemma 3.2.

Assuming, $a<{K}^{2}$, the $H_{k}$ coefficients, defined in (2.11) are given
by:

$\displaystyle H_{k}$ $\displaystyle=$
$\displaystyle\frac{2}{b-a}\int_{a}^{b}\max(\sqrt{y}-{K},0)\cos\left(k\pi\frac{y-a}{b-a}\right){\rm
d}y$ $\displaystyle=:$
$\displaystyle\frac{2}{b-a}I_{1}(k)-{K}\frac{2}{b-a}I_{2}(k),$

where two integral terms $I_{1}(k)$ and $I_{2}(k)$ are given by:

$I_{1}(k):=\left\\{\begin{array}[]{l
l}a_{0}\Big{[}\sin({k_{2}})\left(C_{F}\left(a_{1}\right)-C_{F}\left(a_{2}\right)\right)+\cos({k_{2}})\left(S_{F}\left(a_{2}\right)-S_{F}\left(a_{1}\right)\right)\Big{]}\\\
+a_{5}\Big{[}\sqrt{b}\sin(a_{3})-{K}\sin\left(a_{4}\right)\Big{]},&k\neq 0,\\\
&\\\ \frac{2}{3}(b^{\frac{3}{2}}-{K}^{3}),&k=0,\end{array}\right.$

and

$I_{2}(k):=\left\\{\begin{array}[]{l
l}a_{5}\left[\displaystyle\sin\left(a_{3}\right)+\sin\left(-a_{4}\right)\right],&k\neq
0,\\\ &\\\ (b-K^{2}),&k=0,\end{array}\right.$

with $k_{1}=\frac{k\pi}{b-a}$, $k_{2}=\frac{ak\pi}{b-a}$,
$a_{0}=\frac{\sqrt{\frac{\pi}{2}}}{{k_{1}}^{3/2}}$,
$a_{1}=\sqrt{bk_{1}\frac{2}{\pi}},$ $a_{2}={K}\sqrt{{k_{1}}\frac{2}{\pi}}$,
$a_{3}=b{k_{1}}-{k_{2}}$, $a_{4}={K}^{2}{k_{1}}-{k_{2}}$,
$a_{5}=\frac{1}{{k_{1}}}$ and where $C_{F}(\cdot)$ and $S_{F}(\cdot)$ are the
so-called Fresnel integrals defined as:

$C_{F}(x)=\int_{0}^{x}\cos(t^{2}){\rm
d}t,\;\;\;S_{F}(x)=\int_{0}^{x}\sin(t^{2}){\rm d}t.$

###### Proof.

After the change of variables, $x=\sqrt{y},$ the integration is
straightforward. ∎

By semi-analytic expressions for pricing options on the underlying $S(t)$, via
Equation (2.11) and on the VIX, in (3.13), we can calibrate the randomized AD
models. The pricing formulae depend on the integration ranges $[a,b]$ that
needs to be chosen carefully. It is advised to choose it based on the
cumulants, which have been derived for VIX derivatives in [22]. We report,
however, that calibration of VIX derivatives may result in rather extreme
parameters giving rise to numerical instabilities. Since we can utilize the
analytically known distribution for $\overline{\mathrm{vix}}$, the pricing may
also be performed by directly integrating the payoff function and employing
the PDF in (3.11):

$\displaystyle V_{\mathrm{vix}}(t)=100\times 2\alpha_{1}{\rm
e}^{-r(T-t)}\int_{\overline{K}}x(x-\overline{K})f_{\chi^{2}(\delta,\bar{\kappa}(t,T))}\left(\alpha_{1}(x^{2}-\alpha_{2})\right){\rm
d}x,$

with $\alpha_{1}$ and $\alpha_{2}$ in (3.11) and $\delta$, $\bar{\kappa}(t,T)$
defined in (3.5). With one of the model parameters stochastic, the RAnD
pricing equation reads:

$\displaystyle V_{\mathrm{vix}}(t)=100\times 2\alpha_{1}{\rm
e}^{-r(T-t)}\sum_{n=1}^{N}\omega_{n}\int_{\overline{K}}x(x-\overline{K})f_{\chi^{2}(\delta,\bar{\kappa}(t,T))}\left(\alpha_{1}(x^{2}-\alpha_{2});\theta_{n}\right){\rm
d}x,$

where $\theta_{n}$ in $f_{\chi^{2}(\cdot,\cdot)}\left(\cdot;\theta_{n}\right)$
indicates a particular realization of the model parameter and $\omega_{n}$
corresponds to its weight.

###  3.2 Market data experiment: calibration to S&P500 and VIX options

This section illustrates how the randomized Bates (RAnD Bates) model can fit
implied volatilities for S&P 500 and VIX. In the experiment, we consider three
randomly selected dates: 02/02/2022, 13/05/2022 and 14/07/2022, for which the
calibration exercise was performed. In addition, we have considered short
expiry options of one month in all the cases- a challenging task for
conventional stochastic volatility models. The calibration was performed in
two stages; in the first step, the pure Bates model was calibrated to European
S&P 500 options, and such attained parameters form an initial guess for the
RAnD Bates model. Then, the optimization was performed with the target
function, including implied volatilities on S&P 500 and VIX. In the
calibration, an additional weight for the ATM level was used. In this
calibration stage, we have used the randomized vol-vol parameter $\gamma$ and
randomized it with a uniformly distributed random variable,
$\gamma\sim\mathcal{U}([\hat{a},\hat{b}])$.

The calibration results are presented in Figures 6,7,8. For all three cases,
the results were compared against the Bates model (indicated by the red line).
We conclude that jumps in the Bates model allow for a good fit for European
options on the index (illustrated in the left panels). However, it is also
clear that the Bates model does not have sufficient flexibility to fit VIX
options (right panels). The standard Bates model gives results significantly
outside the bid-ask spreads for all three studied cases.

The calibrated parameters for the RAnD Bates model are reported in Table 4.
The results are excellent, i.e., for all three cases, we can calibrate to both
products simultaneously, and the implied volatilities for the RAnD Bates model
for VIX were well within the bid-ask spreads.

Figure 6: Calibration results of the RAnD Bates model. The implied
volatilities for S&P and ViX were obtained on 02/02/2022. Dotted lines
indicate bid-ask spreads. Left: S&P, Right: VIX.

Figure 7: Calibration results of the RAnD Bates model. The implied
volatilities for S&P and ViX were obtained on 14/07/2022. Dotted lines
indicate bid-ask spreads. Left: S&P, Right: VIX.

Figure 8: Calibration results of the RAnD Bates model. The implied volatilities for S&P and ViX were obtained on 13/05/2022. Dotted lines indicate bid-ask spreads. Left: S&P, Right: VIX.  Table 4: Parameters determined in calibration of S&P and VIX Calibrated RAnD Bates parameters |
---|---
date | $\kappa$ | $v_{0}$ | $\bar{v}$ | $\rho$ | $\mu_{J}$ | $\sigma_{J}$ | $\lambda$ | $\gamma$
02/02/2022 | 0.5 | $0.170^{2}$ | 0.23 | -0.65 | $-0.25$ | 0.05 | 0.25 | $\gamma\sim\mathcal{U}([0.01,2.3])$
13/05/2022 | 0.14 | $0.267^{2}$ | 0.28 | -0.8 | $-0.25$ | 0.02 | 0.1 | $\gamma\sim\mathcal{U}([0.002,2.1])$
14/07/2022 | 0.5 | $0.250^{2}$ | 0.10 | -0.85 | $-0.25$ | 0.05 | 0.15 | $\gamma\sim\mathcal{U}([0.05,1.4])$

It is important to remember that the calibration of the randomized models,
especially for models like the Bates model with multiple parameters for
stochastic volatility and jumps, is not trivial. For that reason, we have used
the standard Bates model parameters as the initial guess for the RAnD
framework- otherwise, one may encounter the well-known problem of the
possibility of ending up in a local minimum.

###  3.3 Hedging under the RAnD method

Using the pricing Equation (2.11), one can find the hedge parameters $\Delta$
and $\Gamma$ explicitly [26]. However, it is a rather involved task to compute
the sensitivity to the model parameters, or the parameters of the randomizer
need to be computed.

Typically, when the market data changes, a model needs to be re-calibrated,
requiring a trader to hedge each parameter. Applying the RAnD method allowed
us to make use of randomized processes; it led us to a marginal density for
the asset price as a convex combination of affine densities (see Remark 2).
Conversely, the option price can be represented as the same convex combination
among the corresponding option prices. Moreover, due to the linearity of the
derivative operator, the same convex combination applies to all option Greeks.
In particular, this ensures that starting from analytically tractable affine
densities, one deals with a model that preserves the analytical tractability.
Let us consider a scenario in which a model parameter $\theta$ follows a
specific distribution with a parameter $\hat{a}$. Under such a randomized
model, the sensitivity of the derivative price, $V(t_{0})$, should be
calculated to the parameter, $\hat{a}$, that drives the parameter’s
randomness. The RAnD framework uses the randomized ChF, which must be
differentiated to compute the sensitivities. By differentiation of the
pricing, Equation (2.11), for $\bar{A}:=\bar{A}({u};\tau,\theta_{n})$ and
$\bar{B}:={\bf{\bar{B}^{\rm T}}}({u};\tau,\theta_{n})$, we find:

$\frac{\partial V(t_{0})}{\partial\hat{a}}={\rm
e}^{-r(T-t_{0})}\sideset{}{{}^{\prime}}{\sum}_{k=0}^{N_{c}-1}\Re\left[\exp\left(-ik\pi\frac{a}{b-a}\right)\frac{\partial}{\partial\hat{a}}\phi_{\bf
X}\left(\frac{k\pi}{b-a};t_{0},T\right)\right]\cdot H_{k}+\epsilon_{c_{2}},$
(3.14)

with:

$\displaystyle\frac{\partial}{\partial\hat{a}}\phi_{\bf
X}\left(u;t_{0},T\right)=\sum_{n=1}^{N}\phi_{{\bf
X}|\vartheta={\theta_{n}}}(u)\left[\frac{\partial\omega_{n}}{\partial\hat{a}}+\omega_{n}\left(\frac{\partial\bar{A}}{\partial\theta_{n}}+{\bf
X}(t)\frac{\partial\bar{B}}{\partial\theta_{n}}\right)\frac{\partial\theta_{n}}{\partial\hat{a}}\right],$
(3.15)

where $\phi_{{\bf X}|\vartheta={\theta_{n}}}(u)={\rm e}^{\bar{A}({\bf
u};\tau,\theta_{n})+{\bf{\bar{B}^{\rm T}}}({\bf
u};\tau,\theta_{n}){\bf{X}}(t)}$, with $H_{k}$ the payoff coefficient
(independent of the model parameters). Error $\epsilon_{c_{1}}$ is a function
of $\epsilon_{c_{1}}$ in (2.11) and the quadrature error, $\epsilon_{N}$, in
(2.9). The expression in (3.15) involves derivatives of the ChF with respect
to the model parameters $\partial\bar{A}/\partial\theta_{n}$ and
$\partial\bar{B}/\partial\theta_{n}$, which may to be derived analytically,
and the sensitivity to the quadrature pairs, $\\{\omega_{n},\theta_{n}\\}$,
$\partial\theta_{n}/\partial\hat{a}$ and $\partial\omega_{n}/\partial\hat{a}$.
Since the computation of these pairs is not analytic, it is recommended to
compute these derivatives numerically, with, for example, finite differences:

$\frac{\zeta(\vartheta(\hat{a}+\delta_{\hat{a}}))-\zeta(\vartheta(\hat{a}-\delta_{\hat{a}}))}{2\delta_{\hat{a}}}\approx\left\\{\frac{\partial\omega_{n}}{\partial\hat{a}},\frac{\partial\theta_{n}}{\partial\hat{a}}\right\\},$

where $\vartheta(\hat{a})$ indicates the dependence of the random variable and
parameter $\hat{a}$, $\delta_{\hat{a}}$ is the shock size and
$\zeta(\vartheta):\mathbb{R}\rightarrow\\{\omega_{n},\theta_{n}\\}_{n=1}^{N}$
is defined in A.1. Due to the applied finite difference shocks to $\hat{a}$,
an additional bias will be introduced. We expect, however, this error to be of
acceptable magnitude.

###  3.4 Convergence improvements for Monte Carlo simulation with randomized
parameters

Under the RAnD framework, the randomness of the model parameters is not
modelled by an AD SDE, but it is enforced “externally” and consequently
affects the simulation with Monte Carlo methods and the convergence. In this
section, we discuss how to perform Monte Carlo sampling and propose an
efficient variance reduction technique that will improve the convergence of
the method.

Monte Carlo simulation of the randomized model implies that every simulated
path of the underlying process ${\bf X}(t)$ will be randomized in a sense that
will depend on a realization of $\vartheta$, $\theta_{n}$. Suppose we consider
a discretization of some process ${\bf X}(t)$ from (2.1) with $M$ Monte Carlo
paths given by:

$\displaystyle\small n:\;{{x}}_{n}(t+\Delta t;\theta_{n})$ $\displaystyle=$
$\displaystyle{{x}}_{n}(t;\theta_{n})+\int_{t}^{t+\Delta
t}\bar{\mu}(s,{{x}}_{n}(s);\theta_{n}){\rm d}s+\int_{t}^{t+\Delta
t}\bar{\sigma}(s,{\bf{x}}_{n}(s);\theta_{n}){\rm d}{\widetilde{{w}}_{n}}(s)$
(3.16) $\displaystyle+$ $\displaystyle\int_{t}^{t+\Delta
t}{J}_{n}(s;\theta_{n})^{\rm T}{\rm d}{x}_{\mathcal{P},n}(s),$

with ${\widetilde{{w}}_{n}}(s)$ and ${x}_{\mathcal{P},n}(s)$ representing the
realization of the Brownian motion and Poisson process, respectively. Then,
depending on the discretization method, the integrals above can be found
explicitly or approximately by either Euler, Milstein or any other higher-
order discretization scheme. Since the randomization of model parameters is
time-invariant, this implies that every discretized Monte Carlo path of (3.16)
will be driven by a different realization of the randomized parameter.
Formally, considering a simulation with $M$ paths, every $n$’th,
$n=1,\dots,M$, path will depend on the $n$’th sample from the parameter
distribution,

The computational time can be reduced by following a similar strategy as in
Section 2.3. We have two possibilities for derivative pricing under the RAnD
framework: directly using all the Monte Carlo paths and evaluate the pricing
equation, $V_{1}(t_{0})$ in (3.17), or we can employ the idea behind the RAnD
method and divide the pricing problem into $N$ sub-problems, where for each
quadrature node $\\{\theta_{n}\\}_{n=1}^{N}$, a separate simulation is
performed and weighted according to $\\{\omega_{n}\\}_{n=1}^{N}$ (see
$V_{2}(t_{0})$ in Equation (3.17)),

$\displaystyle V_{1}(t_{0})\approx\frac{{\rm
e}^{-rT}}{M}\sum_{i=1}^{M}\Pi({x}_{i}(T),T),\;\;\;V_{2}(t_{0})\approx\frac{{\rm
e}^{-rT}}{M}\sum_{n=1}^{N}\omega_{n}\sum_{i=1}^{M}\Pi({x}_{i}(T;\vartheta=\theta_{n}),T).$
(3.17)

The division into sub-problems gives stability to the model; however, the cost
is the quadrature error, affecting the overall pricing error.

Figure 9: Convergence results for an increasing number of Monte Carlo paths.
Left panel: impact of randomization on convergence rate. Right panel: RAnD
method convergence with sub-divided pricing.

In Figure 9, the numerical results for the example case discussed in Section
2.3 are presented. The left panel of Figure 9 confirms that the stability is
affected by the randomized volatility coefficient. In the right panel of
Figure 9, we see stability by dividing the problem into $N$ sub-problems.
Although the division process improves the convergence, it will negatively
impact the computational cost (the pricing needs to be performed $N$ times).
Finally, a particular division strategy will also depend on the pricing
problem.

###  3.5 RAnD method for piece-wise constant parameters

The methodology to obtain the randomized ChF can be extended to multiple or
piece-wise parameters. In such a case, the approach could be considered a
rough approximation for a parameter driven by a stochastic process. Using the
insight that under the affine models functions $\bar{A}(u,\tau)$ and $\bar{\bf
B}(u;\tau)$ in (2.9) are the solution of the Riccati-type ODEs with nonzero
initial conditions $\bar{A}(u,0)$ and $\bar{\bf B}(u,0)$, for a time-grid
$\tau_{n_{1}}\leq\tau_{n_{2}}\dots\leq\tau_{n_{m}}=\tau$, the model parameters
can be evaluated at each interval determined by $T-\tau_{i}$. Piece-wise
constant parameters imply that the characteristic function can be evaluated
recursively, i.e. at the first interval from $[0,\tau_{n_{1}})$, we use the
initial conditions $\bar{{\bf B}}({u},0)=0$ and $\bar{A}({u},0)=0$. When the
corresponding analytic solution is determined (for details see [26], p.521),
we obtain two solutions, $a_{1}$ and $c_{1}$. For interval
$[\tau_{1},\tau_{2})$, we then assign the initial conditions $\bar{{\bf
B}}({u},\tau_{1})=c_{1}$ and $\bar{A}({u},\tau_{1})=a_{1}$. This procedure
will be repeated until the last time step, where the initial values
$c_{n_{m}-1}$ and $a_{n_{m}-1}$ are used to evaluate $\bar{{\bf
B}}({u},\tau_{n_{m}})$ and $\bar{A}({u},\tau_{n_{m}})$.

Piece-wise parameters divided over $n_{m}$ intervals will require separate
handling of each interval. Corollary 3.1 provides the details.

######  Corollary 3.1 (Piece-wise random parameters).

Consider $n_{m}$ independent random variables
$\vartheta_{1},\dots,\vartheta_{n_{m}}$, and a time grid with $n_{m}$
intervals of the time interval $[t,T]$ with $\tau=T-t$. Each of the variables
corresponds to an interval. Then, the corresponding randomized ChF is given
by:

$\displaystyle\phi_{{\bf
X}}({u};t,T)=\sum_{n_{1}=1}^{N}\dots\sum_{n_{m}=1}^{N}\omega_{n_{1}}\dots\omega_{n_{m}}{\rm
e}^{\bar{A}({u};\tau,\theta_{n_{1}},\dots,\theta_{n_{m}})+{\bf{\bar{B}^{\rm
T}}}({u};\tau,\theta_{n_{1}},\dots,\theta_{n_{m}}){\bf{X}}(t)}+\bar{\epsilon}_{N},$
(3.18)

under the same assumption as specified in Theorem 2.1.

The representation above gives us an extension for $n_{m}$ independent model
parameters. The summation is over all parameters and their corresponding
quadrature points; thus, it is subject to the so-called curse of
dimensionality. Therefore, this expression is not desirable for many
stochastic parameters. However, the representation in (3.18) can be reduced
using the sparse grid approach (see [31, 14]). The significant advantage of
such a technique is that the number of grid points does not grow exponentially
with the dimension but only polynomially.

##  4 Conclusion

In this article, we have introduced the RAnD method for efficient computation
of the affine models with random parameters. The proposed framework is generic
and can be applied to any stochastic model, even outside the class of affine
diffusions. As long as the randomizing random variable gives rise to finite,
preferably closed-form, moments, one can price European-style options
efficiently. The heart of the method is formed by a few critical collocation
points to recover the characteristic function. Fast computation of the
characteristic function is possible because the method converges exponentially
in the number of expansion terms. We have shown that the randomization of
stochastic models provides a breeze of fresh air to the class of affine
models. Finally, we have applied the RAnD method to the Bates model and shown
that randomization allows for simultaneous calibration to S&P and VIX options-
a heavily desired feature of modern models.

## References

  * Baldeaux and Badran [2014] J. Baldeaux and A. Badran.  Consistent modelling of VIX and equity derivatives using a 3/2 plus jumps model.  _Applied Mathematical Finance_ , 21(4):299–312, 2014.
  * Bates [1996] D. Bates.  Jumps and stochastic volatility: Exchange rate processes implicit in Deutsche mark options.  _Review of Financial Studies_ , 9(1):69–107, 1996\.
  * Bayer et al. [2016] C. Bayer, P. Friz, and J. Gatheral.  Pricing under rough volatility.  _Quantitative Finance_ , 16(6):887–904, 2016\.
  * Brigo and Mercurio [2002] D. Brigo and F. Mercurio.  Lognormal-mixture dynamics and calibration to market volatility smiles.  _International Journal of Theoretical and Applied Finance_ , 5(04):427–446, 2002.
  * Bulirsch and Stoer [2002] R. Bulirsch and J. Stoer.  _Introduction to numerical analysis_ , volume 3.  Springer, 2002.
  * Callegaro et al. [2019] G. Callegaro, L. Fiorin, and M. Grasselli.  Quantization meets Fourier: a new technology for pricing options.  _Annals of Operations Research_ , 282(1):59–86, 2019.
  * Carr and Wu [2007] P. Carr and L. Wu.  Stochastic skew in currency options.  _Journal of Financial Economics_ , 86(1):213–247, 2007.
  * Duffie et al. [2000] D. Duffie, J. Pan, and K. Singleton.  Transform analysis and asset pricing for affine jump-diffusions.  _Econometrica_ , 68:1343–1376, 2000.
  * Favard [1935] J. Favard.  Sur les polynomes de Tchebicheff.  _C. R. Acad. Sci., Paris_ , 200:2052–2053, 1935.  ISSN 0001-4036.
  * Fouque and Saporito [2018] J.-P. Fouque and Y. F. Saporito.  Heston stochastic vol-of-vol model for joint calibration of VIX and S&P 500 options.  _Quantitative Finance_ , 18(6):1003–1016, 2018\.
  * Gatheral et al. [2020] J. Gatheral, P. Jusselin, and M. Rosenbaum.  The quadratic rough Heston model and the joint S&P 500/VIX smile calibration problem.  _arXiv preprint arXiv:2001.01789_ , 2020.
  * Golub and Welsch [1969] G. H. Golub and J. H. Welsch.  Calculation of Gauss quadrature rules.  _Mathematics of Computation_ , 23(106):221–230, 1969.
  * Goutte et al. [2017] S. Goutte, A. Ismail, and H. Pham.  Regime-switching stochastic volatility model: estimation and calibration to vix options.  _Applied Mathematical Finance_ , 24(1):38–75, 2017.
  * Grzelak [2022] L. A. Grzelak.  Sparse grid method for highly efficient computation of exposures for xVA.  _Applied Mathematics and Computation_ , 434:127446, 2022\.
  * Grzelak et al. [2019] L. A. Grzelak, J. Witteveen, M. Suarez-Taboada, and C. W. Oosterlee.  The stochastic collocation Monte Carlo sampler: highly efficient sampling from expensive distributions.  _Quantitative Finance_ , 19(2):339–356, 2019\.
  * Guyon [2020] J. Guyon.  The joint S&P 500/VIX smile calibration puzzle solved.  _Risk, April_ , 2020.
  * Heston [1993] S. Heston.  A closed-form solution for options with stochastic volatility with applications to bond and currency options.  _Rev. Fin. Studies_ , 6:327–343, 1993.
  * Horvath et al. [2020] B. Horvath, A. Jacquier, and P. Tankov.  Volatility options in rough volatility models.  _SIAM Journal on Financial Mathematics_ , 11(2):437–469, 2020.
  * Jacquier and Roome [2015] A. Jacquier and P. Roome.  Black-Scholes in a CEV random environment.  _arXiv preprint arXiv:1503.08082_ , 2015.
  * Jacquier and Shi [2019] A. Jacquier and F. Shi.  The randomized Heston model.  _SIAM Journal on Financial Mathematics_ , 10(1):89–129, 2019.
  * Jacquier et al. [2018] A. Jacquier, C. Martini, and A. Muguruza.  On VIX futures in the rough Bergomi model.  _Quantitative Finance_ , 18(1):45–61, 2018.
  * Jing et al. [2021] B. Jing, S. Li, and Y. Ma.  Consistent pricing of VIX options with the Hawkes jump-diffusion model.  _The North American Journal of Economics and Finance_ , 56:101326, 2021.
  * Kokholm and Stisen [2015] T. Kokholm and M. Stisen.  Joint pricing of VIX and SPX options with stochastic volatility and jump models.  _The Journal of Risk Finance_ , 2015.
  * Lemaire et al. [2022] V. Lemaire, T. Montes, and G. Pagès.  Stationary Heston model: calibration and pricing of exotics using product recursive quantization.  _Quantitative Finance_ , 22(4):611–629, 2022\.
  * Mechkov [2016] S. Mechkov.  “Hot-Start” initialization of the Heston model.  _Risk.net, 2016, Available at SSRN 2674858_ , 2016.
  * Oosterlee and Grzelak [2019] C. W. Oosterlee and L. A. Grzelak.  _Mathematical modeling and computation in finance_.  World Scientific, 2019.
  * Pagès and Sagna [2015] G. Pagès and A. Sagna.  Recursive marginal quantization of the Euler scheme of a diffusion process.  _Applied Mathematical Finance_ , 22(5):463–498, 2015.
  * Papanicolaou and Sircar [2014] A. Papanicolaou and R. Sircar.  A regime-switching Heston model for VIX and S&P 500 implied volatilities.  _Quantitative Finance_ , 14(10):1811–1827, 2014\.
  * Rudd et al. [2017] R. Rudd, T. A. McWalter, J. Kienitz, and E. Platen.  Fast quantization of stochastic volatility models.  _arXiv preprint arXiv:1704.06388_ , 2017.
  * Sepp [2008] A. Sepp.  Pricing options on realized variance in the Heston model with jumps in returns and volatility.  _Journal of Computational Finance_ , 11(4):33–70, 2008.
  * Smolyak [1963] S. Smolyak.  Quadrature and interpolation formulas for tensor products of certain classes of functions.  _Dokl. Akad. Nauk SSSR_ , 4:240–243, 1963.

##  Appendix A Additional materials and proofs

###  A.1 Moments and optimal collocation pairs:
$\zeta(\vartheta):\mathbb{R}\rightarrow\\{\omega_{n},\theta_{n}\\}_{n=1}^{N}$

This section discusses the mechanism for the construction of orthogonal
polynomials. Later, this methodology is applied to establish the
$\zeta(\vartheta):\mathbb{R}\rightarrow\\{\omega_{n},\theta_{n}\\}_{n=1}^{N}$
that are used in the RAnD method.

A sequence of orthogonal polynomials $\\{p_{i}\\}_{i=0}^{N}$, with
$\deg(p_{i})=i$, is said to be orthogonal in $L^{2}$ with respect to PDF
$f_{\vartheta}(\vartheta)$ of $\vartheta$, if the following holds,

$\mathbb{E}\left[p_{i}(\vartheta)p_{j}(\vartheta)\right]=\int_{\mathbb{R}}p_{i}(x)p_{j}(x)f_{\vartheta}(x){\rm
d}x=\delta_{i,j}\mathbb{E}\left[p_{i}^{2}(\vartheta)\right],\;\;\;i,j=0,\dots,N,$
(A.1)

with $\mathbb{R}$ the support of $\vartheta$, $\delta_{i,j}$ the Kronecker
delta.

An important property of orthogonal polynomials is their definition in terms
of a recurrence relation, given in the theorem below.

######  Theorem A.3 (Recurrence in orthogonal polynomials).

For any given density function $f_{\vartheta}(\cdot)$, a unique sequence of
monic orthogonal polynomials ${p}_{i}(x)$ exists, with
$\text{deg}({p}_{i}(x))=i$, which can be constructed as follows,

$\displaystyle{p}_{i+1}(x)=(x-\alpha_{i}){p}_{i}(x)-\beta_{i}{p}_{i-1},\;\;\;i=0,\dots,N-1,$
(A.2)

where ${p}_{-1}(x)\equiv 0$, ${p}_{0}(x)\equiv 1$ and where $\alpha_{i}$ and
$\beta_{i}$ is the recurrence coefficients,

$\displaystyle\alpha_{i}=\frac{\mathbb{E}[\vartheta{p}_{i}^{2}(\vartheta)]}{\mathbb{E}[{p}_{i}^{2}(\vartheta)]},\;\mbox{for
}\;i=0,\dots,N-1,\;\;\;\;\beta_{i}=\frac{\mathbb{E}[{p}_{i}^{2}(\vartheta)]}{\mathbb{E}[{p}_{i-1}^{2}(\vartheta)]},\;\mbox{for
}\;i=1,\dots,N-1,$ (A.3)

with $\beta_{0}=0$.

###### Proof.

The proof can be found in [9]. ∎

Parameters $\alpha_{i}$ and $\beta_{i}$ are entirely determined in terms of
the moments of a random variable $\vartheta$. For many densities (weight
functions in the integration in (A.4)) the three-term recurrence relation in
(A.2) of the corresponding orthogonal polynomials has been determined. In a
scenario when the probability density function $f_{\vartheta}(\cdot)$ is not
known explicitly, or its evaluation is computationally expensive. In such
circumstances, it is desirable to express $\alpha_{i}$ and $\beta_{i}$ in
(A.3) in terms of the moments of $\vartheta$ [12]. Let us consider the
monomials $m_{i}(\vartheta)=\vartheta^{i}$, and define $\mu_{i,j}$ as follows,

$\mu_{i,j}=\mathbb{E}\left[{m}_{i}(\vartheta){m}_{j}(\vartheta)\right]=\int_{\mathbb{R}}x^{i+j}f_{\vartheta}(x){\rm
d}x=\mathbb{E}[\vartheta^{i+j}],\;\;\;i,j=0,\dots,N.$ (A.4)

From all moments $\mu_{i,j}$ we construct the so-called Gram matrix
$M=\\{\mu_{i,j}\\}_{i,j=0}^{N}$, which is symmetric and contains all moments
$\\{1,\mathbb{E}[\vartheta^{1}],\dots,\mathbb{E}[\vartheta^{2N}]\\}$. Since
matrix $M$ is, by definition, positive definite [12], we decompose $M=R^{\rm
T}R$, by the Cholesky decomposition of $M.$

The following step relates the Cholesky upper triangular matrix $R$ to the
orthogonal polynomials. This relationship has been established in [12] and is
given by,

$\displaystyle\alpha_{j}=\frac{r_{j,j+1}}{r_{j,j}}-\frac{r_{j-1,j}}{r_{j-1,j-1}},\;\;\;j=1,\dots,N,\;\;\;\text{and}\;\;\;\beta_{j}=\left(\frac{r_{j+1,j+1}}{r_{j,j}}\right)^{2},\;\;\;j=1,\dots,N-1,$
(A.5)

with $r_{0,0}=1$ and $r_{0,1}=0$ and where $r_{i,j}$ is the $(i,j)$-th element
of matrix $R$. This relation gives us the manifestations for $\alpha_{j}$ and
$\beta_{j}$ when the matrix of moments has been computed.

The recurrence coefficients $\alpha_{i}$ and $\beta_{i}$ are often called the
Darboux coefficients. The choice of these coefficients guarantees that the
generated polynomials ${p}_{n}(\vartheta)$ are orthogonal.

The next step is to relate the coefficients $\alpha_{i}$ and $\beta_{i}$ to
the zeros of the orthogonal polynomials ${p}_{n}(\vartheta)$, $n=0,\dots,N$,
utilizing the so-called the eigenvalue method, presented in the theorem below.

######  Theorem A.4 (Eigenvalue method).

The zeros ${\theta_{n}}$, $n=1,\dots,N$, of the orthogonal polynomial
${p}_{N}(\vartheta)$ are the eigenvalues of the symmetric tridiagonal matrix,

$\displaystyle\widehat{J}:=\left(\begin{array}[]{ccccc}\alpha_{1}&\sqrt{\beta_{1}}&0&0&0\\\
\sqrt{\beta_{1}}&\alpha_{2}&\sqrt{\beta_{2}}&0&0\\\ &\ddots&\ddots&\ddots&\\\
0&0&\sqrt{\beta_{N-2}}&\alpha_{N-1}&\sqrt{\beta_{N-1}}\\\
0&0&0&\sqrt{\beta_{N-1}}&\alpha_{N}\\\ \end{array}\right),$ (A.11)

i.e., $\boldsymbol{\theta}=(\theta_{1},\dots,\theta_{N})^{\rm T}$ is a vector
of eigenvalues satisfying $\widehat{J}{\bf v}=\theta_{i}{\bf v}$ with
$i=1,\dots,N$ for any real vector ${\bf v}$, with $\alpha_{i}$ and $\beta_{i}$
being the coefficients of the three-term recurrence relation (A.2).

###### Proof.

The proof can be found in [12]. ∎

Based on the coefficients $\alpha_{i}$ and $\beta_{i}$, the collocation points
$\theta_{n}$ are the eigenvalues of matrix (A.11). Eigenvalue calculation for
a tridiagonal matrix in Theorem A.4 is performed by, e.g. the Lanczos
algorithm.

The final step is to determine the corresponding weights, $\omega_{n}$. This
can be done based on the Theorem A.5 below,

######  Theorem A.5 (Optimal weights $\omega_{n}$).

Let ${\bf v}^{(n)}:=[v_{1}^{(n)},\dots,v_{N}^{(i)}]^{\rm T}$ be an eigenvector
of $\widehat{J}$ in (A.11) for the eigenvalue $\theta_{n}$, $\widehat{J}{\bf
v}^{(n)}=\theta_{n}{\bf v}^{(n)}$. Then, the weights $\omega_{n}$ are given
by: $\omega_{n}=({\bf v}^{(n)}_{1})^{2},\;\;\;n=1,\dots,N.$

###### Proof.

Proof can be found in [5] (p.179). ∎

The pseudo algorithm to compute the pairs
$\\{\omega_{n},\theta_{n}\\}_{n=1}^{N}$ is provided in Algorithm 1 below.

// Construction of the matrix $M$ with size $(N+1\times N+1)$ :

for  _$i=1\ldots N+1$_ do

for  _$j=1\ldots N+1$_ do

$M[i,j]=\mathbb{E}[\vartheta^{i+j-2}]$

// Calculate upper diagonal Cholesky matrix from $M$:

Result: $M=R^{\rm T}R$

// Find $\boldsymbol{\alpha}$ and $\boldsymbol{\beta}$ from matrix $R$,
according to Equation (A.5)

$\alpha[1]=R[1,2]/R[1,1]$

$\beta[1]=(R[2,2]/R[1,1])^{2}$

for  _$i=2:N-1$_ do

$\alpha[i]=R[i,i+1]/R[i,i]-R[i-1,i]/R[i-1,i-1];$

$\beta[i]=(R[i+1,i+1]/R[i,i])^{2};$

$\alpha[N]=R[N,N+1]/R[N,N]-R[N-1,N]/R[N-1,N-1];$

// Construction of the matrix $\widehat{J}$, according to Equation (A.11):

$J=zeros(N,N),J[1,1]=\alpha[1],J[1,2]=\sqrt{\beta(1)}$

for  _$i=2\ldots N-1$_ do

$J[i,i-1]=\sqrt{\beta[i-1]},J[i,i]=\alpha[i],J[i,i+1]=\sqrt{\beta[i]}$

$J[N,N]=\alpha[N],J[N,N-1]=\sqrt{\beta[N-1]};$

// Find the collocation points ${\theta}_{n}$ and the corresponding weights
${\omega}_{n}$

${\bf\theta}$ is computed from ${\bf x}=EigenValues(J)$

${\bf\omega}$ is obtained by squaring the first row from eigenvalues vector

Result: $\\{\omega_{n},\theta_{n}\\}_{n=1}^{N}$

Algorithm 1 Calculation of optimal quadrature weights and nodes
$\\{\omega_{n},\theta_{n}\\}_{n=1}^{N}$, The accompanying Python and MATLAB
codes can be found at https://github.com/LechGrzelak/Randomization.

###  A.2 Impact of randomizers on implied volatilities

In this section, we illustrate the impact of randomizers on the implied
volatilities under the Black-Scholes model. In particular we consider
$\sigma\sim~{}\Gamma(\hat{a},\hat{b})$ and
$\sigma\sim\hat{a}\chi^{2}(\hat{b},\hat{c})$. The results are reported in
Figure 10.

Figure 10: Implied volatility surface for the RAnD BS model for different
randomizers of $\sigma.$

###  A.3 ChF for the Bates model

For the Bates model, the complex valued functions $\bar{C}(u,\tau)$,
$\bar{A}_{H}(u,\tau)$ are known explicitly and are given by:

$\displaystyle\bar{C}(u,\tau)$ $\displaystyle=$ $\displaystyle\frac{1-{\rm
e}^{-D\tau}}{\gamma^{2}(1-g{\rm e}^{-D\tau})}\left(\kappa-\gamma\rho
iu-D\right),$ $\displaystyle\bar{A}_{H}(u,\tau)$ $\displaystyle=$
$\displaystyle
r(iu-1)\tau+\frac{\kappa\bar{v}\tau}{\gamma^{2}}\left(\kappa-\gamma\rho
iu-D\right)-\frac{2\kappa\bar{v}}{\gamma^{2}}\log\left(\frac{1-g{\rm
e}^{-D\tau}}{1-g}\right),$ $\displaystyle\bar{A}_{B}(u,\tau)$ $\displaystyle=$
$\displaystyle\bar{A}_{H}({u},\tau)-\lambda iu\tau\left({\rm
e}^{\mu_{J}+\frac{1}{2}\sigma_{J}^{2}}-1\right)+\lambda\tau\left({\rm
e}^{iu\mu_{J}-\frac{1}{2}\sigma_{J}^{2}u^{2}}-1\right).$

for $D=\sqrt{(\kappa-\gamma\rho iu)^{2}+(u^{2}+iu)\gamma^{2}}$ and
$g=\displaystyle\frac{\kappa-\gamma\rho iu-D}{\kappa-\gamma\rho iu+D}$ and
where $\tau=T-t.$

###  A.4 Proof of Lemma 3.1

Using the dynamics of the Bates model in Equation (3.3) and the definition of
the VIX (3.8) we have:

$\displaystyle\log\frac{S(T)}{S(t)}=-\lambda\mathbb{E}[{\rm
e}^{J}-1](T-t)-\frac{1}{2}\int_{t}^{T}v(u){\rm d}u+\int_{t}^{T}\sqrt{v(u)}{\rm
d}W(u)+\int_{t}^{T}J{\rm d}X_{\mathcal{P}}(u),$ (A.12)

with $\mathbb{E}_{t}[{\rm e}^{J}]={\rm e}^{\mu_{J}+\frac{1}{2}\sigma_{J}^{2}}$
and all the model parameters as specified in (3.3). By taking the
expectations, we find,

$\displaystyle\mathbb{E}_{t}\left[\log\frac{S(T)}{S(t)}\right]=-\lambda\left({\rm
e}^{\mu_{J}+\frac{1}{2}\sigma_{J}^{2}}-\mu_{J}-1\right)(T-t)-\frac{1}{2}\int_{t}^{T}\mathbb{E}_{t}[v(u)]{\rm
d}u.$ (A.13)

The expectation, conditioned on time $t$, of the CIR process is known
analytically and it is given by:

$\mathbb{E}_{t}[v(T)]=\bar{c}(t,T)(\delta+\bar{\kappa}(t,T)),$
$\bar{c}(t,T)=\frac{1}{4\kappa}\gamma^{2}(1-{\rm
e}^{-\kappa(t-T)}),\;\;\delta=\frac{4\kappa\bar{v}}{\gamma^{2}},\;\;\bar{\kappa}(t,T)=\frac{4\kappa
v(t){\rm e}^{-\kappa(T-t)}}{\gamma^{2}(1-{\rm e}^{-\kappa(T-t)})}.$ (A.14)

By integration over the expectation in Equation (A.13) and substitution to
definition of VIX in (3.8) we find:

$\displaystyle{{\mathrm{vix}}}^{2}(t,T)=100^{2}\times\frac{-2}{T-t}\mathbb{E}_{t}\left[\log\frac{S(T)}{S(t)}\right]=a(t,T)v(t)+b(t,T)+c,$
(A.15)

with

$a(t,T)=\frac{1-{\rm
e}^{-\kappa(T-t)}}{\kappa(T-t)},\;\;\;b(t,T)=\bar{v}\left(1-a(t,T)\right),\;\;c=2\lambda\left({\rm
e}^{\mu_{J}+\frac{1}{2}\sigma_{J}^{2}}-\mu_{J}-1\right).$

Since the CIR process $v(t)|v(t_{0}),\;t>0$, follows a scaled, with
$\bar{c}(t_{0},t)$, non-central chi-square distribution,
$\chi^{2}(\delta,\bar{\kappa}(t_{0},t))$, where $\delta$ is the “degrees of
freedom” parameter and the noncentrality parameter is $\bar{\kappa}(t_{0},t)$,
i.e.,

$\displaystyle
v(t)|v(t_{0})\sim\bar{c}(t_{0},t)\chi^{2}\left(\delta,\bar{\kappa}(t_{0},t)\right),$
$\displaystyle t>t_{0}.$ (A.16)

Therefore the CDF for $\mathrm{VIX}$ is given by:

$\displaystyle\mathbb{P}\left[\mathrm{vix}(t,T)\leq x\right]$ $\displaystyle=$
$\displaystyle\mathbb{P}\left[\sqrt{a(t,T)v(t)+b(t,T)+c}\leq x\right]$
$\displaystyle=$
$\displaystyle\mathbb{P}\left[v(t)\leq\frac{1}{a(t,T)}\left(x^{2}-b(t,T)-c\right)\right]=F_{v(t)}\left(\frac{1}{a(t,T)}\left(x^{2}-b(t,T)-c\right)\right).$

Since the unconditional CDF for $v(t)$ is known explicitly, we have

$F_{v(t)}(x)=\mathbb{Q}[v(t)\leq
x]=\mathbb{Q}\left[\chi^{2}\left(\delta,\bar{\kappa}(t_{0},t)\right)\leq\frac{x}{\bar{c}(t_{0},t)}\right]=F_{\chi^{2}(\delta,\bar{\kappa}(t_{0},t))}\left(\frac{x}{\bar{c}(t_{0},t)}\right),$
(A.17)

thus finally we have:

$\displaystyle
F_{\mathrm{vix}}(x)=F_{\chi^{2}(\delta,\bar{\kappa}(t_{0},t))}\left(\frac{x^{2}-b(t,T)-c}{a(t,T)\bar{c}(t_{0},t)}\right).$

By differentiation, the probability density function reads:

$\displaystyle f_{\mathrm{vix}}(x)$ $\displaystyle:=$ $\displaystyle\frac{{\rm
d}}{{\rm
d}x}F_{\chi^{2}(\delta,\bar{\kappa}(t_{0},t))}\left(\frac{x^{2}-b(t,T)-c}{a(t,T)\bar{c}(t_{0},t)}\right)$
$\displaystyle=$
$\displaystyle\frac{2x}{a(t,T)\bar{c}(t_{0},t)}f_{\chi^{2}(\delta,\bar{\kappa}(t_{0},t))}\left(\frac{x^{2}-b(t,T)-c}{a(t,T)\bar{c}(t_{0},t)}\right)=:2\alpha_{1}xf_{\chi^{2}(\delta,\bar{\kappa}(t_{0},t))}\left(\alpha_{1}(x^{2}-\alpha_{2})\right),$

with $\alpha_{1}=\frac{1}{a(t,T)\bar{c}(t_{0},t)}$ and $\alpha_{2}=b(t,T)+c$.

Having all the ingredients at hand we derive the ChF for
$\mathrm{vix}^{2}(t,T)$. By definition we find:

$\displaystyle\phi_{\mathrm{vix}^{2}(t,T)}(u)$ $\displaystyle=$
$\displaystyle\mathbb{E}_{t_{0}}\left[{\rm
e}^{iu(a(t,T)v(t)+b(t,T)+c)}\right]$ $\displaystyle=$ $\displaystyle{\rm
e}^{iu(b(t,T)+c)}\mathbb{E}_{t_{0}}\left[{\rm e}^{iua(t,T)v(t)}\right]={\rm
e}^{iu(b(t,T)+c)}\mathcal{M}_{v(t)}(iua(t,T)),$

where $\mathcal{M}_{v(t)}(iua(t,T))$ is the moment-generating function of
$v(t)$, known explicitly (see [26] p.314):

$\mathcal{M}_{v(t)}(u)=\mathbb{E}_{t_{0}}\left[{\rm
e}^{uv(t)}\right]=\left(\frac{1}{1-2u\bar{c}(t_{0},t)}\right)^{\frac{1}{2}\delta}\exp\left(\frac{u\bar{c}(t_{0},t)\bar{\kappa}(t_{0},t)}{1-2u\bar{c}(t_{0},t)}\right).$

###  A.5 Implied volatilities for RAnD BS model with uniform randomizer

Table 5: Maximum implied volatility, $\text{IV}(\cdot)$, error computed against Monte Carlo results with 10 million samples, with $K=S_{0}{\rm e}^{0.1\sqrt{T}\delta}$, $\delta=[-3,-2,-1,-0.5,0,0.5,1,2,3]$. $\text{error}=\max_{i}{\big{|}\text{IV}(t_{0},K_{i})-\text{IV}_{\text{ref}}(t_{0},K_{i})\big{|}}$ |
---|---
$\sigma\sim\mathcal{U}([a,b])$ | $N=2$ | $N=3$ | $N=4$ | $N=5$ | $N=6$ | $N=7$ | $N=8$ | $N=9$
$T=1d$ | $0.16\;\%$ | $0.04\;\%$ | $0.02\;\%$ | $0.01\;\%$ | $0.01\;\%$ | $0.01\;\%$ | $0.01\;\%$ | $0.01\;\%$
$T=1w$ | $0.16\;\%$ | $0.04\;\%$ | $0.01\;\%$ | $0.01\;\%$ | $0.00\;\%$ | $0.00\;\%$ | $0.00\;\%$ | $0.00\;\%$
$T=2w$ | $0.16\;\%$ | $0.03\;\%$ | $0.02\;\%$ | $0.02\;\%$ | $0.02\;\%$ | $0.02\;\%$ | $0.02\;\%$ | $0.02\;\%$
$T=1m$ | $0.15\;\%$ | $0.03\;\%$ | $0.01\;\%$ | $0.01\;\%$ | $0.00\;\%$ | $0.01\;\%$ | $0.01\;\%$ | $0.01\;\%$
$T=3m$ | $0.15\;\%$ | $0.03\;\%$ | $0.01\;\%$ | $0.01\;\%$ | $0.01\;\%$ | $0.01\;\%$ | $0.01\;\%$ | $0.01\;\%$
$T=6m$ | $0.15\;\%$ | $0.03\;\%$ | $0.01\;\%$ | $0.00\;\%$ | $0.00\;\%$ | $0.00\;\%$ | $0.00\;\%$ | $0.00\;\%$
$T=12m$ | $0.15\;\%$ | $0.03\;\%$ | $0.01\;\%$ | $0.01\;\%$ | $0.01\;\%$ | $0.01\;\%$ | $0.01\;\%$ | $0.01\;\%$

###  A.6 Impact of randomization on Bates model parameters

The impact of different randomizers on implied volatilities under the Bates
model is presented in Figure 11 and Figure 12.

Figure 11: Impact of randomized parameters on implied volatlities. Left:
randomized long-term vol, $\bar{v}$. Right: randomized jump’s size vol,
$\sigma_{J}$.

Figure 12: Impact of randomized parameters on implied volatlities. Left:
randomized speed of mean reversion, $\kappa$. Right: randomized jump’s
intensity, $\lambda$.