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

#  Basket Options with Volatility Skew: Calibrating a Local Volatility Model
by Sample Rearrangement

NICOLA F. ZAUGG<EMAIL_ADDRESS>LECH A. GRZELAK<EMAIL_ADDRESS>Mathematical Institute, Utrecht University, Utrecht, the Netherlands
Financial Engineering, Rabobank, Utrecht, the Netherlands  Capital Markets
Technology, swissQuant Group AG, Zürich, Switzerland

###### Abstract

The pricing of derivatives tied to baskets of assets demands a sophisticated
framework that aligns with the available market information to capture the
intricate non-linear dependency structure among the assets. We describe the
dynamics of the multivariate process of constituents with a copula model and
propose an efficient method to extract the dependency structure from the
market. The proposed method generates coherent sets of samples of the
constituents process through systematic sampling rearrangement. These samples
are then utilized to calibrate a local volatility model (LVM) of the basket
process, which is used to price basket derivatives. We show that the method is
capable of efficiently pricing basket options based on a large number of
basket constituents, accomplishing the calibration process within a matter of
seconds, and achieving near-perfect calibration to the index options of the
market.

###### keywords:

Basket Options, Local Volatility, Correlation Structure, Copula, Rearrangement
Algorithms

††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

Basket derivatives are a common class of exotic financial instruments. These
derivatives depend on the performance of a linear combination of multiple
underlying assets, referred to as the constituents of the basket. Constituents
typically include equities, currencies, or commodities, although
theoretically, they can be any financial asset. Basket products take various
forms, ranging from vanilla options to intricate structured products that
involve more advanced payoffs and additional interest rate components. They
enable trading participants to gain exposure to multiple assets with a single
trade, thereby carrying a significant amount of correlation risk.

Basket derivatives are exotic derivatives and are commonly traded over-the-
counter (OTC), meaning they are usually illiquid instruments. Some exceptions
are derivatives on stock indices, commodity indices, or interest rate spreads
(hereafter denoted as “index derivatives”), which are also classified as
basket products. Vanilla options on such baskets often exhibit high liquidity
and are actively traded on exchanges. For instance, the equity index options
volume exceeded the single stock option volume by more than 400% in 2022 [6].
Considering a substantial portion of basket products are traded over-the-
counter, the development of an efficient and consistent risk-neutral pricing
framework becomes crucial for traders trying to offer a fair price to the
market and hedge the risks accordingly [15]. Financial modelers encounter two
primary challenges in implementing such frameworks. The first challenge
involves addressing the inherent high dimensionality of the problem. Because
each constituent contributes to the overall performance of the derivative,
numerical approximations rapidly become computationally expensive as their
calibration and evaluation time increases substantially. An effective pricing
framework frequently incorporates efficient approximations and dimensionality
reductions to enhance performance as needed.

Secondly, a pricing framework has to be calibrated to the available market
data in a consistent way. Once calibrated, the model should replicate the
market price of all liquid financial instruments in the scope of the model. It
requires enough flexibility to price constituents’ implied volatility skew as
well as the index volatility skew, for which liquid market data is available.

Suppose that an index $I(t)$ consists of $N$ constituents
$S_{1}(t),S_{2}(t),\dots S_{N}(t)$. The price process of the index at time $t$
is then given by

$I(t)=\sum_{n=1}^{N}w_{n}S_{n}(t),$ (1.1)

where $w_{n}$ are the weights of the index. Since the prices of the
constituents can be scaled individually, let us, without loss of generality,
assume that $w_{n}=1,\;\forall n\leq N$, such that the index is the sum of all
constituents. Suppose that $N_{I}=\mathbb{N}_{\leq N}$ is the set of all
integers below or equal to $N$. A basket $B(t)$ is formed from any subset
$N_{B}\subset N_{I}$ of constituents and is also equal to the sum of its
(possibly weighted) constituents

$B(t)=\sum_{n\in N_{B}}S_{n}(t).$ (1.2)

In a risk-neutral setting, the constituents $S_{n}(t)$ are random processes on
a probability space $(\Omega,\mathcal{F}(t),\mathbb{Q})$, where the dynamics
of $S_{n}(t)$ are determined by the market prices of vanilla options on the
constituents.

To price a derivative whose payoff depends on Equation 1.2, we require
information about the distribution of the process $B(t)$. Since the marginal
distribution functions of $S_{n}(t)$ are determined by the individual
constituents’ options, the difficulty in modeling the distribution of $B(t)$
lies solely in the joint distribution of the constituents. While the essence
of such a pricing framework is _not_ to price derivatives on index $I(t)$, the
availability of liquid index options provides certain market-implied
information about the joint distribution between the individual constituents.
Extracting this information and embedding it in the model allows us to price
basket options on any subset of constituents consistently.

A practical approach to modeling basket distributions is the so-called
_correlated local volatility model_ [20]. The model assumes correlated local
volatility diffusion processes for the individual constituents, where the
correlations are based on historical data or calibrated to index ATM (At-The-
Money) volatilities. The basket process is then simulated using Monte-Carlo
methods or derived from some analytical approximations (See, for instance,
[14]).

It is widely acknowledged that such an approach suffers from a key pitfall [4,
7, 11]. As mentioned earlier, to properly capture the volatility skew of the
index, the probability distribution of $I(t)$ must align with the option-
implied risk-neutral probability distribution. However, the probability
distribution assumed by the model is usually inconsistent with its implied
distribution by the market. This inconsistency means that the index option
smile/skew cannot be faithfully replicated [7, 11]. The primary issue lies in
the fact that the correlated Brownian motions, which model the linear
dependency between the assets, are insufficient to explain the complex
dependency structure among the constituents. Consequently, this leads to
inconsistencies in the smile/skew, as the calibration procedure can only match
ATM volatilities and the model will lead to large mispricing in volatile
market conditions [2]. A detailed discussion of this issue will be provided in
Section 2.

The most notable class of models aimed at resolving the issue of the index
skew is the class of so-called local correlation models, extensively discussed
in the literature [11, 16, 8]. In these models, the correlation between
constituents at time $t$ depends on the state $S_{n}(t)$ or $I(t)$. While
these models effectively address the problem of non-constant correlation, they
face significant drawbacks due to their high computational complexity.
Specifically, each constituent of the index needs to be simulated using the
complex non-constant correlation in a numerical scheme, which leads to
inefficiencies. Efforts to resolve this have been made in [10], which promise
lower computational costs, but in turn come at the cost of imperfect
calibration to the index skew.

Parallel to local correlation models are stochastic correlation models,
grounded in correlation as a Jacobi process and further enriched with jumps
[21]. While these models provide insights into the precise nature of the
correlation skew among assets, their practical usability is limited due to
lower flexibility with similar computational complexity, compared to the local
correlation model.

While local and stochastic correlation models can reprice the index options by
enhancing the flexibility of the model with extra parameters, at the cost of
computational complexity, they only partially capture the underlying issue:
correlation is not sufficient to model the joint distribution. Aiming to
address this issue more fundamentally, copula models are proposed as an
alternative to correlation models (e.g., [17, 19, 12]). Copula models leverage
a mathematical tool to describe complex dependencies among random variables,
capturing the intricate nature of joint distributions. While these models are
mathematically intriguing, explicit copula models are limited to bivariate
copulae or lack dynamic form and provide no straightforward method for
calibration other than a brute-force search for parameters [2].

A feature shared by all mentioned approaches so far is that they propose a
multivariate model, meaning that each constituent asset price is modeled
explicitly by a stochastic differential equation (SDE). For pricing
derivatives on the process $B(t)$, however, it would suffice to model the
dynamics of $B(t)$ with a one-dimensional SDE, for instance, by deriving a
local volatility model (LVM)

$\frac{{\rm d}B(t)}{B(t)}=r{\rm d}t+\sigma_{LV}(t,B(t)){\rm d}W(t),\;\;t\geq
0,$ (1.3)

where the local volatility function $\sigma_{LV}(t,B(t))$ is given by Dupires
formula111Alternatively, an equivalent formulation using implied volatility
can be used. See [13, Section 4.3.1] [13],

$\sigma_{LV}^{2}(t,k)=\frac{\frac{\partial}{\partial
t}\left[e^{-rt}\int_{k}^{\infty}(y-k)f_{B(t)}(y){\rm
d}y\right]+rk\left(\int_{-\infty}^{k}f_{B(t)}(y){\rm
d}y-1\right)}{\frac{1}{2}k^{2}f_{B(t)}(k)},$ (1.4)

with $f_{B(t)}$ the probability density function (PDF) of the basket $B(t)$ at
time $t$ and $r$ the prevailing interest rate. While the calibration of this
model still requires a suitable multidimensional model to derive $f_{B(t)}$
from the component processes, this model of $B(t)$ has an advantage: The
estimation of $f_{B(t)}$ at a fixed time $t$ can be done statically, meaning
that $f_{B(t)}$ can be determined _independently_ for a time discretization
$t\in\\{t_{1},t_{2},\dots,t_{E}\\}$ with $E$ time steps. The probability
densities for each $t$ can then be interpolated and combined to determine the
local volatility function in continuous time. Such a static calibration method
was proposed by Grzelak et al. [7] where the authors aim to estimate
$f_{B(t)}$ by constructing a set of samples of $B(t)$ for all maturity times
$t$ available in the market. The sample sets are initialized on a constituent
level with an appropriate sampling scheme to enforce a dependency between the
constituents and then combined to create a sample of the basket $B(t)$. The
probability density function $f_{B(t)}$ is then estimated from the empirical
distribution of the samples. The authors thus show that the problem of
calibrating the local volatility function of Equation 1.4 is reduced to
constructing (static) sample sets of $B(t)$ using an appropriate sampling
scheme.

A particular sampling scheme that is suitable for this type of problem are
rearrangement algorithms. Rearrangement algorithms are used to construct
samples of multivariate random variables when the marginal distributions are
given, but only limited information about the joint probability distribution
is available. These algorithms have found success in mathematical finance
before, most notably in the context of risk management [5]. Bernard et al. [3]
showed that the algorithms are applicable to the index/basket option setting
as the constituent PDF $f_{S_{n}(t)}$ are given by the constituents option
market. The constituent samples are first initialized marginally and are then
rearranged to ensure that the resulting empirical distribution of the index
samples (obtained by summing up the rearranged samples) is “admissible”,
meaning that it is close to the probability distribution obtained from the
market. The empirical distribution of $B(t)$ is then available by summing up
the appropriate constituent samples. Rearrangement algorithms implicitly
extract a dependency structure from the available market information by
finding an optimal arrangement, such that the samples reflect the available
information provided by the cumulative distribution function (CDF) $F_{I(t)}$
of the index.

In this paper, we will describe the theoretical and practical foundation of
utilizing rearrangement algorithms in the context of basket derivatives. Using
a copula model we first describe a risk-neutral pricing framework where each
constituent process is marginally given by a local volatility model. We then
show how rearrangement algorithm are used to extract a market-implied
dependency structure to calibrate the copula of the processes. Building on the
theoretical foundation, we introduce our own algorithm which we call the
Iterative-Sort-Mix algorithm (ISM). This algorithm generates constituent
samples for an index with relatively low computation requirements. In the
final step, the samples of the algorithm are used to calibrate the local
volatility model of the basket $B(t)$ given by (1.3), which can then be used
to price vanilla and exotic payoffs on the basket $B(t)$. We assess the
effectiveness of the algorithm using real market data for index options on the
Dow Jones Industrial Average (DJIA) with 30 underlying constituents. The
algorithm successfully prices the volatility skew for all available time to
maturities, establishing its suitability as a pricing model for basket
derivatives.

The paper is structured as follows: in Section 2, we first describe the
general multivariate copula model and explain the important aspects of a
suitable dependency structure, particularly highlighting the inadequacy of
correlated LVMs. Furthermore, we discuss the problem of underdetermination
when calibrating a copula model to the market. In Section 3, we provide the
mathematical foundation for a rearrangement algorithm to generate the samples
required to estimate $f_{B(t)}$, and show how to apply the theory from Section
2 to rearrangement algorithms. In Section 4, we introduce the “Iterative Sort-
Mix” algorithm, an implementation of a rearrangement algorithm, which can be
used to price basket options. Finally, we apply the algorithm to actual data
in Section 5 and we draw our conclusions in Section 6.

##  2 Copula Models: Joint Distribution of Constituents

###  2.1 Model Description

In this paper, we study the joint dynamics of $N$ random processes (called
constituent processes)

$\left(S_{1}(t),S_{2}(t),\dots,S_{N}(t)\right),$ (2.1)

on a probability space $(\Omega,\mathcal{F}(t),\mathbb{Q})$, where
$\mathbb{Q}$ is the risk neutral measure. We will model the individual
marginal processes $S_{n}(t)$ with a local volatility model given by the local
volatility function $\sigma_{n,LV}(t,S_{n}(t))$ and the constant interest rate
$r$ as

$\frac{{\rm d}S_{n}(t)}{S_{n}(t)}=r{\rm d}t+\sigma_{n,LV}(t,S_{n}(t)){\rm
d}W_{n}(t).$ (2.2)

The Brownian motions $W_{n}(t)$ are dependent on each other and their
dependency structure will be defined later. Given a subset $N_{B}$ of the $N$
constituents, the price process of the basket is obtained by summing the
individual constituents of the basket.

$B(t)=\sum_{n\in N_{B}}S_{n}(t).$ (2.3)

The process $B(t)$ can also be modeled by a “direct” stochastic differential
equation, given by a local volatility model

$\frac{{\rm d}B(t)}{B(t)}=r{\rm d}t+\sigma_{LV}(t,B(t)){\rm d}W(t),\;\;t\geq
0.$ (2.4)

As we hinted in the introduction, the key to consistent basket option pricing
lies in modeling the dependency structure among constituent assets. While the
dynamics of Equation 2.2 determine the cumulative distribution function of
constituent assets, the sum’s CDF is dictated by the joint distribution of the
constituents and, consequently, by the dependency structure between the
individual assets. Suppose that $F_{S_{n}(t)}$ is the probability distribution
function of $S_{n}(t)$. To describe the joint distribution of the
constituents, we will use a (time-dependent) copula, which is given by the
function $C\colon[0,\infty)\times[0,1]^{N}\to[0,1]$, such that

$C(t,u_{1},u_{2},\dots,u_{N})=\mathbb{Q}(F_{S_{1}(t)}(S_{1}(t))<u_{1},F_{S_{2}(t)}(S_{2}(t))<u_{2},\dots,F_{S_{N}(t)}(S_{N}(t))<u_{N}).$
(2.5)

The copula function is assumed to be sufficiently smooth in the space
variables to admit a density copula for all times $t\in[0,\infty)$

$c(t,u_{1},u_{2},\dots,u_{N})=\frac{\partial^{N}C(t,u_{1},u_{2},\dots,u_{N})}{\partial
u_{1}\partial u_{2}\dots\partial u_{N}}.$ (2.6)

We will characterize the set of all copulae of the joint random vector
$(S_{1}(t),S_{2}(t),\dots,S_{N}(t))$ as $\Theta$ , such that $\theta\in\Theta$
determines a unique copula function $C_{\theta}$.

The challenge in pricing basket options lies in selecting a suitable copula
function $C$ that is coherent with the market data. Suppose that the market
has a certain view on the dependence of an index’s constituents. This
dependence is reflected in the market-implied risk-neutral probability
distribution of $I(t)$ at time $t\in[0,\infty)$. The availability of liquid
index option prices allows us to extract this probability distribution, which
we call $F_{I^{\text{Mkt}}(t)}(x)$. The copula $C$ must thus be chosen in
accordance with this market expectation 222This is equivalent to saying that
the model replicates the implied volatility skew of the market, meaning that
vanilla options prices from the model will be equal to the market prices.,
i.e.

$F_{I(t)}(x)=\mathbb{Q}(I(t)<x)=\mathbb{Q}\left(\sum_{n=1}^{N}S_{n}(t)<x\right)=F_{I^{\text{Mkt}}(t)}(x),\hskip
14.22636pt\forall x\in\mathbb{R},\;t\in[0,\infty).$ (2.7)

We call a copula $C$ which agrees with the market-implied distribution an
_admissible copula_.

######  Definition 2.1 (Admissible Copula).

Suppose that $S_{1}(t),S_{2}(t),\dots S_{N}(t)$ have a joint distribution with
copula $C$ such that Equation 2.7 holds. Then, we call the copula function $C$
an admissible copula. We denote the set of all admissible copulae as
$\Theta_{A}\subset\Theta$.

Note that the existence of an admissible copula depends on an assumption of
efficiency in the market. The distribution $F_{I^{\text{Mkt}}(t)}$ is obtained
from index option quotes, while we calibrate the dynamics of the individual
constituents to single constituents’ option quotes. We will assume that the
two sets of option quotes are consistent with each other, implying that there
is no arbitrage in the market.

######  Assumption 1 (Market Efficiency).

We will assume that the index options priced in the market are in accordance
with the vanilla options of its constituents, and thus, that an admissible
copula can be found.

So far the introduced model is a time-continuous model for the constituent and
index processes, meaning that the copula $C(t,\cdot)$ is also defined on a
time continuum $[0,\infty)$. The data from the market for $S_{n}(t)$ and
$I(t)$ to calibrate the model is generally only available for a discrete set
of market maturities $T_{E}=\\{t_{1},t_{2},\dots,t_{E}\\}$. This means that we
can only calibrate an admissible copula to this discrete set of time points
and an interpolation/extrapolation in time after calibrating the model. The
chosen time interpolation then determines the exact shape of the copula on the
time continuum. Such a time interpolation can be done on the implied
volatilities and is standard practice (see for instance [13, 1]). This
approach, however, has an impact on our model. The copula, which is defined on
a time-continuous spectrum, is thus fully determined by its value at the
discrete market maturities $T_{E}$. This means that the following additional
assumption is made:

######  Assumption 2 (Time interpolation).

If $C_{1}$ and $C_{2}$ are two copulae, we have

$C_{1}(t,\cdot)=C_{2}(t,\cdot),\;\forall t\in[0,\infty)\iff
C_{1}(t,\cdot)=C_{2}(t,\cdot),\;\forall t\in T_{E}.$

###  2.2 Intricate Dependency Structures

Before diving further into the copula model, let us first investigate why
choosing a simple dependency structure, such as in the correlated local
volatility model (cLVM) is not sufficient to price index and basket options.
The dependence of the constituents in a cLVM is directly modeled by
correlating the Brownian motions in Equation 2.2. This means that for any
distinct pair $n,m\in N_{I}$, we have

${\rm d}W_{n}(t){\rm d}W_{m}(t)=\rho_{n,m}{\rm d}t,$ (2.8)

where $\rho_{n,m}\in[-1,1]$ is the correlation333Note that the correlation
matrix has to be positive definite. between the Brownian motion. The price
processes are thus driven by correlated Brownian motions with correlation
$\rho_{n,m}$, and we can price index options with a Monte-Carlo pricer by
simulating $S_{n}(t)$ for each $n$.

We observe now, when calibrating the correlations to the index option prices,
that the fitted correlation coefficients of the constituents $\rho_{n,m}$ are
not equal when calibration to options of different strikes [4]. It appears as
if the market is assuming different correlation coefficients of the underlying
assets when considering different strikes. This observation is referred to as
“implied correlation”, which, similar to implied volatility, describes the
strike-dependent correlation coefficient observed in the market. In practice,
the existence of implied correlation means that if a model is calibrated to
ATM prices of the index option, non-ATM options will be mispriced by the
model.

Just as with the existence of an implied volatility skew, the existence of
implied correlation is not due to irrationality in the market, but rather due
to the model’s inability to capture all the complex effects present in the
market. It turns out that the dependence structure of Equation 2.8, which is
solely determined by the correlation coefficient, induces a _linear_
relationship between the assets. This linear relationship means that the
dependency of the assets does not depend on the absolute values. This is
insufficient to model dependencies for constituents. On one hand, analysis of
historical correlations shows that correlation is stronger for negative
returns than for positive returns [21, Figure 1]. Such a structure is also
expected by the market, as the market-implied expectation of the correlation
often exhibits a skew or smile [21, Figure 2].

We visualize the impact of the difference in dependency structures in Figure
1. The dependence between two assets at a fixed time $t$ can be visualized
with a scatter plot. While the Pearson’s correlation coefficient is $0.5$ in
both cases, the left shows a non-linear dependency structure. The skewed
distribution exhibits a strong left-tail correlation, meaning that the
correlation is stronger for lower values.

Figure 1: Example of a bivariate copula, visualized. Left: Stronger
correlation for low values. Right: Uniform correlation. In both cases, the
Pearson correlation is $0.5$.

It is not obvious at first sight why the exact structure of the correlation
matters for pricing options. The impact of the structure is best observed by
examining the generated PDF of the index $I(t)$ under uniform correlation vs
non-uniform correlation. Figure 2 displays the distribution of the sum of the
random variables in Figure 1 and the impact on implied volatility skews given
these correlation structures. The non-uniform correlation structure with a
stronger lower-tail correlation (Figure 1, left) will exhibit a broader lower
tail, resulting in a skew in the implied volatility.

Figure 2: Effect of difference in copula on basket pricing: We plot the PDFs
and the implied volatility skew of the sum of two random variables as in
Figure 1. The copula with a strong left tail exhibits a larger implied
volatility skew.

Given that the market expects a non-linear correlation structure, it is
therefore crucial to have a model that is flexible enough to account for this
expectation to reprice the index options. Copula models allow this flexibility
since copulae can describe any type of dependency between $N$ random
variables.

The obvious way to utilize a copula model in practice is to specify an
admissible copula $C$, which then defines the distribution of $I(t)$ and
$B(t)$ and any time $t$. Such an approach was described for instance by Lucic
[12] using product copulae. The main issue with such an approach is the
calibration of the copula. It turns out it is increasingly difficult to
specify a copula for a high number of constituents, and often, a brute-force
search algorithm is required to find the right parameters.

For this reason, rearrangement algorithms are used as an alternative to
describe the dependency. Instead of defining an explicit copula $C$ which is
suitable to model the correlation, rearrangement algorithms implicitly extract
the dependency structure and generate samples for the random vector (2.1)
under an appropriate admissible copula. These samples are sufficient to
calibrate Equation 2.2, since the cumulative distribution function $F_{B(t)}$
can be inferred directly from these samples. The exact workings of
rearrangement algorithms will be described in Section 3.

###  2.3 Uniqueness of Admissible Copula

We defined an admissible copula $C$ to be any copula for which Equation 2.7 is
given. The condition ensures that the model captures the market expectation of
the copula $C$ given the information available from the index option market.
Since the condition (2.7) can be fulfilled for more than one copula $C$, we
should consider the implications of multiple distinct admissible copulae. At
first sight, the uniqueness might not seem important, as all admissible
copulae lead to the same index probability distribution, and thus also to the
same index option prices. However, while the index option prices are the same
for a model with any admissible copula, the resulting basket probability
distributions given a basket subset $N_{B}$ are not necessarily consistent.
This means that two admissible copulae can imply different basket probability
distributions for the subset of constituents. The source of this issue is in
the fact that $F_{I^{\text{Mkt}}(t)}$ only provides partial information about
the dependency of the constituent random variables. We refer to this as the
problem of _underdetermination_ of copula $C$. To illustrate the problem, let
us examine a toy model of three uniformly distributed assets.

######  Example 2.1.

Let $t\in[0,\infty)$ be fixed and let $S_{1}(t),S_{2}(t),S_{3}(t)$ be
marginally uniformly distributed on $[0,1]$. Suppose in a first case that
$S_{1}(t)=S_{2}(t)$ and $S_{3}(t)$ is independent of them. In this case, the
copula $C_{1}(t)$ is given by

$C_{1}(t,u_{1},u_{2},u_{3})=\mathbb{Q}(S_{1}(t)\leq\max(u_{1},u_{2}))\mathbb{Q}(S_{3}(t)\leq
u_{3})=\max(u_{1},u_{2})u_{3}.$

Suppose now that $S_{1}(t)=S_{3}(t)$ and $S_{2}(t)$ is independent of them. In
this case

$C_{2}(t,u_{1},u_{2},u_{3})=\mathbb{Q}(S_{1}(t)\leq\max(u_{1},u_{3}))\mathbb{Q}(S_{2}(t)\leq
u_{2})=\max(u_{1},u_{3})u_{2}.$

The copula $C_{1}$ and $C_{2}$ are not equal to each other, but it is clear
that in either case, the distribution of $S_{1}(t)+S_{2}(t)+S_{3}(t)$ is
$2U_{1}+U_{2}$, where $U_{1},U_{2}$ are independent uniformly distributed
random variables. The copula $C_{1}$ is thus admissible if and only if $C_{2}$
is admissible. We run into an issue when considering a basket
$N_{B}=\\{1,2\\}$, as the distribution of $S_{1}(t)+S_{2}(t)$ is not the same
under $C_{1}$ and $C_{2}$. In the first case, we have

$B(t)=2U_{1},$

while in the second case,

$B(t)=U_{1}+U_{2}.$

Pricing a derivative on $B(t)$ will yield different results in the cases.

The example clearly shows that an admissible copula is not unique, and the
distinct copulae lead to different dynamics of $B(t)$. This poses the question
of how to determine which copula is the “superior” one, yielding the “right”
basket option price. Given the lack of additional information on the
dependency between the constituents, a reasonable approach is to select a
copula with a symmetry property.

######  Definition 2.2 (Symmetric Copula).

Let $C\colon[0,\infty)\times[0,1]^{N}$ be a copula. Suppose that $\pi\colon
N_{I}\to N_{I}$ is a permutation of the vector indices, and
$P_{\pi}\colon[0,1]^{N}\to[0,1]^{N}$ the function which rearranges a vector
according to $\pi$, i.e

$P_{\pi}(u_{1},u_{2},\dots,u_{N})=\left(u_{\pi(1)},u_{\pi(2)},\dots,u_{\pi(N)}\right).$

The copula $C$ is called symmetric iff:

$C(t,u_{1},u_{2},\dots,u_{N})=C(t,P_{\pi}(u_{1},u_{2},\dots,u_{N})),$

for all $t\in[0,\infty)$.

A symmetric copula is a copula where the dependency structure between any two
distinct constituents is the same. In other words, we do not favor the
dependency between any pair within the constituents of the index, since we do
not have any information on such dependency. This is in contrast to the two
copulae chosen Example 2.1, where $S_{1}(t)$ and $S_{2}(t)$ are strongly
dependent in the first case, while $S_{3}(t)$ is independent of them (and
vice-versa in Case 2).

######  Assumption 3.

We consider an admissible copula $C\in\Theta_{A}$ as “valid” if the copula is
symmetric. We denote the set of all symmetric copulae as $\Theta_{A}^{S}$.

Note that most basket option models implicitly only consider symmetric copulae
by construction, although such terminology is not explicitly introduced (for
instance [7, 12, 21]). The concept of symmetry is also discussed in Bernard et
al. [3] in the context of “entropy”.

##  3 Rearrangement Algorithms

In the preceding section, we defined the copula model for a multivariate
process for constituents of an index or basket. We demonstrated that the
difficulty of such a model is the choice of the copula $C$, such that
$(S_{1}(t),S_{2}(t),\dots,S_{N}(t))$ has a joint probability distribution with
an admissible and symmetric copula. In the remaining section, we aim to
utilize the theoretical model and convert it into a practical pricing
framework for any applicable payoff on $B(t)$. In particular, we aim to
calibrate the “direct” local volatility model for $B(t)$ given by Equation
2.4.

Rather than directly specifying a copula for the multivariate model, we will
extract the dependency structure of the multivariate model by using a
rearrangement sampling algorithm. This sampling algorithm will provide samples
for the constituents given the marginal distribution $F_{S_{n}}(t)$ of the
constituents, and the index CDF $F_{I}(t)$. The calibration of the dynamic of
$B(t)$ is then achieved by estimating the probability density function
$f_{B}(t)$ from the samples and approximating the probability distribution
using the empirical distribution. Figure 3 shows a complete overview of the
process.

Figure 3: Rearrangement algorithm for basket options: Shows how the
rearrangement algorithm is used to calibrate the local volatility model. In
the first step, the constituents’ option quotes are converted to risk-neutral
PDFs. Similarly, we obtain the index PDF from the index option quotes. The
rearrangement algorithm then generates samples for the constituents. These
samples are used to estimate $f_{B(t)}$, which in turn is used to calibrate
the local volatility model of $B(t)$.

The sampling strategies we will consider involve rearrangement algorithms. In
this section, we will formally introduce a generic rearrangement algorithm and
demonstrate how it can be utilized to generate the necessary samples.
Subsequently, in Section 4, we will introduce a specific algorithm designed
for calibrating basket models.

######  Remark 1.

The rearrangement algorithm provides samples for $B(t)$ for every $t\in
T_{E}$. We can use these samples to price payoffs on $B(t)$ which only depend
on the PDF of $B(t)$ at a single $t\in T_{E}$, such as vanilla European
options. Any path-dependent payoff requires a stochastic model of $B(t)$. The
local volatility model of Equation 2.4 is a natural choice given that the
local variance function $\sigma^{2}_{LV}$ can be estimated from the static
samples.

###  3.1 General Sample Notation

We begin by introducing notation for collections of samples for constituents,
baskets, and indices. Suppose that, for an asset $S_{n}(t)$ at some time $t\in
T_{E}$, we draw $M$ samples. We denote the obtained samples as

$s_{n}(t)=\left[s^{1}_{n}(t),s^{2}_{n}(t),\dots
s_{n}^{M}(t)\right]^{T}\in\mathbb{R}_{+}^{M}.$ (3.1)

We call the vector $s_{n}(t)$ a _constituent sample vector_ if the asset
$S_{n}(t)$ belongs to an index. Since an index consists of multiple
constituents, we require multiple constituent sample vectors to construct
samples for an index. These constituent sample vectors are then collected an
_index sample matrix_ 444We will often refer to it simply as “sample matrix”.
, denoted as the collection $[s]_{M}$,

$[s]_{M}=\left[s_{n}(t)\colon
n=1,2,\dots,N\right]=\begin{bmatrix}s_{1}^{1}(t)&s_{2}^{1}(t)&\cdots&s_{N}^{1}(t)\\\
s_{1}^{2}(t)&s_{2}^{2}(t)&\cdots&s_{N}^{2}(t)\\\
\vdots&\vdots&\ddots&\vdots\\\
s_{1}^{M}(t)&s_{2}^{M}(t)&\cdots&s_{N}^{M}(t)\end{bmatrix}\in\mathbb{R}_{+}^{M\times
N}.$ (3.2)

The index sample matrix thus contains all sample vectors from all the
constituents. Since the index $I(t)=\sum_{n\leq N}S_{n}(t)$ is a sum of all
the assets in the index, we can define the samples of the index as the sum of
the constituent samples. We define $i(t)$ as the sum of the constituent sample
vectors:

$i(t)=\sum_{n\leq N}s_{n}(t)\in\mathbb{R}_{+}^{M}.$ (3.3)

In the same fashion, we find the basket samples for an arbitrary basket
$N_{B}$ by summing up the corresponding constituent samples:

$b(t)=\sum_{n\in N_{B}}s_{n}(t)\in\mathbb{R}_{+}^{M}.$ (3.4)

Lastly, we introduce the empirical distributions given a sample matrix
$[s]_{M}$. Empirical distribution functions count the number of samples
observed in a certain region given a sample matrix. For individual
constituents, we define the marginal empirical distribution as

$\hat{F}_{[s]_{M};n}(x)=\frac{1}{M}|\\{m\leq M\colon s_{n}^{m}(t)\leq
x\\}|,\hskip 8.5359ptx\in\mathbb{R}_{+},$ (3.5)

where $|A|$ denotes the cardinality of a set $A$. Furthermore, the empirical
distribution for the basket and index is defined as

$\hat{F}_{[s]_{M};B}(x)=\frac{1}{M}|\\{m\leq M\colon b^{m}(t)\leq x\\}|,\hskip
8.5359ptx\in\mathbb{R}_{+},$ (3.6)

and respectively

$\hat{F}_{[s]_{M};I}(x)=\frac{1}{M}|\\{m\leq M\colon i^{m}(t)\leq x\\}|,\hskip
8.5359ptx\in\mathbb{R}_{+}.$ (3.7)

We can also generalize the concept of a copula into its empirical counterpart.
The empirical copula counts all the samples that lie below a certain boundary
in the $[0,1]^{N}$ sample space.

$\displaystyle\hat{C}_{[s]_{M}}(u)=\frac{1}{M}\left|\\{m\leq
M\colon\left(F_{S_{1}(t)}(s^{m}_{1}(t)),F_{S_{2}(t)}(s^{m}_{2}(t)),\dots,F_{S_{N}(t)}(s^{m}_{N}(t))\right)\in
A(u)\\}\right|,$ (3.8)

with $u\in[0,1]^{N}$ and
$A(u)=[0,u_{1}]\times[0,u_{2}]\times\dots\times[0,u_{N}]$.

Having defined the empirical copula, it is possible to define the symmetric
property of a sample matrix. We call a matrix $[s]_{M}$ a symmetric sample
matrix if $\hat{C}_{[s]_{M}}(u)=\hat{C}_{[s]_{M}}(P_{\pi}(u))$ for any
permutation $P_{\pi}$.

###  3.2 General Rearrangement Algorithms

Let the time $t\in T_{E}$ be fixed. A rearrangement algorithms are used to
construct multivariate samples of a random vector, for which:

  1. i)

The marginal probability distributions $F_{S_{n}(t)}$ are given,

  2. ii)

The probability distribution $F_{I(t)}$ of the sum
$I(t)=\sum_{n=1}^{N}S_{n}(t)$ is given.

Let $\Sigma$ denote the collection of the above involved (input) probability
distributions. We define a rearrangement algorithm as a map
$\mathbf{A}(M;\Sigma,\omega)=[s]_{M}^{*}$, which produces a sample matrix
$[s]_{M}^{*}$ of size $M$ given the sampling outcome $\omega\in\Omega$. The
output $[s]_{M}^{*}$ of the rearrangement algorithm should have the following
properties: Let $\lVert\cdot\rVert_{\infty}$ be the $L_{\infty}$-norm and let
$\epsilon>0$ be fixed. Then,

$\lVert\hat{F}_{[s]_{M}^{*};n}-F_{S_{n}(t)}\rVert_{\infty}<\epsilon,$ (3.9)

for all $n\leq N$.

$\lVert\hat{F}_{[s]_{M}^{*};I}-F_{I(t)}\rVert_{\infty}<\epsilon.$ (3.10)

In this case, we refer to $[s]_{M}^{*}$ as an _$\epsilon$ -admissible sample
matrix_.

A general rearrangement algorithm $\mathbf{A}$ works based on the following
principle: the algorithm initializes a sample matrix
$[s]^{\text{init}}_{M}=[s]^{\text{init}}_{M}(\omega)$ of the $N$ constituents
given a time $t\in T_{E}$. This is accomplished by drawing samples
independently according to the marginal distribution $F_{S_{n}(t)}$ for each
$n\leq N$. Such a sampling procedure is relatively cheap since no dependency
between the random variables is assumed. Subsequently, the algorithm
rearranges the order of the sample vectors $s^{m}_{n}(t)$ for each constituent
$s_{n}(t)$ to obtain a new sample matrix $[\bar{s}]_{M}$. The rearrangement
affects the sample vector $i(t)$ and therefore the empirical distribution
$\hat{F}_{[\bar{s}]_{M};I}$, while leaving the empirical distribution
$\hat{F}_{[\bar{s}]_{M};n}$ for the constituent $n$ unchanged. The ordering
per constituent $s_{n}(t)$ is changed until we find $[s]_{M}^{*}$ such that
the empirical distribution converges to $F_{I(t)}$, which means that condition
(3.10) is satisfied. Since the empirical marginal distributions
$\hat{F}_{[s]_{M}^{*};n}$ are unchanged, the condition (3.9) is satisfied too.
In the next subsection, we will formalize the procedure of rearrangement.

###  3.3 Rearranging the constituent sample vector

Suppose that we obtain a sample matrix $[s]_{M}^{\text{init}}$, drawn from the
marginal distributions of $S_{1}(t),S_{2}(t),\dots,S_{N}(t)$ for some time
$t\in T_{E}$, independently of each other. For any $n$, the sample vector is
the $M$-dimensional vector

$s_{n}(t)=\left[s_{n}^{1}(t),s_{n}^{2}(t),\dots,s_{n}^{M}(t)\right]^{T}\in\mathbb{R}_{+}^{M}.$
(3.11)

To alter the order of the vector, we will define a permutation on the indices
of the elements. A permutation is defined as a bijective map
$\bar{\pi}\colon\mathbb{N}_{M}\to\mathbb{N}_{M}$ on the set
$\mathbb{N}_{M}=\\{1,2,3,\dots,M\\}$. We can apply the permutation to the
sample order to obtain a new constituent sample vector

$\bar{s}_{n}(t)=\left[s_{n}^{\bar{\pi}(1)}(t),s_{n}^{\bar{\pi}(2)}(t),\dots,s_{n}^{\bar{\pi}(M)}(t)\right]^{T}\in\mathbb{R}_{+}^{M}.$
(3.12)

To simplify the notation of the permutation $\bar{\pi}$ when applying to a
sample vector $s_{n}(t)$, we will use notation
$\pi\colon\mathbb{R}^{m}_{+}\to\mathbb{R}^{m}_{+}$ to denote

$\pi(s_{n}(t))=\left[s_{n}^{\bar{\pi}(1)}(t),s_{n}^{\bar{\pi}(2)}(t),\dots,s_{n}^{\bar{\pi}(M)}(t)\right]^{T}.$
(3.13)

Since we can apply permutations to every constituent sample vector in a sample
matrix, we use $N$ permutations $\pi_{1},\pi_{2},\dots,\pi_{N}$ to obtain a
new sample matrix

$[\bar{s}]_{M}=[\pi_{n}(s_{n}(t)):n\leq N].$ (3.14)

From the initial sample matrix $[s]_{M}^{\text{init}}$, we can define a
finite555The collection is finite since the number of unique permutations on
$\mathbb{N}_{M}$ is finite collection of sample matrices by applying all
possible permutations on the set. We define

$\mathcal{M}=\mathcal{M}([s]_{M}^{\text{init}},\omega)=\\{[\bar{s}]_{M}:\pi_{1},\pi_{2},\dots,\pi_{N}\text{
are permutations on }[s]_{M}^{\text{init}}(\omega)\\},$ (3.15)

with $[\bar{s}]_{M}$ as in Equation 3.14.

######  Example 3.1.

For example, consider the sample matrix $[s]_{M}$ shown in Table 1, where the
second column on the right-hand side is shifted up by 1. This means that
$\bar{\pi}_{2}(x)=x+1\mod M$ and $\bar{\pi}_{n}=id$ for all $n\neq 2$.

| Constituents
---|---
Samples |  | 1 | 2 | … | N
1 | $s_{1}^{1}$ | $s_{2}^{1}$ | $\ldots$ | $s_{N}^{1}$
2 | $s_{1}^{2}$ | $s_{2}^{2}$ | $\ldots$ | $s_{N}^{2}$
⋮ | $\vdots$ | $\vdots$ | $\ddots$ | $\vdots$
M | $s_{1}^{M}$ | $s_{2}^{M}$ | $\ldots$ | $s_{N}^{M}$
(a) Sample matrix $[s]_{M}$

| Constituents
---|---
Samples |  | 1 | 2 | … | N
1 | $s_{1}^{1}$ | $s_{2}^{2}$ | $\ldots$ | $s_{N}^{1}$
2 | $s_{1}^{2}$ | $s_{2}^{3}$ | $\ldots$ | $s_{N}^{2}$
⋮ | $\vdots$ | $\vdots$ | $\ddots$ | $\vdots$
M | $s_{1}^{M}$ | $s_{2}^{1}$ | $\ldots$ | $s_{N}^{M}$
(b) Permutated matrix by shifting column 2

Table 1: Example Permutation: In this permutation, only the second column
$s_{2}(t)$ is shifted. Every sample index increases by 1 for this column. The
other columns are unchanged.

###  3.4 Optimization Formulation

The aim of a rearrangement algorithms $\mathbf{A}$ is to construct an
$\epsilon$-admissible sample. We will now formulate this in an optimization
problem over the possible permutation. First, we define the objective function
$L$ such that

$L([s]_{M})=\lVert\hat{F}_{[s]_{M},I}-F_{I(t)}\rVert_{\infty}=\sup_{x\in\mathbb{R}}\left|\hat{F}_{[s]_{M},I}(x)-F_{I(t)}(x)\right|.$
(3.16)

Consequently, a rearrangement algorithm $\mathbf{A}$ is defined by the
minimization of the objective function $L$ over the set
$\mathcal{M}([s]_{M}^{\text{init}})$:

$\mathbf{A}(M;\Sigma,\omega)=\operatorname*{arg\,min}_{[s]_{M}\in\mathcal{M}([s]_{M}^{\text{init}},\omega)}L([s]_{M})=[s]_{M}^{*}.$
(3.17)

The above formulation aims to find the minimum of $L$ in the permutations.
Previously we formulated the conditions of $[s]_{M}^{*}$ in terms of the
arbitrarily small error $\epsilon$. There is thus a requirement for
convergence of the algorithm, which is justified with the following lemma.
Since the proof of the theorem is tedious, we leave it to the appendix.

######  Lemma 3.3.

Let $\epsilon>0$ be fixed and denote $\mathcal{M}$ as the set of all
permutations of an initial sample matrix $[s]_{M}^{\text{init}}$ of size $M$.
There exists an $M$, such that $\mathcal{M}$ almost surely contains an
$\epsilon$-admissible sample matrix.

###### Proof.

A complete proof is found in C. ∎

The lemma proves that it is possible for the rearrangement algorithm to
converge to a solution and that the error $L$ can be brought arbitrarily low
by increasing the sample size. There is, however, still an important open
question to consider: Are the multiple ways to arrange the samples to obtain
$\epsilon$-admissibility? Remember that from Example 2.1, it is known that
there is no unique copula of a random vector that can match the distribution
function $F_{I(t)}$. It turns out that the algorithm suffers from the same
problem of underdetermination as the copula model from Section 2. It is
possible to obtain two sample matrices $[s]_{M}$ and $[s]_{M}^{\prime}$ with
similar empirical index probability distributions, although the empirical
copulae diverge drastically. This is a direct result of the convergence of
empirical distributions to the cumulative probability distribution [18], and
the observations from the previous section. The implication of this
observation is that the algorithm possibly generates samples for a different
random vector $(S^{\prime}_{1}(t),S^{\prime}_{2}(t),\dots,S^{\prime}_{N}(t))$,
such that

$\sum_{n=1}^{N}S_{n}(t)\stackrel{{\scriptstyle
d}}{{=}}\sum_{n=1}^{N}S^{\prime}_{n}(t)$ (3.18)

The issue of underdetermination can be tackled by adding penalization for
undesired properties to the sample sets in the optimization. In Section 2 we
define a “valid” copula to be admissible and symmetric. Suppose we define a
measure of sample symmetry $\mathcal{S}([s]_{M})>0$666A symmetry measurement
can be based on the empirical copula., such that
$\mathcal{S}([s]_{M})>\mathcal{S}([s]^{\prime}_{M})$ implies that $[s]_{M}$
closer to being symmetric. Then, we define the optimization of the algorithm
as:

$\mathbf{A}(M;\Sigma,\omega)=[s]_{M}^{*}=\operatorname*{arg\,min}_{[s]_{M}\in\mathcal{M}([s]_{M}^{\text{init}})}L([s]_{M})+\frac{\lambda}{\mathcal{S}([s]_{M})},$
(3.19)

where $\lambda>0$ is a penalization coefficient. The rearrangement algorithm
is thus valid if $[s]_{M}^{*}=\mathbf{A}(M,\Sigma,\omega)$ is
$\epsilon$-admissible, i.e. converges to the market-implied probability
distribution, and is symmetric. In this case, $[s]_{M}^{*}$ can be considered
a sample set of a copula $C\in\Theta_{A}^{S}$.

###  3.5 Rearrangement Algorithms for Basket Derivatives

Rearrangement algorithms are well applicable in the context of basket
derivatives. The marginal distributions for the rearrangement algorithm are
the constituent distributions obtained from the local volatility models of
(2.2) and the probability distribution from the sum is the index probability
distribution. In this case, we can apply a rearrangement algorithm for any
$t\in T_{E}$ to obtain samples of $(S_{1}(t),S_{2}(t),\dots,S_{N}(t))$ and
thereby implicitly extracting a copula $C\in\Theta_{A}$. Since we require
samples for all times $t\in T_{E}$, we create a collection of rearrangement
algorithms.

$\mathbf{A}(M,t;\Sigma,\omega)=\\{\mathbf{A}(M;\Sigma(t),\omega):t\in
T_{E}\\}$ (3.20)

Since the samples are created individually for each $t$ they are static,
meaning that there is no time continuity in the samples. As explained in 1,
the pricing of path-dependent requires a calibration of a dynamic model for
$B(t)$ as explained in Figure 3. We can, however, price European call and put
options since their payoff depends only on the final time $t$ at expiry.

######  Theorem 3.1.

Let $[s]_{M}^{*}$ be an $\epsilon$-admissible sample matrix of constituents
for time $t\in T_{E}$ and $\epsilon>0$. Furthermore, suppose that $F_{I(t)}$
and $\hat{F}_{[s]_{M}^{*};I}$ agree outside an interval
$[a,b]\subset\mathbb{R}^{+}$. Then,

$|V_{\text{Call}}^{MC}([s]_{M}^{*};K,t)-V_{\text{Call}}^{\text{Market}}(K,t)|<\epsilon(b-a),$
(3.21)

where $V_{\text{Call}}^{MC}([s]_{M}^{*};K,t)$ is the Monte Carlo price of a
European call option given the sample matrix $[s]_{M}^{*}$.

###### Proof.

For simplicity, we will assume a constant interest rate of $r=0$. In the case
of non-zero interest rates, a constant discounting factor has to be included.
Since $\hat{F}_{[s]_{M}^{*};I}$ is not differentiable, we use Riemann-
Stieltjes integrals to express the pricing integrals. The market price of a
call option with strike $K$ expiring at $t$ can be written as

$V^{\text{Market}}_{\text{Call}}(t,K)=\int_{-\infty}^{\infty}(x-K)^{+}{\rm
d}F_{I(t)}(x),$

and the Monte-Carlo price as

$V_{\text{Call}}^{MC}([s]_{M}^{*};K,t)=\sum_{j\leq
M}\frac{\left(i^{j}(t)-K\right)^{+}}{M}=\int_{-\infty}^{\infty}(x-K)^{+}{\rm
d}\hat{F}_{[s]_{M}^{*};I}(x),$

where both integrals are Riemann-Stieltjes integrals. By the linearity of the
Riemann-Stieltjes integrals, we have

$\displaystyle|V_{\text{Call}}^{MC}([s]_{M}^{*};K,t)-V_{\text{Call}}^{\text{Market}}(K,t)|$
$\displaystyle=\left|\int_{-\infty}^{\infty}(x-K)^{+}{\rm
d}\hat{F}_{[s]_{M}^{*};I}(x)-\int_{-\infty}^{\infty}(x-K)^{+}{\rm
d}F_{I(t)}(x)\right|$
$\displaystyle=\left|\int_{-\infty}^{\infty}(x-K)^{+}{\rm
d}(\hat{F}_{[s]_{M}^{*};I}-F_{I(t)})(x)\right|$
$\displaystyle=\left|\int_{a}^{b}(x-K)^{+}{\rm
d}(\hat{F}_{[s]_{M}^{*};I}-F_{I(t)})(x)\right|.$

Since $\hat{F}_{[s]_{M}^{*};I}-F_{I(t)}$ is zero at the boundaries,
integration by parts yields

$\displaystyle\left|\int_{a}^{b}(x-K)^{+}{\rm
d}(\hat{F}_{[s]_{M}^{*};I}-F_{I(t)})(x)\right|$
$\displaystyle\leq\left|\int_{a}^{b}\hat{F}_{[s]_{M}^{*};I}(x)-F_{I(t)}(x){\rm
d}x\right|$
$\displaystyle\leq\int_{a}^{b}\left|\hat{F}_{[s]_{M}^{*};I}(x)-F_{I(t)}(x)\right|{\rm
d}x$ $\displaystyle<\epsilon(b-a).$

The proof for put options works the same way. ∎

##  4 Example: Iterative Sort-Mix (ISM) algorithm & Implementation Details

We will now present an implementation of a rearrangement algorithm
$\mathbf{A}(M,t;\Sigma,\omega)$ for basket option pricing as defined in the
previous section. For any time $t\in T_{E}$, the algorithm efficiently creates
a sample matrix $[s]_{M}$ which is used to estimate the probability density
function $f_{B(t)}$ and therefore allows us to price basket options. Note
that, contrary to the previous section, the algorithm is not an optimization,
and thus the pricing error cannot be reduced to an arbitrarily low value. The
algorithm is nevertheless able to provide sufficiently accurate results. In
the empirical section, we will demonstrate that “sufficient” means that the
algorithm performs well enough on empirical data to rapidly and accurately
price basket options.

As is usual for a rearrangement algorithm, the first step of the algorithm is
to generate the initial sample matrix $[s]_{M}^{\text{init}}$. This step is
generic, yet a description will be provided in Section 4.2. The core part of
the algorithm is the rearrangement of the initial matrix, which we will first
explain in a heuristic way. In the following subsections, a detailed procedure
is laid out and a pseudo-code example is provided in A.

###  4.1 Heuristic Explanation

The core principle of the algorithm is not an optimization as suggested
Equation 3.19. Rather, the algorithm iteratively selects and stores a subset
of the samples that are considered suitable, until all samples are selected
out. To see how this concretely works, we first notice that there are two
general ways to arrange the individual vectors $s_{n}(t)$: We can sort the
array from the smallest sample to the largest sample. In this case, we have

$s_{n}^{m}(t)\leq s_{n}^{l}(t)\iff m\leq l,\hskip 14.22636pt\forall n\leq N.$
(4.1)

Alternatively, we can use a random number generator to arbitrarily mix the
individual sample vectors. This is equivalent to applying permutations
$\pi_{n}$ at random and obtaining $\pi(s_{n}(t))$ for each constituent $n$.

We now also note that in the first arrangement, when all vectors are sorted,
the sample correlation between the individual vectors is high. For any pair
$n,l\leq N$, we have that the Pearson sample correlation coefficient between
$s_{n}(t)$ and $s_{l}(t)$ is close to 1. On the other hand, when all sample
vectors are randomly ordered, the correlation is close to 0. We now arrange
each column in the sample matrix $[s]_{M}^{\text{init}}$ in the two
arrangements and each time plot the resulting empirical PDF of the index as a
histogram. This is given in Figure 4, where we also plot the PDF of the
“target PDF” which is the market-observed PDF of the index.

Figure 4: Index probability density function: The graph shows the empirical
density of two sample sets of an index as a histogram. The yellow histogram
shows the samples when they are in the “mixed” arrangement. The brown
histogram shows the “sorted” arrangement, displaying the fat tails due to the
strong correlation of the constituents - Samples based on data from DJIA 1
Year

The rationale behind the algorithm is straightforward: We want to create a
superposition of sorted samples and mixed samples, such that the empirical
probability density function matches the target PDF. To do so, we iteratively
arrange the samples in the two arrangements, and store the samples which “fit
under the target PDF”. Fitting under the target PDF means that the resulting
empirical density function is locally below the target PDF. We then continue
the iteration only with the samples that are not stored and choose the other
arrangement. The algorithm always starts with the “sort” arrangement.

After only a few iterations, the algorithm can no longer make progress, and
the remaining samples are added to the previously stored ones. This results in
a sub-optimal empirical distribution compared to the target. As we will see,
this error is not significant, and the function $f_{B(t)}$ can still be
estimated effectively.

The resulting sample matrix $[s]_{M}^{*}$ is thus $\epsilon$-admissible up to
the error $\epsilon$ of the remaining samples which cannot be fit.
Furthermore, the samples in both sorted and mixed form are symmetric, which
means that the resulting superposition of samples is also symmetric.

Figure 5 shows a graphical example of one iteration of the algorithm 777the
data is from Section 5.

Figure 5: One iteration of sorting - mixing procedure: The algorithm starts
with 10’000 samples. Step 1 (Top Left): Samples are ordered. Step 2 (Top
Right): All valid samples (below target PDF) are removed from the sample
matrix and stored away. Step 3 (Bottom Left): The target is updated, and the
remaining samples are outside of the target range. Step 4 (Bottom Right): The
remaining samples are mixed together and are now ready to be removed again.
Only $\sim$ 400 samples are left after this iteration.

###  4.2 Independent Constituent Sampling

The initial step consists of generating $M\in\mathbb{N}$ samples of the
constituent prices according to their marginal distribution, to initiate the
matrix $[s]_{M}^{\text{init}}$. This is achieved by generating uniform random
variables and composing them with the inverse of the cumulative distribution
function of the constituents. Let $t>0$ be some expiry time, and let $M$ be
the number of required samples. Given the probability distributions
$F_{S_{n}(t)}$ of $S_{n}(t)$, which we obtain from calibrating Equation 2.2 to
the constituents’ vanilla option prices, we generate the samples separately
for each constituent. We will use the following mapping:

$S_{n}(t)\stackrel{{\scriptstyle d}}{{=}}F_{S_{n}(t)}^{-1}\left(U_{n}\right),$
(4.2)

where $U_{n},1\leq n\leq N$ are distinct and independent uniformly distributed
random variables. Therefore, we obtain the samples $s_{n}(t)$ as

$s_{n}(t)=\left[s^{1}_{n}(t),s^{2}_{n}(t),\dots,s^{M}_{n}(t)\right]^{T}=\left[F^{-1}_{S_{n}(t)}(u_{n}^{1}),F^{-1}_{S_{n}(t)}(u^{2}_{n}),\dots,F^{-1}_{S_{n}(t)}(u^{M}_{n})\right]^{T},$
(4.3)

where $u_{n}^{1},u_{n}^{1},\dots,u_{n}^{M}$ are the uniform samples.

###  4.3 Range Discretization and Target Vector

Before the algorithm can arrange the samples, we require a discretization of
the involved PDFs. We aim to obtain a target vector $V\in\mathbb{N}^{k}$ with
$k\in\mathbb{N}$, which specifies the expected number of samples in a subset
of the range of $I(t)$ under the target PDF $f_{I(t)}$.

We begin by defining a discretization of the range of the index. Since the
range of $I(t)$ is $\mathbb{R}_{+}$ and therefore unbounded, we need to reduce
it to an appropriately bounded set $[g_{0},g_{K}]\subset\mathbb{R}_{+}$ to
ensure numerical efficiency. We choose the bounds $g_{0}<g_{K}$ for the
interval with the following interpretation: For $M$ samples from an admissible
copula, we expect exactly one sample below $g_{0}$ and one sample above
$g_{k}$. These values are given by the probability distribution function

$g_{0}=F^{-1}_{I(t)}\left(\frac{1}{M}\right)\hskip 14.22636pt\text{ and
}\hskip 14.22636ptg_{K}=F^{-1}_{I(t)}\left(\frac{M-1}{M}\right).$ (4.4)

We then define an equidistant grid $g_{0}<g_{1}<g_{2}<\dots<g_{K}$, which
define $K$ bins $G_{k}=(g_{k-1},g_{k}]$888For completeness, we extend $G_{0}$
with $g_{0}$ to $[g_{0},g_{1}]$ to partition the interval 999The parameter $K$
is a factor for the quality and computational complexity of the algorithm and
should be chosen in relation to $M$, for instance $\frac{M}{K}=10$. We thus
have

$[g_{0},g_{K}]=\bigcup_{0\leq k\leq K}G_{k}.$ (4.5)

For each bin, we can now calculate the expected samples to fall into each bin
under $F_{I(t)}$. This value is given by the cumulative distribution, rounded
to an integer

$\Delta F_{k}=\big{\lceil}M\;\mathbb{Q}(I(t)\in
G_{k})\big{\rfloor}=\big{\lceil}M\left(F_{I(t)}(g_{k})-F_{I(t)}(g_{k-1})\right)\big{\rfloor}\in\mathbb{N},$
(4.6)

where $\lceil x\rfloor$ denotes $x$ rounded to its nearest integer. We capture
the discrete PDF with the _target vector_

$V=\left[\Delta F_{1},\Delta F_{2},\dots,\Delta
F_{K}\right]\in\mathbb{N}^{K}.$ (4.7)

Furthermore, we introduce a discretization of the loss function on the PDF,
given by $\hat{\ell}$. For a sample matrix $[s]_{M}\in\mathcal{M}$, we define
the bin count $c([s]_{M})$ with

$c_{k}([s]_{M})=|\\{1\leq m\leq M:i^{m}\in G_{k}\\}|,$ (4.8)

with $1\leq k\leq K$. The discrete error is now defined as the difference
between expected and actual discrete PDF:

$\hat{\ell}([s]_{M})=\frac{\sum_{k=1}^{K}|c_{k}([s]_{M})-\Delta
F_{k}|}{2M}\in[0,1].$ (4.9)

We divide by two since a sample that cannot be matched to a bin will be
counted once twice in the numerator.

###  4.4 Sorting - Mixing Procedure

With the target vector $V$ defined and samples $s_{n}(t)$ initiated, we will
now describe the details of the core of the algorithm, where the different
sample arrangements are iterated. The main algorithm thus runs through
iterations of first sorting and then mixing samples. Suppose we have a sample
matrix $[s]_{M}$ with constituent sample vectors $s_{n}(t)$ for all $n\leq N$
and index vector $i(t)=\sum_{n=1}^{N}s_{n}(t)$. An iteration of the sorting-
mixing procedure starts by arranging the constituent sample vectors in sorted
order and obtaining the index sample vector $i(t)$. Now, for every bin
$G_{k}$, we collect the vector positions $m$ 101010By ”vector positions” we
mean the index $m$ of $i^{m}(t)$ of the vector $i(t)$. We refrain from calling
it the “index” of the vector due to the ambiguity with the index $I(t)$.  of
the samples that have values $i^{m}(t)$ in $G_{k}$. We define the set of
positions of the samples in $G_{k}$ as

$R_{k}=\\{1\leq m\leq M:i^{m}(t)\in G_{k}\\}.$ (4.10)

The sample count $c_{k}$, which was defined previously, is thus given by
$c_{k}=|R_{k}|$. We want to determine a set of positions $\bar{R_{k}}\subset
R_{k}$, which we can store away and remove from the sample matrix. The count
$c_{k}$ does not match the expected count, given by $\Delta F_{k}$. We
therefore want to remove _at most_ $\Delta F_{k}$ samples and we might have to
reduce the set $R_{k}$ of samples to remove. If $c_{k}>\Delta F_{k}$, we need
to reduce the set $R_{k}$ to a smaller set by $\bar{R}_{k}$. We do this by
randomly selecting $\Delta F_{k}$ elements from $R_{k}$. If $c_{k}\leq\Delta
F_{k}$, all the samples from $R_{k}$ can be removed and we have
$\bar{R}_{k}=R_{k}$. In both cases, we now have

$\bar{c}_{k}=|\bar{R}_{k}|=\min(c_{k},\Delta F_{k}).$ (4.11)

To visualize the involved quantities, we see that the height of the brown
histogram in Figure 5 is given by $c_{k}$, while the height of the yellow
columns is $\bar{c}_{k}$.

We can now store all samples with indices in $\bar{R}_{k}$ and reduce the
sample vectors. The updated constituent sample vectors $s_{n}(t)^{\prime}$
contain all elements $s_{n}^{m}(t)$ such that $m\notin\bigcup_{k\leq
K}\bar{R_{k}}$. The length of $s_{n}(t)^{\prime}$ is thus $M-\sum_{k\leq
K}\bar{c}_{k}$ for all constituents.

Lastly, the target vector needs to be updated too, since we no longer require
that many samples per bin. Since we remove $\bar{c}_{k}$ samples, we update as

$\Delta F_{k}^{\prime}=\Delta F_{k}-\bar{c}_{k}.$ (4.12)

To finish one iteration of the algorithm, we need to repeat this procedure in
a random arrangement. We thus now “shuffle” the elements of
$s_{n}(t)^{\prime}$ and then repeat the derivation of $R_{k},\bar{R}_{k}$, et
cetera, removing and storing more samples. After sorting and shuffling once,
we conclude one iteration. Figure 6 shows how the final PDF of the samples
matches the target PDF.

Figure 6: Final PDF after 5 iterations. The match of the empirical PDF and
target PDF means that options are priced correctly. This can be seen on the
implied volatility graph on the right (Example DJIV 1Y)

######  Example 4.1.

We demonstrate the workings of the algorithm on a small dataset with 2 assets
and 10 samples. The samples are marginally uniformly distributed on $[0,1]$.
The index samples fall into $[0,2]$, which we split into $K=3$ equally large
bins with boundaries

$\left\\{0,\frac{2}{3},\frac{4}{3},2\right\\}.$

The bins are setup such that111111Since the $[0,2]$ is already bounded, there
is no need to obtain $g_{0}$ and $g_{K}$ by the quantiles.

$G_{1}=\left[0,\frac{2}{3}\right],\quad
G_{2}=\left(\frac{2}{3},\frac{4}{3}\right],\quad
G_{3}=\left(\frac{4}{3},2\right].$

Furthermore, suppose we obtain the from the index option quotes the following
discrete probability density function for $I$:

$f_{I}(x)=\begin{cases}0.3,\text{ if }x\in G_{1},\\\ 0.5,\text{ if }x\in
G_{2},\\\ 0.2,\text{ if }x\in G_{3}.\end{cases}$

We then calculate the target vector for $M=10$ samples as

$\Delta F_{k}=\lceil M\cdot\mathbb{Q}(I\in G_{k})\rfloor,\quad
k\in\\{1,2,3\\}$

from which we obtain

$V=\left[\Delta F_{1},\Delta F_{2},\Delta F_{3}\right]=\left[3,5,2\right].$

Algorithm: We start the algorithm by creating samples for $S_{1}$ and $S_{2}$
independently. This step is shown in a) in Table 2. Then, in b) we sort both
column vectors smallest to largest and sum the columns to obtain the index
$I$. For this arrangement we obtain

$c=[c_{1},c_{2},c_{3}]=[4,3,3],$

which means that the amount of samples to be removed are

$\bar{c}=\min(c,\Delta F)=[3,3,2].$

We therefore select three samples from bin 1, three samples from bin 2 and two
samples from bin 3, which is shown in c), indicated by the green color that a
sample is picked out. For bin 1 and 3 we select them at random since
$c_{1/3}>\Delta F_{1/3}$. Since two samples cannot be selected, the loss at
this stage is given by $\ell([s]^{*}_{M})=\frac{2}{10}$. We update the target
vectors according to $\Delta F_{k}^{\prime}=\Delta F_{k}-\bar{c}_{k}$ to
obtain

$\Delta F_{k}^{\prime}=[0,2,0].$

The remaining samples are now entering the mixing stage. In d), we mix the
samples for $S_{1}$ and $S_{2}$ column-wise and aggregate to obtain an new
value for $I$. Based on these values we obtain the new $c$ as $c=[0,2,0]$.
Since this is exactly the remaining target vector, we select both new samples
and the algorithm is finished. The remaining error is thus
$\ell([s]_{M}^{*})=0$ and the selected samples can be used.

$S_{1}$ | $S_{2}$
---|---
0.29 | 0.31
0.14 | 0.47
0.26 | 0.17
0.44 | 0.07
0.05 | 0.01
1.00 | 0.69
0.31 | 0.83
0.76 | 0.49
0.72 | 0.41
0.04 | 0.76
(a) Initial data sampled as two independent uniform random variables

| $S_{1}$ | $S_{2}$ | $I$
---|---|---|---
| 0.04 | 0.01 | 0.05
| 0.05 | 0.07 | 0.12
| 0.14 | 0.17 | 0.31
$G_{1}$ | 0.26 | 0.31 | 0.57
| 0.29 | 0.41 | 0.70
| 0.31 | 0.47 | 0.78
$G_{2}$ | 0.44 | 0.69 | 1.13
| 0.72 | 0.76 | 1.48
| 0.76 | 0.83 | 1.59
$G_{3}$ | 1.00 | 0.76 | 1.76
(b) Sort data column-wise, compute $I=S_{1}+S_{2}$ and split into bins based
on index value

| $S_{1}$ | $S_{2}$ | $I$
---|---|---|---
| 0.04 | 0.01 | 0.05
| 0.05 | 0.07 | 0.12
| 0.14 | 0.17 | 0.31
$G_{1}$ | 0.26 | 0.31 | 0.57
| 0.29 | 0.41 | 0.70
| 0.31 | 0.47 | 0.78
$G_{2}$ | 0.44 | 0.69 | 1.13
| 0.72 | 0.76 | 1.48
| 0.76 | 0.83 | 1.59
$G_{3}$ | 1.00 | 0.76 | 1.76
(c) Select $\bar{c}_{k}$ samples per $G_{k}$. Two samples cannot be attributed
to the right bin.

| $S_{1}$ | $S_{2}$ | $I$
---|---|---|---
| 0.14 | 0.83 | 0.97
$G_{2}$ | 0.76 | 0.17 | 0.93
(d) The remaining two samples are mixed and binned.

Table 2: One iteration on toy data

##  5 Numerical Results

###  5.1 Setup

It remains to demonstrate the effectiveness of the ISM algorithm on actual
market data. In this section, we will consider an experiment on actual market
data for the Dow Jones Industrial Average (DJIA) and its 30 constituents as of
12th August 2021121212An accompanying Python repository of the implementation
is available at https://github.com/NFZaugg/BasketOptionsRearrangement. The
market data consists of implied volatility surfaces and spot prices for the 30
assets and the index. Table 3 shows an overview of the available data and the
selection for the parameters.

To extract the marginal distribution of the constituents as well the index
distribution we calibrate SABR (Stochastic Alpha-Beta-Rho) parameters [9] for
each of the 30 constituents and each maturity $T$, fixing parameter $\beta$ to
$0.9$. The choice of SABR is not relevant as long as it provides a smooth
price of $S_{n}(t)$, such that we can estimate $F_{S_{n}}(t)$. The first few
of the calibrated parameters are shown in Table 4. The same model is
calibrated for the index options which enables the estimation of
$F_{I^{\text{Mkt}}(t)}$.

Data Information
---
Date | 2021/08/12
Number of constituents $(N)$ | 30
Constituent weights | 0.066
Vol Surface Moneyness | 0.8, 0.85, 0.9, 0.95, 0.975, 1, 1.025, 1.05, 1.1, 1.15, 1.2
Vol Surface Maturity | 3m, 6m, 1y, 1.25y, 2y

Hyperparameters
---
Number of Samples $(M)$ | 20’000
Number of Bins | 1’400
Number of Iterations | 10
Constituents Marginal Model | SABR, fixed $\beta=0.9$
Index model | SABR, fixed $\beta=0.9$

Table 3: Experiment setup

We generate a matrix of $M$ samples of $N$ independent uniform random
variables and use Equation 4.2 to transpose those into random variables with
the proper marginal distribution function.

With the distribution of the index, we calculate the target vector (4.7) with
2’000 equidistant bins, yielding the discrete PDF we are aiming to match.
Having collected all the input, we run the algorithm to rearrange the samples.

Name | $\beta$ | $\alpha$ | $\rho$ | $\gamma$
---|---|---|---|---
United Health | 0.9 | 0.47 | -0.52 | 1.20
Home Depot | 0.9 | 0.45 | -0.28 | 1.47
Goldman Sachs | 0.9 | 0.50 | -0.33 | 1.46
Microsoft Corp | 0.9 | 0.43 | -0.43 | 1.49
… | … |  | … |
DJIA | 0.9 | 0.31 | -0.64 | 1.91
Table 4: First 4 calibrated parameters, $T=3m$

###  5.2 Index Repricing

Figure 7 shows the implied volatility of the index option market prices vs the
prices obtained by the model. The model can reprice the index option from the
market well. The remaining error is split between calibration inaccuracies of
$F_{I(t)}$ and the discrete error, which is displayed in Table 5. The errors,
which do not exceed 0.2% in the implied volatilities per moneyness, are well
within the bid and ask quotes of the options (ranging from 0.2% to 4%). We
conclude that the model is fit to price basket options according to the
market-implied correlation skew.

Figure 7: Market vs model implied volatilities - DJIA options. Note: The increased mispricing of 1.5y and 2y mainly stem from numerical inaccuracies in $F_{I(t)}$, which assumes significant mass at $I(t)=0$ Maturity | 3m | 6m | 1y | 1.5y | 2y
---|---|---|---|---|---
Discrete Error $\hat{\ell}(S)$ | $2.09\%$ | $1.54\%$ | $1.6\%$ | $1.92\%$ | $2.27\%$
Table 5: Discrete error per maturity

###  5.3 Performance

Given the results from the previous subsection, we analyze the time spent on
each step in the calibration of the model. Table 6 shows that for each
maturity, the total time spent was between 5 and 7 seconds. Most of the time
was spent on calibrating constituents and generating independent samples. The
arranging of the samples took about 1/3 of the time. Given that this algorithm
can be run in parallel, the LVM can thus be calibrated in below 7 seconds with
20’000 samples.

Action | 3m | 6m | 1y | 1.5y | 2y
---|---|---|---|---|---
Calibrating Constituents | 1.76 | 1.60 | 1.76 | 1.83 | 2.09
Calibrating Index (Bid Ask Mid) | 0.19 | 0.19 | 0.2 | 0.2 | 0.19
Sampling Constituents | 1.44 | 1.52 | 1.63 | 1.56 | 1.79
Preparing Target | 0.07 | 0.08 | 0.08 | 0.07 | 0.08
Arranging Samples | 1.27 | 1.56 | 1.72 | 1.92 | 2.1
Total | 5.73s | 5.95s | 6.39s | 6.58s | 6.25s
Table 6: Time Breakdown (s) for example of 30 constituents and 20’000 samples

###  5.4 Greeks Calculation

The rapid calibration of the model allows for fast Greeks calculation for the
basket derivatives with a finite difference approximation. Suppose that
$\Theta_{LVM}=\Theta_{LVM}(\sigma,\sigma_{I})$ is the local variance function
of (1.4) after calibration given the input market data $\sigma$ and
$\sigma_{I}$. Here, $\sigma=(\sigma_{1},\sigma_{2},\dots,\sigma_{N})$ and
$\sigma_{I}$ are the implied volatility curves of the constituent and index,
respectively.

We express the price of a derivative
$V\left(B(0);\Theta_{LVM}\right)=V\left(B(0)\right)$ as a pricing function of
a basket derivative given the parameters and a spot price $B(0)$. The delta
and gamma of $V$, the first and second derivatives with respect to the spot
price, are approximated using finite differences, for instance

$\displaystyle\Delta_{B}=\frac{\partial V\left(B(0)\right)}{\partial B(0)}$
$\displaystyle\approx\frac{V\left(B(0)+\epsilon\right)-V\left(B(0)\right)}{\epsilon},$
(5.1) $\displaystyle\gamma_{B}=\frac{\partial^{2}V\left(B(0)\right)}{\partial
B(0)^{2}}$
$\displaystyle\approx\frac{V\left(B(0)-\epsilon\right)-2V\left(B(0)\right)+V\left(B(0)+\epsilon\right)}{\epsilon^{2}},$
(5.2)

for some bump size $\epsilon>0$. The sensitivity to the individual
constituents is given by the chain rule 131313Note that if the basket is
weighted, the delta needs to be weighted too

$\displaystyle\frac{\partial V\left(B(0)\right)}{\partial S_{n}(0)}$
$\displaystyle=\Delta_{B}\cdot\frac{\partial B(0)}{\partial
S_{n}(0)}=\Delta_{B}.$ (5.3)

Furthermore, one can also calculate the sensitivities to the constituent
implied volatility curves. Let vega $v_{n}$ be the sensitivity of the basket
derivative to the $n$-th constituent implied volatility curve,

$\displaystyle v_{n}=\frac{\partial
V\left(B(0)\right)}{\partial\sigma_{n}}\approx\frac{V\left(B(0);\Theta_{LVM}(\sigma+\mathbb{\epsilon}_{n},\sigma_{I})\right)-V\left(B(0);\Theta_{LVM}(\sigma,\sigma_{I})\right)}{\epsilon},$
(5.4)

where $\mathbb{\epsilon}_{n}$ is the vector
$(0,0,\dots,0,\epsilon,0,\dots,0)$. To obtain
$V\left(B(0);\Theta_{LVM}(\sigma_{n}+\epsilon,\sigma_{I})\right)$ we need to
re-calibrate the model, as $\Theta_{LVM}$ are changing due to the shift in
implied volatility. Fortunately, the re-calibration is a surprisingly
computationally cheap procedure. As a first step, we re-calibrate the
volatility model for the constituent $n$. This provides us with an updated
probability distribution function

$F_{S_{n}(t)}^{\prime}(x)=\mathbb{Q}\left(S_{n}(t;\sigma_{n}+\epsilon)<x)\right),\;x\in\mathbb{R},$
(5.5)

where $S_{n}(t;\sigma_{n}+\epsilon)$ indicates the constituent price process
at $t$ under its new distribution due to the shift in volatility. To maintain
the same dependency structure as before the shift, we have to adjust the
samples consistently. In the initial calibration, for a fixed $t\in T_{E}$, we
obtained a sample vector for a uniform random variable
$(u_{n}^{1},u_{n}^{1},\dots,u_{n}^{M})$. Furthermore, we used a permutation
$\pi_{n}$ to rearrange the samples. We can now simply obtain the samples with
the new distribution as

$\\{{F_{S_{n}(t)}^{\prime}}^{-1}(u_{n}^{m}):m\leq M)\\},$

and order them in the same way as the old samples

$\\{{F_{S_{n}(t)}^{\prime}}^{-1}(\pi_{n}(u_{n}^{m})):m\leq M\\}.$

Important is that we reuse the same samples for the uniform random variable
$u_{n}^{m}$, as the permutation is determined based on these values.
Alternatively, the same can be obtained by applying the function

$x\mapsto{F_{S_{n}(t)}^{\prime}}^{-1}\left(F_{S_{n}(t)}\left(x\right)\right),$
(5.6)

to the pre-shift samples of $S_{n}(t)$. We thus only adjust the samples for
the $n$-th constituent and re-calibrate the local variance function for every
$t\in T_{E}$. The total computational cost is low as no new sample ordering is
required. Results for $\Delta,\gamma,v$ are shown in Table 7, Table 8. Other
Greeks are obtained similarly.

Strike (relative to ATM) | 0.8 | 0.9 | 1 | 1.1 | 1.2
---|---|---|---|---|---
$\Delta$ | 0.95 | 0.88 | 0.60 | 0.09 | 0.01
$\gamma$ | 1e-7 | 1e-6 | 0.03 | 1e-10 | 1e-10
Table 7: $\Delta,\gamma$ of DJIA 3m European option at spot price for various strikes  Constituent Name | $v_{i}$
---|---
United Health | 4.63%
Home Depot | 3.68%
Goldman Sachs | 3.64%
Microsoft Corp | 3.02%
Salesforce Inc | 2.32 %
Table 8: Vega of DJIA 3m European option for first 5 constituents at the money
(Upwards finite difference, $\epsilon=0.01$ (1% in implied volatility))

##  6 Conclusions

Pricing basket options consistently with index options poses a significant
challenge. Creating a model that accurately reproduces the index skew is often
challenging to calibrate efficiently and, therefore, undesirable for practical
purposes. The difficulty arises from the intricate nature of modeling the
dependency structure between constituents, which is crucial for accurate
pricing and is often not straightforward to extract from available market
data.

Copula functions, describing the joint distribution of random variables mapped
to a uniform space, emerge as valuable tools for modeling such dependency
structures. Calibrating a copula model to the available market data remains
difficult, due to two main issues. Firstly, the inherent multidimensionality
of the problem makes calibration of a multidimensional model computationally
expensive. Furthermore, the problem of underdetermination of the copula $C$
can lead to errors in the pricing of basket options, even if the model can
replicate index options without errors. In this paper, we define a copula $C$
to be a valid pricing model for basket options if the copula is symmetric and
replicates the index PDF.

Recognizing the challenges of directly modeling copulas, we utilize
rearrangement algorithm to imply an admissible and symmetric copula for the
dependency of the constituents. These algorithms work by initializing samples
according to the known marginal distribution and then pragmatically
rearranging them to reflect the dependency structure to match the index
probability distribution. Extra constraints based on the empirical symmetry
are added to infer a symmetric copula.

We develop a particular rearrangement algorithm, called the Iterative-Sort-Mix
algorithm, that automatically calibrates the local volatility function using a
simple rearrangement technique, sidestepping the need for an expensive
multivariate parameter optimization. Furthermore, the algorithm calculates the
Greeks, in particular sensitivities to the constituents’ implied vol curves,
with ease and avoids additional expensive calculations.

The viability of this approach is tested on historical market data, revealing
its effectiveness in calibrating a model for pricing baskets with up to 30
constituents within seconds and with limited errors. The algorithm thus
exhibits good performance and is suitable for practical purposes.

Furthermore, based on the theoretical considerations, the approach encourages
the development of novel and enhanced rearrangement algorithms with robust
convergence properties. A promising direction is the use of artificial
intelligence to address the high-dimensional problem of rearranging the
samples in an appropriate order.

## References

  * Andreasen and Huge [2011] J. Andreasen and B. Huge.  Volatility interpolation.  _Risk_ , 24(3):76, 2011.
  * Bedendo et al. [2010] M. Bedendo, F. Campolongo, E. Joossens, and F. Saita.  Pricing multiasset equity options: How relevant is the dependence function?  _Journal of Banking & Finance_, 34(4):788–801, 2010.
  * Bernard et al. [2021] C. Bernard, O. Bondarenko, and S. Vanduffel.  A model-free approach to multivariate option pricing.  _Review of Derivatives Research_ , 24:135–155, 2021.
  * Branger and Schlag [2004] N. Branger and C. Schlag.  Why is the index smile so steep?  _Review of Finance_ , 8(1):109–127, 2004.
  * Embrechts et al. [2013] P. Embrechts, G. Puccetti, and L. Rüschendorf.  Model uncertainty and var aggregation.  _Journal of Banking & Finance_, 37(8):2750–2764, 2013.
  * Euronext [2023] Euronext.  Euronext fact book 2022, 2023.  URL https://live.euronext.com/en/resources/statistics.
  * Grzelak et al. [2023] L. A. Grzelak, J. Jablecki, and D. Gatarek.  Efficient pricing and calibration of high-dimensional basket options.  _International Journal of Computer Mathematics_ , 2023.  doi: 10.1080/00207160.2023.2266051.
  * Guyon [2016] J. Guyon.  Calibration of local correlation models to basket smiles.  _The Journal of Computational Finance_ , 20(2), 2016.  doi: 10.21314/jcf.2016.326.
  * Hagan et al. [2002] P. S. Hagan, D. Kumar, A. Lesniewski, and D. E. Woodward.  Managing smile risk.  _Wilmott Magazine_ , September:84–108, 2002.
  * Koster and Oeltz [2019] F. Koster and D. Oeltz.  A pairwise local correlation model.  _Journal of Computational Finance_ , 2019.
  * Langnau [2010] A. Langnau.  A dynamic model for correlation.  _Risk_ , 23(4):74, 2010.
  * Lucic [2012] V. Lucic.  Correlation skew via product copula.  In _Financial Engineering Workshop, Cass Business School_ , 2012.
  * Oosterlee and Grzelak [2019] C. W. Oosterlee and L. A. Grzelak.  _Mathematical Modeling and Computation in Finance: With exercises and Python and MATLAB computer codes_.  World Scientific Publishing Europe Ltd, 2019.
  * Pirjol [2023] D. Pirjol.  Smile-consistent basket skew.  _Risk_ , 2023.
  * Qu [2010] D. Qu.  Pricing basket options with skew.  _Wilmott magazine_ , pages 58–64, 2010.
  * Reghai [2010] A. Reghai.  Breaking correlation breaks.  _Risk_ , 23(10):92, 2010.
  * Salmon and Schleicher [2007] M. H. Salmon and C. Schleicher.  Pricing multivariate currency options with copulas.  In _Copulas: From Theory to Application in Finance_. Risk Books, 2007.
  * Tucker [1959] H. G. Tucker.  A generalization of the Glivenko-Cantelli theorem.  _The Annals of Mathematical Statistics_ , 30(3):828–830, 1959.
  * Van den Goorbergh et al. [2005] R. W. Van den Goorbergh, C. Genest, and B. J. Werker.  Bivariate option pricing using dynamic copula models.  _Insurance: Mathematics and Economics_ , 37(1):101–114, 2005.
  * Xu and Zheng [2010] G. Xu and H. Zheng.  Basket options valuation for a local volatility jump–diffusion model with the asymptotic expansion method.  _Insurance: Mathematics and Economics_ , 47(3):415–422, 2010.
  * Zetocha [2015] V. Zetocha.  Correlation skew via stochastic correlation and jumps.  _Risk_ , 28(12):1–6, 2015.

##  Appendix A Algorithm Pseudo-Code

We express the algorithm as a pseudo-code:

/* As an input, we have initialized samples, the discretized target vector,
and the corresponding bins of the discretization */

Input: sampleMatrix: $(M\times N)$, targetVector $(g\times 1)$, bins
$((g+1)\times 1)$

1

Output: outputMatrix: $(M\times N)$

2 Def  _Main(_sampleMatrix,targetVector,bins_)_:

3 remainingMatrix $\leftarrow$ sampleMatrix;

4 outputMatrix = new emptyMatrix;

5 while  _$\left(\frac{remainingMatrix.Dim1}{M}\geq\epsilon\right)$_ do

/* Part 1: Sort, pick, update */

6 Sort(remainingMatrix);

7 outputMatrix $\leftarrow$ SelectValidSamples(outputMatrix, remainingMatrix,
targetVector, bins);

/* Part 2: Mix, pick, update */

8 Mix(remainingMatrix);

9 outputMatrix $\leftarrow$ SelectValidSamples(outputMatrix, remainingMatrix,
targetVector, bins);

10

11 end while

12 return outputMatrix;

13

14Def  _Sort(_matrix_)_:

/* This method sorts each column of the matrix in ascending order */

15 for  _column in matrix_ do

16 column $\leftarrow$ sortVectorAscending(column)

17 end for

18 return matrix;

19

20

21Def  _Mix(_matrix_)_:

/* This method mixes each column of the matrix */

22 for  _column in matrix_ do

23 column $\leftarrow$ mixColumnRandomly(column)

24 end for

25 return matrix;

26

27

28Def  _SelectValidSamples(_outputMatrix, remainingMatrix,target,bins_)_:

/* This method selects all samples which are deemed valid from
"remainingMatrix" and stores them in "outputMatrix" */

/* First we get the index prices from the samples */

29 indexPrices = sumOverRows(remainingMatrix);

30 for  _bin in bins_ do

/* Then we select all samples which fall into a bin */

31 rowsInBins = remainingMatrix.where(indexPrices in bin);

32

33 tccNow we store at most N samples of rowsInBins N = targetVector(bin);

34 selectedRows = selectMaxRows(rowsInBins,N);

35 outputMatrix.add(selectedRows);

36 matrix.remove(selectedRows);

37 target(bin) $\leftarrow$ target(bin) - count(selectedRows)

38 end for

39 return outputMatrix;

40

41

42

Algorithm 1 Iterative Sort-Mix

##  Appendix B Pricing baskets with assets not contained in $N_{I}$

The proposed methodology focuses on pricing derivatives for a basket of assets
that are contained in the index. However, it can be extended to handle
derivatives on baskets of constituents not part of the index but presumed to
exhibit a similar correlation structure.

Consider a stock $S_{e}(t)$ where $e>N$. Suppose we aim to price derivatives
on a basket $N_{B}^{+}=N_{B}\cup{e}$, with $N_{B}\subset\mathbb{N}_{\leq N}$
representing a specific basket subset. Since there is no available information
about the correlation structure of $S_{e}$ with another asset in $N_{B}$,
coherent pricing of these derivatives is not possible. Yet, if we anticipate
comparable correlation characteristics between $S_{e}$ and the other basket
constituents (for instance, if the stocks share the same industry), we can
employ the calibrated copula from the index to also include the asset $e$.

To implement this, consider $s_{n}(t)$ as samples of a constituent (following
the order obtained after a certain step, e.g., step ii)), with $n\notin N_{B}$
to ensure these samples are not used in the basket calibration. Applying the
probability distribution function of $n$ to these samples allows us to obtain
uniformly distributed samples.

$u_{n}(t)=F_{S_{n}(t)}(s_{n}(t)).$ (B.1)

As the ordering of $u_{n}(t)$ remains consistent with other samples, the
copula remains unchanged. By projecting these samples using the probability
distribution of $S_{e}$, we derive samples represented as

$s_{e}(t)=F_{S_{e}(t)}(u_{n}(t)),$ (B.2)

which are now used to obtain the basket samples

$b(t)=\sum_{n\in N_{B}^{+}}s_{n}(t).$ (B.3)

This approach thus allows for the construction of baskets with up to $N$
external components.

##  Appendix C Proof of 3.3

The proof of 3.3 is the theoretical justification of a rearrangement
algorithm. Although conceptually not difficult, the proof is technically
involved. We will show the proof for a simplified model of two assets since
the generalization to multiple assets is trivial. Consider a model of two
assets $S_{1},S_{2}$. Since we only consider a fixed time, we drop the time
indication. The probability distribution of the sum of $S_{1}$ and $S_{2}$ is
given by $F_{S_{1}+S_{2}}(x)$.

For the rearrangement algorithm, we start with two independent random
variables, which are equal in distribution to $S_{1}$ and $S_{2}$. We denote
them as $\bar{S}_{1}$, $\bar{S}_{2}$ such that $F_{S_{1}}=F_{\bar{S}_{1}}$ and
$F_{S_{2}}=F_{\bar{S}_{2}}$. The proof of the theorem is simple: We first
assume that we can obtain samples of $S_{1}$ and $S_{2}$ directly. Since these
samples have all the desired properties, we will show that we can find an
ordering $\pi$, such that $\pi(\bar{s}_{1}),\pi(\bar{s}_{2})$ is very close to
the samples $s_{1},s_{2}$ obtained from $S_{1}$ and $S_{2}$

The first observation we make is that if $[s]_{M}$ is a sample of size $M$,
the empirical distributions $\hat{F}_{\bar{S}_{n},[s]_{M}}$ and
$\hat{F}_{S_{n},[s]_{M}}$ converge uniformly almost surely, as $M$ goes to
infinity. This is a consequence of the Glivenko-Cantelli theorem since both
empirical distributions converge uniformly a.s. to the same probability
distribution.

This means that we can fix $\epsilon$ and with probability 1, samples of
$[s]_{M}$ of size $M$ (of all 4 random variables) will be such that

$\left|\hat{F}_{S_{n},[s]_{M}}(x)-\hat{F}_{\bar{S}_{n},[s]_{M}}(x)\right|<\epsilon,\hskip
14.22636pt\forall x\in\mathbb{R},n\in\\{1,2\\}.$ (C.1)

Let us now define a permutation on $\bar{s}_{n}$, which will be used for this
purpose. We first define the ordering per”mutation
$o_{s_{n}}\colon\mathbb{N}_{M}\to\mathbb{N}_{M}$ as the permutation which
ordered a sample vector in ascending order:

$o_{s_{n}}(s_{n})=(s_{n}^{(1)},s_{n}^{(2)},\dots,s_{n}^{(M)}),$ (C.2)

where $s_{n}^{(1)}$ and $s_{n}^{(1)}$ are the smallest and largest sample of
$s_{n}$, respectively For both $n\in\\{1,2\\}$, we define $\pi_{n}$ as:

$\pi_{n}(\bar{s}_{n})=o^{-1}_{s_{n}}o_{\bar{s}_{n}}(\bar{s}_{n}).$ (C.3)

The proof now consists of showing that the following quantity can be bounded
based on $\epsilon$ for all $x$, so that it converges to $0$ as $M\to\infty$

$E(x):=\left|\hat{F}_{S_{1}+S_{2},[s]_{M}}(x)-\hat{F}_{\bar{S}_{1}+\bar{S}_{2},\pi([s]_{M})}(x)\right|.$
(C.4)

Based on the bound of the marginal distributions, we will now make the
following assumption, which leads to the worst possible bound of $E(x)$:

$\hat{F}_{S_{n},[s]_{M}}(x)+\epsilon=\hat{F}_{\bar{S}_{n},[s]_{M}}(x),\hskip
14.22636pt\forall x\in\mathbb{R},n\in\\{1,2\\}.$ (C.5)

Without loss of generality, we chose
$\hat{F}_{S_{n},[s]_{M}}(x)<\hat{F}_{\bar{S}_{n},[s]_{M}}(x)$. We find

$\hat{F}_{S_{1}+S_{2},[s]_{M}}(x)=\int_{0}^{\infty}\frac{1}{M}\left|\\{m\leq
M:s^{m}_{1}\leq y,s^{m}_{2}\leq x-y\\}\right|{\rm d}y$ (C.6)

and

$\hat{F}_{\bar{S}_{1}+\bar{S}_{2},\pi([s]_{M})}(x)=\int_{0}^{\infty}\frac{1}{M}\left|\\{m\leq
M:\bar{s}^{m}_{1}\leq y,\bar{s}^{m}_{2}\leq x-y\\}\right|{\rm d}y.$ (C.7)

The difference between the two quantities depends on the size of the set under
the integral:

$E(x)=\int_{0}^{\infty}\frac{1}{M}\left|\\{m\leq M:s^{m}_{1}\leq
y,s^{m}_{2}\leq x-y\\}\right|-\left|\\{m\leq M:\bar{s}^{m}_{1}\leq
y,\bar{s}^{m}_{2}\leq x-y\\}\right|{\rm d}y$ (C.8)

This term under the integral is of the following form $|A\cap B|-|C\cap D|$,
where $C\subset A$ and $D\subset B$ due to the assumption above. With the
general equation $|A\cap B|+|A\cup B|=|A|+|B|$, we then obtain:

$\displaystyle|A\cap B|-|C\cap D|$ $\displaystyle=|A|+|B|-|A\cup
B|-(|C|+|D|-|C\cup D|)$ (C.9) $\displaystyle=|A|-|D|+|B|-|C|-(|A\cup B|-|C\cup
D|)$ (C.10)

We observe now that since $(|A\cup B|-|C\cup D|)$ is certainly positive, it is
also bounded by

$(|A\cup B|-|C\cup D|)\leq|A|-|D|+|B|-|D|.$ (C.11)

For this reason, we can obtain the previous result that

$|A\cap B|-|C\cap D|\leq|A|-|D|+|B|-|D|<2\epsilon.$ (C.12)

Therefore, combining all the results, we have that

$E(x)\leq\int_{0}^{\infty}\frac{1}{M}\left||A\cap B|-|C\cap
D|\right|\leq\int_{0}^{T}2\epsilon.$ (C.13)

Since the samples $s_{1},s_{2},\bar{s}_{1},\bar{s}_{2}$ are almost surely
bounded, we obtain that

$E(x)\leq 2\epsilon B$ (C.14)

for some $B<\infty$