Efficient pricing and calibration of high-dimensional basket options

Abstract

This paper studies equity basket options – i.e. multi-dimensional derivatives whose payoffs depend on the value of a weighted sum of the underlying stocks – and develops a new and innovative approach to ensure consistency between options on individual stocks and the index comprising them. Specifically, we show how to resolve a well-known problem that when individual constituent distributions of an equity index are inferred from the single-stock option markets and combined in a multi-dimensional local/stochastic volatility model, the resulting basket option prices will not generate a skew matching that of the options on the equity index corresponding to the basket.

Keywords:

• Basket options
• index skew
• local volatility
• stochastic volatility
• collocation methods

Mathematic Subject Classifications:

• 91G20
• 91G30
• 91G50
• 60G07

1. Introduction

Basket options are contingent claims in which the underlying is a group of assets – typically equities (single stocks as well as equity indices/ETFs), commodities or currencies. As with standard options, the holder of a basket option has the right, but no obligation, to buy (call) or sell (put) the group at a specified strike price which itself is based on the weighted value of the component assets. Basket options can be traded on their own, but commonly feature implicitly in structured products like equity basket linked notes (ELNs) or equity basket certificates of deposit, where return to investors depends on the percentage change in the value of stocks or stock indices in the basket, with partial or total capital protection. For example, on 16 April 2021, JP Morgan priced $\mathrm{\$}$9.6 million of notes linked to a basket of 30 unequally weighted U.S. stocks with exposure to infrastructure.¹ The structure is effectively a combination of a European-style basket option and a hypothetical zero coupon bond of the issuing entity.

The key feature of basket options, and any structured notes based on them, is their inherently multidimensional risk profile which presents two sets of related challenges. First, there is the obvious computational problem of handling a potentially large number of correlated risk factors driving the basket, which in turn requires involved Monte Carlo simulation and calls for the deployment of some dimension-reduction technique, especially in light of the fact that many such structures exhibit early to exercise or other exotic features (e.g. Asian or lookback options).

Second, and perhaps more importantly, there is the problem of ensuring consistency between model dynamics chosen for the basket and its constituents, on the one hand, and market pricing on the other. Typically, stocks or stock indices/ETFs included in the basket will have their own vanilla options markets and the basket collection itself can also be traded on its own account. Thus, ensuring consistency entails not only fitting the implied volatility smiles/skews for individual assets in the basket – which can be done with excellent precision using e.g. the local volatility model or the Heston model – but also their covariance structure, such as to match closely the implied volatility smile for the market index corresponding to the basket.

Market practice is to use the local volatility model for individual basket members and historical correlations estimates between underliers to relate the respective risk factors in a Monte Carlo simulation.² In addition, the correlation matrix would normally need to be adjusted to ensure that it is positive semi-definite (e.g. by flooring of the eigenvalues at zero or even adopting a flat correlation level). However, this procedure does not automatically ensure that the model price of a basket option is consistent with the market price of the option on the index corresponding to the basket. And even if, as per market practice, the correlation coefficient is chosen specifically to match the price of at-the-money index options, the procedure generates an implied volatility smile/skew for the basket that is significantly less pronounced than that observed for the index options.

Consider for example options on the EURO STOXX Food and Beverage Index (Bloomberg ticker SX3E). The index currently comprises 11 large cap European equities from the food and beverage sector, such as French Danone, Dutch Heineken and Irish Kerry Group, each of which has its own listed options market. Matching 1Y ATM implied volatility for the index (15.51$\mathrm{\%}$) requires setting $\rho =0.22$ (data as of 1 April 2022). Repricing the basket with this flat correlation parameter for different moneyness levels generates an implied volatility skew of about 3.6 vol points,³ which is less than half of the skew observed in the market (7.5 vol points; cf. Figure 1, LHS). Similarly, fitting correlation to 1Y ATM vol does not allow to reproduce the term structure of index implied volatilities observed in the market (Figure 1, RHS).

Figure 1. Local volatility model calibrated to Euro STOXX Food and Beverages Index (SX3E) and a corresponding basket of its constituents; LHS: skew for 1Y index options; RHS: implied volatility term structure of index options; data as of 1 April 2022; correlation parameter $\rho =0.22$ calibrated to 1Y ATM implied SX3E volatility.

The above heuristic argument shows that a covariance structure implied by the modelling choice, like the multi-dimensional local volatility (or perhaps Heston) model, will not in general match the market well. This pattern becomes even more pronounced when handling really large baskets of tens or hundreds of dimensions, as the generated implied volatilities degenerates to almost a flat line, resembling implied volatilities produced by a simple Black-Scholes model. These shortcomings have been known in the literature for quite some time, as a number of papers have documented that when a constant correlation is picked to match the price of the ATM implied volatility of the index, it generates a skew that is smaller by roughly a factor of two than the market skew – as indeed shown in Figure 1 (see [3–5]).

One could argue at this point that since the implied volatilities for the index are generally available, there is no need to model individual underlying assets, and the problem of pricing a basket corresponding exactly to the index is a theoretical exercise with little practical relevance. Note, however, that by building a model for the basket that matches the covariance structure implied by index options (while simultaneously calibrating to individual stocks' options) we can easily price products whose value depends on any linear combination of any subset of stocks included in the broader index. Hence, in essence, we propose a tool for transferring market-implied covariance patterns between liquid and illiquid baskets. This allows us, for example, to price a basket of utilities stocks included in FTSE 100 or health care stocks in Dow Jones, while making sure that the prices of such illiquid baskets will be consistent with the covariance structures implied by prices of the respective liquid index options.

Several approaches have recently been proposed in the literature to tackle the problem of pricing basket options and their skews. Perhaps most notably, so called local correlation models – which come in various shapes and forms – allow correlation between stocks to depend on their prices, just like in standard local volatility models ([19,20,29]). Local correlation offered the first consistent remedy to the problem of insufficient skewness referred to above (cf. Figure 1; [3]), however not without a cost. Local correlation models tend to be slow, difficult to implement and lead to the usual problems of ensuring a positive semi-definite correlation matrix, all of which limits their practical usefulness in a truly multidimensional setting.

A related family of models proposes stochastic instantaneous correlation, usually introduced via Jacobi processes either in a pure form or with jumps (see e.g. [1,32]). Such models aim not only to fit the market prices for correlation products but also, more specifically, to address the so-called correlation skew, i.e. an empirically documented tendency for cross-correlations to rise when the market falls. Again, however, the added complexity leads to a further increase in dimensionality and the associated computational burden.

Still others propose to enforce the desired level of the skew through non-Gaussian copulas (e.g. [25,26]), however, the resulting model is ‘black boxy’ in nature and lacks dynamics, which means it can only be applied to products with a single maturity, as observed in [32]. Moreover, finding a multi-dimensional copula that isn't Gaussian can be challenging in practice, and anyway the need to calibrate the model through brute-force Monte Carlo makes the whole approach slow and costly.

Finally, [7] proposes an analytical framework to handle stochastic correlation via Wishart processes while [21] generalizes Merton's jump-diffusion option pricing approach to the multi-asset case, however without offering a comprehensive solution that would match both the basket skew as well as marginals of the constituent stocks.

In this study, we present a highly efficient and robust numerical model capable of handling large baskets with dimensions up to 500, or even more. Our model allows for calibration that is both quick and precise. By ‘quick', we mean that the calibration does not require search algorithms involving expensive computations, as it is based on the inversion of a marginal distribution. By ‘precise', we mean that the method guarantees perfect calibration to marginal distributions by construction, similar to Dupire's well-known non-parametric local volatility model [8]. Our local volatility model is constructed by mapping the basket onto a set of proxy variables, ensuring that both the marginal distributions for individual stocks and the covariance structure of the basket correspond to that of the traded index. To achieve excellent fitting of individual stocks' marginal distributions, we employ the Stochastic Collocation Monte Carlo Sampler (cf. [14]), a computationally inexpensive method for approximating complex distributions. Our model makes no assumptions about the models driving the respective stocks in the basket.

We then control the overall covariance structure of thus reconstituted basket via a judicious choice of ‘kernel’ processes, whereby – similarly as in [21] – we utilize Merton's jump-diffusion model [27] and show that the mean and volatility of Poisson jumps, together with cross-correlations of individual kernel processes, produce a rich enough dependence structure to fit both ATM index volatilities and skew with excellent precision.⁴ The proposed approach can be extended into a fully stochastic volatility setting; however, as we shall see, there appears to be no need for the added complexity as the classic Merton diffusion already delivers acceptable results. The framework presented in this article is in the spirit of the variance gamma-induced skew for the basket presented in [22,23]. However, our approach is of a local volatility type. It, therefore, does not require heavy calibration procedures for each constituent asset. For the sake of completeness, it is worth mentioning that model-independent approaches to basket option pricing are presented in [6,17].

The proposed framework offers highly efficient and almost perfect calibration to individual assets in the large baskets while controlling the basket's skew and thus facilitating efficient calibration to index options. Our approach does not ‘impose’ a model structure on the constituent assets, as in [23], but focuses on modelling the underlying distributions utilizing polynomial expansion, facilitating significant speed gain. As noted by Linders and Stassen in [23], calibrating each asset based on the Variance Gamma model requires substantial computational effort due to its reliance on Fourier-based techniques. Furthermore, the well-documented limitation of the Variance-Gamma model in terms of its difficulty in calibrating to the entire volatility surface, particularly for long-term expiries, remains a challenge. Moreover, the framework provides analytical pricing equations for the basket options, enabling fast calibration.

The rest of the paper proceeds as follows. Section 2 presents the general idea of the model and discusses the stochastic collocation approach to basket reconstitution. Section 3 discusses calibration strategy followed by a discussion of numerical exercises in Section 4. Finally, Section 5 draws conclusions.

2. The modelling framework

We consider a collection of N assets⁵ $S_1(t),S_2(t),\dots ,S_N(t)$ and define a basket, $B(\cdot )$, as consisting of $w_i$ shares of each individual stock $S_n(t)$, such that at any time t its price is given by: (1) $\begin{array}{r}B(t)=\sum _{i=1}^N{\omega }_i{S}_i(t),{\omega }_j\in \mathbb{R}.\end{array}$(1) In what follows, for modelling purposes, we shall distinguish between a basket $B(\cdot )$ as a portfolio of stocks, and a stock index, with price process $I(t)$, representing the underlying instrument for equity index options (e.g. an ETF), with realistic values of N ranging from 30 to 50 (DAX Index) up to hundreds (S&P 500 Index) or even thousands constituents (Russell 3000 Index).

Our approach to valuing structured payoffs on a basket $B(t)$ will be endowed with three critical features: (i) low dimensionality to facilitate quick and efficient pricing; (ii) consistency with the values of options on individual stocks making up the basket; and (iii) consistency with market prices of options on the stock index corresponding to the basket.

We pursue features (i)-(iii) above by building a one-dimensional local volatility process⁶: (2) $\begin{array}{r}\mathrm{d}\overline{B}(t)=r\overline{B}(t)\mathrm{d}t+{\sigma }_{LV}(t,\overline{B}(t))\overline{B}(t)\mathrm{d}W(t),\end{array}$(2) with r being an interest rate and ${\sigma }_{LV}$ the local volatility function, defined as follows: (3) $\begin{array}{r}{\sigma }_{LV}^2(t,k)=\frac{\frac{\mathrm{\partial }}{\mathrm{\partial }t}[{\mathrm{e}}^{-rt}{\int }_k^{\mathrm{\infty }}(y-k)f_{B(t)}(y)\mathrm{d}y]+rk(F_{B(t)}(k)-1)}{\frac{1}{2}{k}^2{f}_{B(t)}(k)},\end{array}$(3) where $F_{B(t)}(\cdot )$ is the CDF of the basket $B(t)$ and $f_{B(t)}(\cdot )>0$ the corresponding PDF. We shall require that the approximate model for the basket $\overline{B}(t)$ generates an implied volatility smile/skew matching the implied volatility surface observed in the market for options on the corresponding stock index. This will ensure that prices of index options derived from our basket model will by definition match perfectly (cf. [9,10,15]). Furthermore, it will also guarantee consistency between model-derived prices of baskets on any subset $\{S_{i_1},\dots ,S_{i_K}\}\subseteq \{S_1,\dots ,S_N\}$ of the N stocks included in the index and the prices of liquid index options.

In practical terms, our approach entails two critical steps. First, for each $T_j$ we reconstitute the basket $B(T_j)$ by projecting each process $S_i$ – assumed to be known only through its marginal distributions/options prices – on a polynomial of synthetic variables $g_{i,j}(X_i(T_j))$, such that (4) $\begin{array}{r}{S}_i(T_j)\approx {\mathrm{e}}^{g_{i,j}(X_i(T_j))},\mathrm{and}B(T_j)\approx \sum _{i=1}^N{\omega }_i{\mathrm{e}}^{g_{i,j}(X_i(T_j))}.\end{array}$(4) The construction of $g_{i,j}(\cdot )$ will be described in detail below, but its main role will be to ensure the fit of marginal distributions for any choice of $X_i(T_j)$. The reconstitution of the basket paves the way for the second stage of the modelling procedure, namely imposing a multivariate dependence structure on the surrogate kernel processes such as to match the index implied volatility skew, or equivalently, matching marginal distributions between the basket and the index for each maturity point $T_j$, i.e. ensuring that: $B(T_j)\stackrel{d}{=}I(T_j).$ This step will be handled through parameters of the ‘kernel’ processes $X_j(t)$, which we assume to be driven by Merton-type jump diffusions. This procedure will allow us to both calibrate the model to individual stocks, and control/modify the covariance structure of all the underlying assets.

2.1. Covariance structure of the basket and the ‘Leaking correlation’ problem

At this stage it may be worthwhile to pause for a moment and provide some further motivation for why controlling basket dependence structure is so tricky and a judicious choice of kernel processes is required. After all, it could be argued that the model price of a basket option can be fitted to a desired market level by altering the correlations ${\rho }_{i,j}$ between the Brownian motions driving individual stocks $S_i,S_j$. Such correlation coefficients would then essentially become model inputs describing the stochastic nature of the underlying assets. However, as discussed in [2], correlations can be expected to ‘leak’, i.e. dissolve, in a Monte Carlo setting, so that effective correlations between the respective basket constituents ends up weaker than assumed. Clearly, this will limit the extent to which we may be able to control the basket smile generated by our model, leading to erroneous pricing of derivatives away from at-the-money levels.

To illustrate the correlation leakage problem, consider two stock exchange processes $S_1(t)$ and $S_2(t)$ driven by the Heston model with correlated Brownian motions.⁷ In the simulation of the underlying processes, we consider a simple Euler discretization as in [28]. In Figure 2 the cases for the Feller's condition⁸ are presented for a range of correlation values. We see that the imposed correlations are not preserved in time, even when Feller's condition is satisfied. As a potential remedy to ‘leaking correlation’, we could attempt to enforce the desired dependence patterns via copula. This, however, does not guarantee success either. To see why, consider the previous example of the Heston model but now – to make the exposition clearer – assume that the two stocks are independent, i.e. $\mathrm{d}{W}_{1,1}(t)\mathrm{d}{W}_{2,1}(t)=0\cdot \mathrm{d}t$. In such a setup the marginal, cumulative distribution function of each stock can be recovered using Fourier transform (cf. [28]), and the covariance structure can then be imposed using a copula. The basket, $B(T)$, at given maturity T is thus computed via: $\begin{array}{r}B(T)=F_{S_1(T)}^{-1}(U_1)+F_{S_2(T)}^{-1}(U_2)=S_1(T)+S_2(T),\end{array}$where the uniformly distributed random variables $U_1$ and $U_2$ are joined using a Gaussian copula ${\mathcal{C}}_{\mathcal{R}}.$ Using such a model we can now price a call option and compute the corresponding implied volatility for a range of copula parameters.

Figure 2. Correlation, $(\log S_1(t),\log S_2(t))$, as a function of time. Left: Feller condition satisfied; Right: Feller condition not satisfied.

Figure 3 shows that varying the copula correlation has a significant impact on the level of corresponding implied volatilities. However, it is also clear that while varying the covariance structure alone might help in matching, say ATM index volatilities, Gaussian copula is unlikely to generate a desired level of skewness – a problem we already alluded to in the introduction. One potential remedy would be to consider more elaborate copula functions (as e.g. proposed in [25]), however we opt for a different approach. As has already been hinted above, and as we shall explain in detail below, choosing Merton's jump diffusion dynamics [27] for the synthetic variables mapping individual stocks offers great flexibility in fitting both ATM index volatilities and skew with excellent precision.

Figure 3. Implied volatilities for a basket $B(T)$ for d = 2 driven by the Heston model.

2.2. Basket reconstruction via projection with collocation

As discussed before, our ELV model relies on the availability of the marginal distribution of each asset, $S_i(t)$, at every time $T_j$. Such a distribution may be given explicitly via a specific form of the SDE governing the evolution of $S_i$. Alternatively, the relevant distributions can be inferred from options prices on $S_i$ – either obtained from the market or generated by some unknown model – leaving us completely agnostic as to the shape or form of the asset-pricing model. In this section we follow the latter more general approach and derive the distribution through the following well known relation: (5) $\begin{array}{r}{F}_{S_i(T_j)}(y)=1+\frac{\mathrm{\partial }{V}_c(t_0,T_j,S_i(t_0),K))}{\mathrm{\partial }K}{|}_{K=y},f_{S_i(T_j)}(y)=\frac{{\mathrm{\partial }}^2{V}_c(t_0,T_j,S_i(t_0),K))}{\mathrm{\partial }{K}^2}{|}_{K=y},\end{array}$(5) where $F_{S_i(T_j)}(y)$ is the CDF of $S_i(\cdot )$ observed at time $T_j$ and $f_{S_i(T_j)}(y)$ is the corresponding PDF. We shall take such distributions as given and ultimately try to reconstruct the marginal distributions of the basket $B(T_j)$ defined as the weighted sum of all individual assets. This step will be ultimately performed by ‘coupling’ the individual marginals via correlated kernel processes.

However, our first goal is to project each basket constituent, $S_i(T_j)$, on a polynomial, $g_{i,j}^m(\cdot )$, of synthetic variables $X_i(T_j)$, such as to ensure a perfect fit to the respective marginal distributions for every set of model parameters $X_i(T_j)$. To this end we resort to the so-called Stochastic Collocation Monte Carlo Sampler (SCMC) discussed in [14,24]. Using the SCMC method we can express $\log S_i(T_j)$ as: (6) $\begin{array}{r}\log S_i(T_j)\approx g_{i,j}^m(X_i(T_j))=\sum _{k=0}^{m-1}{s}_{i,j,k}{\ell }_k(x_{i,j,k},X_i(T_j)),\end{array}$(6) where ${\ell }_k(\cdot )$ are the Lagrange basis function evaluated at the collocation points $x_{i,j,k}$, $k=0,\dots ,m-1$ based on the kernel variable $X_i(T_j)$. Alternatively, the representation above can be re-expressed using a polynomial representation, (7) $\begin{array}{r}{g}_{i,j}^m(X_i(T_j))=\sum _{k=0}^{m-1}{\hat{\alpha }}_{i,j,k}{X}_i^k(T_j)={\hat{\alpha }}_{i,j,0}+{\hat{\alpha }}_{i,j,1}{X}_i(T_j)+\cdots +{\hat{\alpha }}_{i,j,m-1}{X}_i^{m-1}(T_j),\end{array}$(7) where coefficients ${\hat{\alpha }}_{i,j,k}$ are known explicitly [13]. The expansion in (6(6) $\begin{array}{r}\log S_i(T_j)\approx g_{i,j}^m(X_i(T_j))=\sum _{k=0}^{m-1}{s}_{i,j,k}{\ell }_k(x_{i,j,k},X_i(T_j)),\end{array}$(6) ) ensures that the CDFs of $S_i(T_j)$ and $g_{i,j}^m(X_i(T_j))$ agree on the so-called collocation points, i.e. $\begin{array}{r}{F}_{\log S_i(T_j)}(s_{i,j,k})=F_{g_m(X_i(T_j))}(s_{i,j,k}),k=0,\dots ,m-1,\end{array}$where $s_{i,j,k}=F_{\log S_i(T_j)}^{-1}(F_X(x_{i,j,k}))$ and where $X_i$ is the collocating variable.

The choice regarding variables $X_i(T_j)$, although in principle free, is crucial for fitting to the marginal distribution of $S_i(T_j)$ and controlling the basket's covariance structure. The most appealing choice regarding variables $X_i$ is a normal distribution [14], which would imply sampling from multivariate Gaussian copula, thus a computationally cheap, numerical procedure. However, although Gaussian random variables are sufficient for finding optimal ${\hat{\alpha }}_{i,j,k}$ in (6(6) $\begin{array}{r}\log S_i(T_j)\approx g_{i,j}^m(X_i(T_j))=\sum _{k=0}^{m-1}{s}_{i,j,k}{\ell }_k(x_{i,j,k},X_i(T_j)),\end{array}$(6) ) guaranteeing excellent fit to marginal distribution of $S_i(T_j)$ such a model will generate low skew of the basket and therefore the final calibration to the index will not be satisfactory (see the remarks in Section 2.1). As presented later in the article, the non-trivial covariance structure of the underlying assets is necessary for controlling the implied volatility skew. This is confirmed by empirical studies [3].

In the ELV model presented in this article, the kernel that is used is Merton's jump diffusion model [27], which under the risk-neutral measure consists of a Brownian motion and a compound Poisson process, which is defined by: (8) $\begin{array}{r}{X}_i(t)=X_i(t_0)+(r-\frac{1}{2}{\sigma }_i^2-{\xi }_{p,i}\mathbb{E}[{\mathrm{e}}^{J_i}-1])t+{\sigma }_i{W}_i(t)+\sum _{k=1}^{X_{\mathcal{P},i}(t)}{J}_{k,i},\end{array}$(8) with $X_i(t_0)=\log S_i(t_0)$, for ${\sigma }_i>0$, Brownian motion, $W_i(t)$, Poisson process $X_{\mathcal{P},i}(t)$, $t\ge 0$ with parameter ${\xi }_{p,i}$ such $\mathbb{E}[X_{\mathcal{P},i}(t)]={\xi }_{p,i}t$ and where the jump sizes $J_{k,i}$ are i.i.d. given by $J_{k,i}\sim \mathcal{N}({\mu }_{J,i},{\sigma }_{J,i}^2)$. For a given index i the jumps, $J_{k,i}$, Brownian motion, $W_i(t)$, and Poisson process, $X_{\mathcal{P},i}(t)$ are assumed to be independent.

Remark

Alternative dynamics for $X_i$

Alternatively, one may also consider the more advanced process structure $X_i(\cdot )$ by adding stochastic volatility. Here, however, we will focus on a simplified model, and we will show that such a structure is sufficient for excellent calibration of the market data, even for large baskets consisting of 30 assets.

Once the collocating variable has been chosen, we need to establish the mapping procedure from the asset $\log S_i(T_j)$ onto the surrogate variable $X_i(\cdot ).$ In order to achieve accurate approximations one needs to firstly determine the so-called collocation points, $x_{i,j}$, $j=1,\dots ,N$, of ‘surrogate’ variable $X_t(\cdot )$- these collocation points are based on moments of $X_i(t)$. This ensures that the collocation points are the zeros of the orthogonal polynomial corresponding to the distribution of $X_i(\cdot )$, and we can establish the connection with the computation of integrals by Gauss quadrature. Typically, for N-collocation points $2N$ moments of $X_i(\cdot )$ are needed. However, for $X_i(\cdot )$ defined in (8(8) $\begin{array}{r}{X}_i(t)=X_i(t_0)+(r-\frac{1}{2}{\sigma }_i^2-{\xi }_{p,i}\mathbb{E}[{\mathrm{e}}^{J_i}-1])t+{\sigma }_i{W}_i(t)+\sum _{k=1}^{X_{\mathcal{P},i}(t)}{J}_{k,i},\end{array}$(8) ) one may, due to the presence of jumps, expect an effect often called ‘moment explosion’ which translates to huge moments causing significant numerical instability of the method. As an alternative approach, we use the collocation points $z_k$, $k=0,\dots ,m-1$ corresponding to the standard normal random variable⁹ $Z\stackrel{\mathrm{d}}{=}\mathcal{N}(0,1)$. Such a strategy implies, for a given time $T_j$, the following relation between $\log S_i(T_j)$, $X_i(\cdot )$ and Z: (9) $\begin{array}{rl}{x}_{i,j,k}& =F_{X_i(T_j)}^{-1}(F_Z(z_k)),\end{array}$(9) (10) $\begin{array}{rl}\log s_{i,j,k}& =F_{\log S_i(T_j)}^{-1}(F_{X_i(T_j)}(x_{i,j,k})),\end{array}$(10) which illustrates how to avoid computation of the collocation points of $X_i(T_j)$ using moments of $X_i(T_j).$ It is important to note that although variable $X_i(T_j)$ is not directly used to compute the collocation points of $S_i(T_j)$ the approximating function $g_{i,j}^m(X_i(T_j))$ is explicitly defined in terms of $X_i$ and points $x_{i,j,k}$ and $\log s_{i,j,k}$. The details on the optimality of such ‘substitution’ can be found in [31].

Remark

Grid stretching

When dealing with heavy-tailed distributions (of leptokurtic type or the distribution is highly skewed), as here proposed for $X_i(T_j)$, to avoid numerically unstable inversions $F_{X_i(T_j)}^{-1}(a)$ for either $a\to 0$ or $a\to 1$ it is recommended to use the so-called ‘grid stretching’ technique [14] which facilitates stable inversions, especially for a high number of the collocation points, m.

Availability of the ChF in closed form allows finding, via Fourier inversion [28], the corresponding CDF, $F_{X_i(T_j)}(x)$, and thus finding of the corresponding collocation points, $x_{i,j,k}$. (11) $\begin{array}{r}{\varphi }_{X_i(t)}(u)=\mathbb{E}\left[{\mathrm{e}}^{iu{X}_i(t)}\right]=\exp (iu(X_i(t_0)+{\overline{\mu }}_{i}t)-\frac{1}{2}{\sigma }_i^2{u}^{2}t+{\xi }_{p,i}t{\mathrm{e}}^{iu{\mu }_{J,i}-\frac{1}{2}{u}^2{\sigma }_{J,i}^2}),\end{array}$(11) with ${\overline{\mu }}_i=r-\frac{1}{2}{\sigma }_i^2-{\xi }_{p,i}({\mathrm{e}}^{{\mu }_{J,i}+\frac{1}{2}{\sigma }_{J,i}^2}-1).$

Then, by utilizing the COS method [28], one is able to determine the corresponding CDF: (12) $\begin{array}{r}{F}_{X_i(T_j)}(y)\approx \sum _{k=0}^{N_c-1}{\overline{F}}_k\psi (a,b,y),{\overline{F}}_k=\frac{2}{b-a}\mathrm{\Re }\{{\varphi }_{X_i(T_j)}\left(\frac{k\pi }{b-a}\right)\cdot \exp (-i\frac{ka\pi }{b-a})\},\end{array}$(12) with (13) $\begin{array}{r}\psi (a,b,y)=\{\begin{array}{ll}\frac{b-a}{k\pi }\sin \left(\frac{k\pi (y-a)}{b-a}\right),& \mathrm{for}k=1,2,\dots ,N_c-1,\\ (x-a),& \mathrm{for}k=0,\end{array}\end{array}$(13) where $N_c$ indicates the number of the expansion terms, a and b are the domain parameters, typically determined in terms of cummulants.

In the next step, we build a multivariate dependence structure between all the marginal distributions using a system of correlated SDEs. The basket, $B(T_j)$, yields, in terms of the latent variables, the following form: (14) $\begin{array}{rl}B(T_j)& =\sum _{i=1}^N{S}_i(T_j)\approx \sum _{i=1}^N\exp (g_{i,j}^m(X_i(T_j)))\\ & =\sum _{i=1}^N\exp (\sum _{k=0}^{m-1}{\hat{\alpha }}_{i,j,k}{X}_i^k(T_j)),\end{array}$(14) with $X_i(T_j)$ defined in (8(8) $\begin{array}{r}{X}_i(t)=X_i(t_0)+(r-\frac{1}{2}{\sigma }_i^2-{\xi }_{p,i}\mathbb{E}[{\mathrm{e}}^{J_i}-1])t+{\sigma }_i{W}_i(t)+\sum _{k=1}^{X_{\mathcal{P},i}(t)}{J}_{k,i},\end{array}$(8) ). Although we have initially assumed that each $X_i$ has its own independent jumps $J_{k,i}$ controlled by $X_{{\mathcal{P}}_i(t)}$, our numerical studies have shown that without sacrificing much in terms of skew fitting possibilities, one can choose a sparse form of the model with $X_{{\mathcal{P}}_i(t)}$ being the same for all the assets, i.e. $\begin{array}{r}{X}_{{\mathcal{P}}_i(t)}\equiv X_{\mathcal{P}(t)}.\end{array}$In such form, the dependence between the individual stocks in the basket is controlled via two main elements: correlated Brownian motions $W_i(T_j)$ and common Poisson process $X_{\mathcal{P}}(T_j)$ for all Merton's processes, $X_i.$ Thus, the basket is reconstructed using a combination of marginal distributions for each asset. The basket PDF, $f_{B(t)}$, can then be obtained by Monte Carlo. It should be stressed that the sampling from multivariate normal is computationally inexpensive. The expensive part involves the inversion of marginal distributions, but as explained above, we handle this using collocation, mitigating the computational cost quite considerably.

Importantly, the proposed methodology does not rely on a calibration of the process for marginal distributions, but solely on the computation of ${\hat{\alpha }}_{i,j,k}$, for $i=1,\dots ,N$, $j=1,\dots ,N_T,$ $k=1,\dots ,m$, and inexpensive samples from the multidimensional distribution of $(X_i(\cdot ),\dots ,X_N(\cdot ))$.

Remark

Marginal distribution

We would like to stress that due to the collocation method and the CDF mapping procedure, determined coefficients ${\hat{\alpha }}_{i,j,k}$ will ensure agreement between the target CDF, $F_{\log S_i(T_j)}(\cdot )$ and the surrogate, $F_{X_i(T_j)}(\cdot ).$ Once the marginal distributions are calibrated for any configurations of model parameters of $X_i(\cdot )$, the smile/skew of the index needs to be matched. Later, in this article, we will show that the model proposed allows for intuitive control of different volatility shapes of the basket, therefore facilitating calibration to the index.

3. Basket calibration to the index

Given that the ELV model, by its construction, guarantees fit to marginal distributions one still needs to calibrate the basket to index volatilities. The brute force strategy of such optimization can be performed using Monte Carlo simulation, i.e. for different model parameters, the samples of $X_i(T_j)$ are fed to basket equation in (14(14) $\begin{array}{rl}B(T_j)& =\sum _{i=1}^N{S}_i(T_j)\approx \sum _{i=1}^N\exp (g_{i,j}^m(X_i(T_j)))\\ & =\sum _{i=1}^N\exp (\sum _{k=0}^{m-1}{\hat{\alpha }}_{i,j,k}{X}_i^k(T_j)),\end{array}$(14) ) and the option prices are computed. Such routine is typically sub-optimal as it would require multiple iterations and samplings for $X_i(\cdot ).$ To speed up the process and save computational burden, we opt for a different approach, based on moment matching.

Recall that in its most general form the model might appear to have an overwhelming number of parameters to be calibrated. Fortunately, as we have already hinted above, many simplifications can be resorted to without sacrificing calibration precision, while improving tractability and sparsity.

Firstly, we note that the ELV model's heart lies in the idea of construction of the local volatility function; therefore, only the marginal distribution for the basket $B(T_j)$ is relevant. This implies that the calibration can be performed sequentially, or parallel, for every time $T_j$, $j=1,\dots ,N_T.$ Secondly, we assume that the individual stocks share not only the type of the kernel process but also the driving parameters. Specifically, we assume henceforth that all the assets have a common counting process $X_{\mathcal{P}}(T_j)$, independent i.i.d. jumps $J\equiv J_{k,i}\sim \mathcal{N}({\mu }_J,{\sigma }_J^2)$ and correlated Brownian motions $\mathrm{d}{W}_i(t)\mathrm{d}{W}_j(t)=\rho \mathrm{d}t,$ which leads to the following simplified representation: $\begin{array}{r}{X}_i(T_j)=X_i(t_0)+(r-\frac{1}{2}{\sigma }^2-{\xi }_p\mathbb{E}[{\mathrm{e}}^J-1])T_j+\sigma W_i(T_j)+\sum _{k=1}^{X_{\mathcal{P}}(T_j)}{J}_{k,i}.\end{array}$Thus, in order to calibrate the basket $B(T_j)$ to an index we need, for every $T_j$, to determine the set of optimal model parameters σ, ${\xi }_p$, ${\mu }_J$ and ${\sigma }_J$, driving the model's covariance structure.

Each of the jump parameters has a different effect on the shape of the implied volatility curve, i.e. ${\sigma }_J$ has a significant impact on the curvature, ${\xi }_p$ controls the overall level of the implied volatility, whereas ${\mu }_J$ influences the implied volatility slope (the skew).

Having already gained much in terms of tractability, we can further save the computational burden by performing calibration via the moment matching procedure. Given that the density of the index, $f_I(y)$, exists and can be calculated from options prices, as per Equation (5(5) $\begin{array}{r}{F}_{S_i(T_j)}(y)=1+\frac{\mathrm{\partial }{V}_c(t_0,T_j,S_i(t_0),K))}{\mathrm{\partial }K}{|}_{K=y},f_{S_i(T_j)}(y)=\frac{{\mathrm{\partial }}^2{V}_c(t_0,T_j,S_i(t_0),K))}{\mathrm{\partial }{K}^2}{|}_{K=y},\end{array}$(5) ), the objective is to match moments from the index, $I(t)$, and the basket $B(t)$. It is important to note that the polynomial coefficients ${\alpha }_{i,j,k}$ are implicitly a function of model parameters, i.e. (15) $\begin{array}{r}{\alpha }_{i,j,k}:={\alpha }_{i,j,k}(\sigma ,{\xi }_p,{\mu }_J,{\sigma }_J),\end{array}$(15) therefore any change in the model parameters requires re-computation of these coefficients. Since this procedure only relies on the computation of a few inversions of the CDF in Equation (12(12) $\begin{array}{r}{F}_{X_i(T_j)}(y)\approx \sum _{k=0}^{N_c-1}{\overline{F}}_k\psi (a,b,y),{\overline{F}}_k=\frac{2}{b-a}\mathrm{\Re }\{{\varphi }_{X_i(T_j)}\left(\frac{k\pi }{b-a}\right)\cdot \exp (-i\frac{ka\pi }{b-a})\},\end{array}$(12) ), it is computationally cheap.

Since the first moment is trivial let us start with the second moment. By definition we have: $\begin{array}{r}\mathbb{E}[B^2(T_j)]=\sum _{i_1=1}^N\sum _{i_2=1}^N{\omega }_{i_1}{\omega }_{i_2}\mathbb{E}[S_{i_1}(T_j)S_{i_2}(T_j)].\end{array}$One the other hand, assuming that we are able to determine the correlation $\hat{\rho }(S_{i_1},S_{i_2})$ between $S_{i_1}(T_j)$ and $S_{i_2}(T_j)$ (this correlation will differ from ${\rho }_{i_1,i_2}\mathrm{d}t=\mathrm{d}{W}_{i_1}(t)\mathrm{d}{W}_{i_2}(t)$ between the corresponding Brownian motions), we find: $\begin{array}{r}\hat{\rho }(S_{i_1},S_{i_2})\stackrel{\mathrm{def}}{=}\frac{\mathbb{E}[S_{i_1}(T_j)S_{i_2}(T_j)]-\mathbb{E}[S_{i_1}(T_j)]\mathbb{E}[S_{i_2}(T_j)]}{{\sigma }_{i_1}{\sigma }_{i_2}},\end{array}$where ${\sigma }_{i,1}$ and ${\sigma }_{i,2}$ are standard deviations of $S_{i_1}(T_j)$ and $S_{i_2}(T_j)$ respectively (these quantities can be calculated directly using (5(5) $\begin{array}{r}{F}_{S_i(T_j)}(y)=1+\frac{\mathrm{\partial }{V}_c(t_0,T_j,S_i(t_0),K))}{\mathrm{\partial }K}{|}_{K=y},f_{S_i(T_j)}(y)=\frac{{\mathrm{\partial }}^2{V}_c(t_0,T_j,S_i(t_0),K))}{\mathrm{\partial }{K}^2}{|}_{K=y},\end{array}$(5) )). Given that the marginal distributions of $S_{i_1}$ and $S_{i_2}$ are available, both: first moments and the standard deviations can be easily obtained, by for example proper integration of the market-implied densities in (5(5) $\begin{array}{r}{F}_{S_i(T_j)}(y)=1+\frac{\mathrm{\partial }{V}_c(t_0,T_j,S_i(t_0),K))}{\mathrm{\partial }K}{|}_{K=y},f_{S_i(T_j)}(y)=\frac{{\mathrm{\partial }}^2{V}_c(t_0,T_j,S_i(t_0),K))}{\mathrm{\partial }{K}^2}{|}_{K=y},\end{array}$(5) ). Then, we have: (16) $\begin{array}{r}\mathbb{E}[S_{i_1}(T_j)S_{i_2}(T_j)]=\hat{\rho }(S_{i_1},S_{i_2}){\sigma }_{i_1}{\sigma }_{i_2}+\mathbb{E}[S_{i_1}(T_j)]\mathbb{E}[S_{i_2}(T_j)],\end{array}$(16) therefore the second moment of the basket becomes: (17) $\begin{array}{r}\mathbb{E}[B^2(T_j)]=\sum _{i_1=1}^N\sum _{i_2=1}^N{\omega }_{i_1}{\omega }_{i_2}(\hat{\rho }(S_{i_1},S_{i_2}){\sigma }_{i_1}{\sigma }_{i_2}+\mathbb{E}[S_{i_1}(T_j)]\mathbb{E}[S_{i_2}(T_j)]).\end{array}$(17) Since the correlation, $\hat{\rho }(S_{i_1},S_{i_2})$, measures the linear relationship between assets $S_{i_1}$ and $S_{i_2}$, we will compute this coefficient based on a fewer number of the collocation points ,¹⁰ m = 2 which yields the following representation: (18) $\begin{array}{rl}\hat{\rho }(S_{i_1},S_{i_2})& \approx \rho ({\mathrm{e}}^{{\hat{\alpha }}_{i_1,j,0}+{\hat{\alpha }}_{i_1,j,1}{X}_{i_1}(T_j)},{\mathrm{e}}^{{\hat{\alpha }}_{i_2,j,0}+{\hat{\alpha }}_{i_1,j,1}{X}_{i_2}(T_j)})\\ & =\rho ({\mathrm{e}}^{{\hat{\alpha }}_{i_1,j,1}{X}_{i_1}(T_j)},{\mathrm{e}}^{{\hat{\alpha }}_{i_2,j,1}{X}_{i_2}(T_j)})=:\overline{\rho }(S_{i_1},S_{i_2}),\end{array}$(18) with $X_i(\cdot )$ defined in (8(8) $\begin{array}{r}{X}_i(t)=X_i(t_0)+(r-\frac{1}{2}{\sigma }_i^2-{\xi }_{p,i}\mathbb{E}[{\mathrm{e}}^{J_i}-1])t+{\sigma }_i{W}_i(t)+\sum _{k=1}^{X_{\mathcal{P},i}(t)}{J}_{k,i},\end{array}$(8) ), with correlated Brownian motions, $\begin{array}{r}\overline{\rho }(S_{i_1},S_{i_2})\stackrel{\mathrm{def}}{=}\frac{\mathbb{E}[{\mathrm{e}}^{{\hat{\alpha }}_{i_1,j,1}{X}_{i_1}(T_j)+{\hat{\alpha }}_{i_2,j,1}{X}_{i_2}(T_j)}]-\mathbb{E}[{\mathrm{e}}^{{\hat{\alpha }}_{i_1,j,1}{X}_{i_1}(T_j)}]\mathbb{E}[{\mathrm{e}}^{{\hat{\alpha }}_{i_2,j,1}{X}_{i_2}(T_j)}]}{\sqrt{\mathbb{V}\mathrm{ar}[{\mathrm{e}}^{{\hat{\alpha }}_{i_1,j,1}{X}_{i_1}(T_j)}]\mathbb{V}\mathrm{ar}[{\mathrm{e}}^{{\hat{\alpha }}_{i_2,j,1}{X}_{i_2}(T_j)}]}}.\end{array}$We immediately note the connection between the moments and the corresponding ChF of $X_i(t)$ defined in (11(11) $\begin{array}{r}{\varphi }_{X_i(t)}(u)=\mathbb{E}\left[{\mathrm{e}}^{iu{X}_i(t)}\right]=\exp (iu(X_i(t_0)+{\overline{\mu }}_{i}t)-\frac{1}{2}{\sigma }_i^2{u}^{2}t+{\xi }_{p,i}t{\mathrm{e}}^{iu{\mu }_{J,i}-\frac{1}{2}{u}^2{\sigma }_{J,i}^2}),\end{array}$(11) ) resulting in, $\begin{array}{r}{\varphi }_{X_i(t)}(u)=\mathbb{E}[{\mathrm{e}}^{iu{X}_i(t)}],\mathrm{thus}\mathbb{E}[{\mathrm{e}}^{a{X}_i(t)}]={\varphi }_{X_i(t)}(-ia),a\in \mathbb{R},\end{array}$resulting in the following representation (19) $\begin{array}{r}\overline{\rho }(S_{i_1},S_{i_2})=\frac{\mathbb{E}[{\mathrm{e}}^{c_1{X}_{i_1}(T_j)+c_2{X}_{i_2}(T_j)}]-{\varphi }_{X_{i_1}}(-i{c}_1){\varphi }_{X_{i_2}}(-i{c}_2)}{\sqrt{({\varphi }_{X_{i_1}}(-2i{c}_1)-{\varphi }_{X_{i_1}}^2(-i{c}_1))({\varphi }_{X_{i_2}}(-2i{c}_2)-{\varphi }_{X_{i_2}}^2(-i{c}_2))}},\end{array}$(19) with $c_1={\hat{\alpha }}_{i_1,j,1}$, $c_2={\hat{\alpha }}_{i_2,j,1}$ and ${\varphi }_{X_i(T_j)}$ defined in Equation  (11(11) $\begin{array}{r}{\varphi }_{X_i(t)}(u)=\mathbb{E}\left[{\mathrm{e}}^{iu{X}_i(t)}\right]=\exp (iu(X_i(t_0)+{\overline{\mu }}_{i}t)-\frac{1}{2}{\sigma }_i^2{u}^{2}t+{\xi }_{p,i}t{\mathrm{e}}^{iu{\mu }_{J,i}-\frac{1}{2}{u}^2{\sigma }_{J,i}^2}),\end{array}$(11) ). The only term that caries the correlation between both processes is the expectation involving correlated Merton's processes. This expectation is given in closed form and it is presented in Lemma 3.1 below.

Lemma 3.1

Expectation for a 2D Merton's model

For a given 2D Merton's model defined as: (20) $\begin{array}{r}{X}_1(t)={\mu }_1+{\sigma }_1{W}_1(t)+\sum _{k=1}^{X_{\mathcal{P}}(t)}{J}_{k,1},X_2(t)={\mu }_2+{\sigma }_2{W}_2(t)+\sum _{k=1}^{X_{\mathcal{P}}(t)}{J}_{k,2},\end{array}$(20) with correlated Brownian motions $\mathrm{d}{W}_1(t)\mathrm{d}{W}_2(t)={\rho }_{1,2}\mathrm{d}t,$ common counting Poisson process $X_{\mathcal{P}}(t)$ with the intensity parameter ${\xi }_p$ and independent identically distributed jumps, $J\equiv J_{k,\cdot }\sim \mathcal{N}({\mu }_J,{\sigma }_J^2)$, the following expression holds: (21) $\begin{array}{rl}\mathbb{E}[{\mathrm{e}}^{a{X}_1(T)+b{X}_2(T)}|\mathcal{F}(t_0)]& =\exp [a{\mu }_1+b{\mu }_2+\frac{1}{2}{\hat{\sigma }}^{2}T+{\xi }_{p}T({\mathrm{e}}^{(a+b){\mu }_J+\frac{1}{2}(a^2+b^2){\sigma }_J^2}-1)]\\ & =:{\omega }_X(a,b),\end{array}$(21) where ${\hat{\sigma }}^2=a^2{\sigma }_1^2+b^2{\sigma }_2^2+2{\rho }_{1,2}ab{\sigma }_1{\sigma }_2$.

Proof.

Proof can be found in Appendix 1.

Utilizing the results above we find the closed form solution for the expectation and therefore the correlation coefficient, $\overline{\rho }(S_{i_1},S_{i_2})$, by simply setting where $c_1={\hat{\alpha }}_{i_1,j,1}$, $c_2={\hat{\alpha }}_{i_2,j,1}$, ${\mu }_i=X_i(t_0)+(r-\frac{1}{2}{\sigma }^2-{\xi }_p\mathbb{E}[{\mathrm{e}}^J-1])T$, and ${\sigma }_1={\sigma }_2=\sigma .$

Once the correlation coefficient, $\overline{\rho }(S_{i_1},S_{i_2})$, is determined one is able to perform the ATM calibration of the basket to the index. (22) $\begin{array}{r}\underset{{\sigma }_J,{\mu }_J,{\xi }_p,\rho }{min}{(\mathbb{E}[B^2(T_j)]-\mathbb{E}[I^2(T_j)])}^2,\mathbb{E}[I^2(T_j)]={\int }_{\mathbb{R}}{x}^2{f}_{I(T_j)}(x)\mathrm{d}x,\end{array}$(22) where the index density $f_{I(T_j)}(x)$ is implied from the option quotes for the index (5(5) $\begin{array}{r}{F}_{S_i(T_j)}(y)=1+\frac{\mathrm{\partial }{V}_c(t_0,T_j,S_i(t_0),K))}{\mathrm{\partial }K}{|}_{K=y},f_{S_i(T_j)}(y)=\frac{{\mathrm{\partial }}^2{V}_c(t_0,T_j,S_i(t_0),K))}{\mathrm{\partial }{K}^2}{|}_{K=y},\end{array}$(5) ) and $\mathbb{E}[B^2(T_j)]$ is defined in (17(17) $\begin{array}{r}\mathbb{E}[B^2(T_j)]=\sum _{i_1=1}^N\sum _{i_2=1}^N{\omega }_{i_1}{\omega }_{i_2}(\hat{\rho }(S_{i_1},S_{i_2}){\sigma }_{i_1}{\sigma }_{i_2}+\mathbb{E}[S_{i_1}(T_j)]\mathbb{E}[S_{i_2}(T_j)]).\end{array}$(17) ). It is important to note that, each change of the model parameters of the kernel process does impact parameters ${\alpha }_{i_1,j,1}$ and ${\alpha }_{i_2,j,1}$, implying that the mapping coefficients, ${\alpha }_{\cdot }$, need to be recomputed at each iteration. This however is extremely cheap operation as it only requires the CDF mappings presented in (10(10) $\begin{array}{rl}\log s_{i,j,k}& =F_{\log S_i(T_j)}^{-1}(F_{X_i(T_j)}(x_{i,j,k})),\end{array}$(10) ).

Following a similar strategy a third moment for the basket can be derived. The detailed derivations for the approximation results are presented in Appendix 2. By definition of the third moment we have: $\begin{array}{rl}\mathbb{E}[B^3]& =\sum _{i_1=1}^N\sum _{i_2=1}^N\sum _{i_3=1}^N{\omega }_{i_1,i_2,i_3}[\hat{\rho }(S_{i_1,i_2},S_{i_3}){\sigma }_{i_1,i_2}{\sigma }_{i_3}+\mathbb{E}[S_{i_1,i_2}]\mathbb{E}[S_{i_3}]]\\ & =\sum _{i_1=1}^N\sum _{i_2=1}^N\sum _{i_3=1}^N{\omega }_{i_1,i_2,i_3}[\hat{\rho }(S_{i_1,i_2},S_{i_3}){\sigma }_{i_1,i_2}{\sigma }_{i_3}+(\hat{\rho }(S_{i_1},S_{i_2}){\sigma }_{i_1}{\sigma }_{i_2}+\mathbb{E}[S_{i_1}]\mathbb{E}[S_{i_2}])\mathbb{E}[S_{i_3}]],\end{array}$with ${\omega }_{i_1,i_2,i_3}:={\omega }_{i_1}{\omega }_{i_2}{\omega }_{i_3}$, and $S_{i_1,i_2}:=S_{i_1}{S}_{i_2}$. The first moments $\mathbb{E}[S_{i_1}]$, $\mathbb{E}[S_{i_2}]$ and $\mathbb{E}[S_{i_3}]$ are known explicitly using the market market data for individual assets. Standard deviations, ${\sigma }_{i_1}$ and ${\sigma }_{i_2}$, of assets $S_{i_1}(T_j)$ and $S_{i_2}(T_j)$ can be computed utilizing the market implied density in (5(5) $\begin{array}{r}{F}_{S_i(T_j)}(y)=1+\frac{\mathrm{\partial }{V}_c(t_0,T_j,S_i(t_0),K))}{\mathrm{\partial }K}{|}_{K=y},f_{S_i(T_j)}(y)=\frac{{\mathrm{\partial }}^2{V}_c(t_0,T_j,S_i(t_0),K))}{\mathrm{\partial }{K}^2}{|}_{K=y},\end{array}$(5) ). Correlation, $\hat{\rho }(S_{i_1,i_2},S_{i_3})$, is approximated, as for the second moment, and yields: $\begin{array}{r}\hat{\rho }(S_{i_1,i_2},S_{i_3})\approx \frac{{\omega }_X(c_1,c_2,c_3)-{\omega }_X(c_1,c_2){\varphi }_{X_{i_3}}(-i{c}_3)}{\sqrt{({\omega }_X(2c_1,2c_2)-{\omega }_X^2(c_1,c_2))({\varphi }_{X_{i_3}}(-2i{c}_3)-{\varphi }_{X_{i_3}}^2(-i{c}_3))}},\end{array}$where: $c_1:={\hat{\alpha }}_{i_1,j,1}$, $c_2:={\hat{\alpha }}_{i_2,j,1}$, $c_3:={\hat{\alpha }}_{i_3,j,1}$ with ${\omega }_X(c_1,c_2,c_3)$ defined in (A5(A5) $\begin{array}{rl}{\omega }_X(a,b,c)& :=\mathbb{E}[{\mathrm{e}}^{a{X}_1+b{X}_2+c{X}_3}]\\ & =\exp (a{\mu }_1+b{\mu }_2+c{\mu }_3+\frac{1}{2}\sigma T+{\xi }_{p}T\mathbb{E}[{\mathrm{e}}^{\hat{J}}-1])\\ & =\exp [a{\mu }_1+b{\mu }_2+c{\mu }_3+\frac{1}{2}\hat{\sigma }T+{\xi }_{p}T({\mathrm{e}}^{(a+b+c){\mu }_J+\frac{1}{2}(a^2+b^2+c^2){\sigma }_J^2}-1)],\end{array}$(A5) ) and ${\omega }_X(c_1,c_2)$ in (21(21) $\begin{array}{rl}\mathbb{E}[{\mathrm{e}}^{a{X}_1(T)+b{X}_2(T)}|\mathcal{F}(t_0)]& =\exp [a{\mu }_1+b{\mu }_2+\frac{1}{2}{\hat{\sigma }}^{2}T+{\xi }_{p}T({\mathrm{e}}^{(a+b){\mu }_J+\frac{1}{2}(a^2+b^2){\sigma }_J^2}-1)]\\ & =:{\omega }_X(a,b),\end{array}$(21) ). Finally, the variance ${\sigma }_{i_1,i_2}^2$ is given explicitly by: (23) $\begin{array}{rl}{\sigma }_{i_1,i_2}^2& :=\mathbb{V}\mathrm{ar}[S_{i_1}{S}_{i_2}]=\mathbb{E}[S_{i_1}{S}_{i_2}]-\mathbb{E}[S_{i_1}]\mathbb{E}[S_{i_2}]\\ & =\hat{\rho }(S_{i_1},S_{i_2}){\sigma }_{i_1(T_j)}{\sigma }_{i_2(T_j)},\end{array}$(23) where $\hat{\rho }(S_{i_1},S_{i_2})$ is given in (18(18) $\begin{array}{rl}\hat{\rho }(S_{i_1},S_{i_2})& \approx \rho ({\mathrm{e}}^{{\hat{\alpha }}_{i_1,j,0}+{\hat{\alpha }}_{i_1,j,1}{X}_{i_1}(T_j)},{\mathrm{e}}^{{\hat{\alpha }}_{i_2,j,0}+{\hat{\alpha }}_{i_1,j,1}{X}_{i_2}(T_j)})\\ & =\rho ({\mathrm{e}}^{{\hat{\alpha }}_{i_1,j,1}{X}_{i_1}(T_j)},{\mathrm{e}}^{{\hat{\alpha }}_{i_2,j,1}{X}_{i_2}(T_j)})=:\overline{\rho }(S_{i_1},S_{i_2}),\end{array}$(18) ) and it is approximated by $\overline{\rho }(S_{i_1},S_{i_2})$ and derived in (19(19) $\begin{array}{r}\overline{\rho }(S_{i_1},S_{i_2})=\frac{\mathbb{E}[{\mathrm{e}}^{c_1{X}_{i_1}(T_j)+c_2{X}_{i_2}(T_j)}]-{\varphi }_{X_{i_1}}(-i{c}_1){\varphi }_{X_{i_2}}(-i{c}_2)}{\sqrt{({\varphi }_{X_{i_1}}(-2i{c}_1)-{\varphi }_{X_{i_1}}^2(-i{c}_1))({\varphi }_{X_{i_2}}(-2i{c}_2)-{\varphi }_{X_{i_2}}^2(-i{c}_2))}},\end{array}$(19) ). Availability of the third moment allows for an extended optimization problem proposed in (22(22) $\begin{array}{r}\underset{{\sigma }_J,{\mu }_J,{\xi }_p,\rho }{min}{(\mathbb{E}[B^2(T_j)]-\mathbb{E}[I^2(T_j)])}^2,\mathbb{E}[I^2(T_j)]={\int }_{\mathbb{R}}{x}^2{f}_{I(T_j)}(x)\mathrm{d}x,\end{array}$(22) ), i.e. (24) $\begin{array}{r}\underset{{\sigma }_J,{\mu }_J,{\xi }_p,\rho }{min}\sum _{k=2}^3{(\mathbb{E}[B^k(T_j)]-\mathbb{E}[I^k(T_j)])}^2,\mathbb{E}[I^k(T_j)]={\int }_{\mathbb{R}}{x}^k{f}_{I(T_j)}(x)\mathrm{d}x,\end{array}$(24) which will not only allow for the ATM fit but also the level of the implied volatility skew can be matched.

The quality of the approximating moment formulae will be discussed in the follow-up section (see Table 5), where a basket consisting of 5 assets will be considered.

4. Numerical results

This numerical section is dedicated to the numerical aspects and implementation details of the proposed model. We start with a step-by-step example outlining the procedure of constructing a 5-dimensional basket including the calibration of the collocation method to market-implied volatilities. We follow up with a discussion of the impact of model parameters on the resulting basket implied volatility and a finally conclude with a fully-fledged calibration of a high-dimensional basket to market data.

4.1. Illustrative 5D example

This section presents a step-by-step procedure for building up a basket using the ELV method. We illustrate in detail the complete process, starting from individual asset calibration to basket construction and the impact of the model parameters on the basket. In this example, we consider a basket consisting of 5 assets¹¹ that are a part of the DJIA 30. In the experiment, we take $T_j=1y$ as of 8/12/2021. The implied volatilities for this set of assets are presented in Table 1.

Table 1. Implied volatilities for 5 assets, for T = 1, observed on 8/12/2021.

asset	0.8	0.85	0.9	0.95	0.975	1	1.025	1.05	1.1	1.15	1.2
stock 1	32.20	31.00	29.98	29.07	28.65	28.25	27.87	27.51	26.85	26.28	25.83
stock 2	29.96	28.87	27.99	27.31	27.03	26.79	26.58	26.41	26.15	25.98	25.90
stock 3	32.23	31.07	30.11	29.34	29.02	28.73	28.48	28.26	27.90	27.65	27.47
stock 4	31.86	30.56	29.39	28.35	27.88	27.43	27.02	26.64	25.99	25.47	25.09
stock 5	34.58	33.84	33.25	32.76	32.55	32.35	32.17	32.00	31.71	31.46	31.26

Notes: Spot values for the underlying assets are: $S_1(t_0)=468.86$, $S_2(t_0)=411.25$, $S_3(t_0)=397.32$, $S_3(t_0)=334.97$, $S_5(t_0)=266.31$.

As the first step of our modelling procedure, we parametrize the market-implied volatilities using the SABR-based parametrization formula [16]. The parameters obtained in this step are presented in Table 2.

Table 2. Implied volatility parametrization using parameterization in [16].

asset	${\beta }_H$	${\alpha }_H$	${\rho }_H$	${\gamma }_H$
stock 1	0.90	0.52	−0.52	0.55
stock 2	0.90	0.48	−0.25	0.65
stock 3	0.90	0.51	−0.30	0.66
stock 4	0.90	0.49	−0.48	0.66
stock 5	0.90	0.56	−0.27	0.46

Notes: The meaning of the parameters is as follows: ${\beta }_H$ -‘exponent beta’, ${\alpha }_H$- ‘initial vol’, ${\rho }_H$- ‘correlation’ and ${\gamma }_H$- ‘vol-vol’. In the calibration procedure ${\beta }_H$ was fixed at 0.9. Interest rates were fixed at r = 0.0.

Once each implied volatility skew/smile is parameterized the next step is use (5(5) $\begin{array}{r}{F}_{S_i(T_j)}(y)=1+\frac{\mathrm{\partial }{V}_c(t_0,T_j,S_i(t_0),K))}{\mathrm{\partial }K}{|}_{K=y},f_{S_i(T_j)}(y)=\frac{{\mathrm{\partial }}^2{V}_c(t_0,T_j,S_i(t_0),K))}{\mathrm{\partial }{K}^2}{|}_{K=y},\end{array}$(5) ) and using the collocation method project the so-called ‘market distribution’ on $X_i$. We set the number of the collocation points to m = 7, and as discussed in Remark in Section 2.2 we use the grid stretching technique to control the tail behaviour. Firstly, we take $Z\sim \mathcal{N}(0,1)$ and we set the corresponding min/max quantiles at $q_1=0.01$ and $q_{m=7}=0.99$. Once the corresponding collocation points $z_k$, $k=1,\dots ,m$ are determined [14] the inverse for $x_{i,k}$ in (10(10) $\begin{array}{rl}\log s_{i,j,k}& =F_{\log S_i(T_j)}^{-1}(F_{X_i(T_j)}(x_{i,j,k})),\end{array}$(10) ) need to take place. Note that in principle each asset, $S_i(T_j)$, may be projected on the same variable $X_i.$ Even if $X_i$'s are correlated (via Brownian motion) it does not affect marginal distributions for $S_i$, but it will be crucial in the basket construction. Table 3 illustrates the details on computation of triples $(z_{i,k},x_{i,k},s_{i,k})\equiv (z_k,x_k,s_{i,k})$, $k=1,\dots ,7,$$i=1,\dots ,5.$

Table 3. Projection details for the collocation method.

	k = 1	k = 2	k = 3	k = 4	k = 5	k = 6	k = 7
$z_k$	−3.7504	−2.3668	−1.1544	0	1.154	2.3668	3.7504
$F_Z(z_k)$	0.0100	0.0710	0.237	0.5000	0.7630	0.929	0.9900
$x_k$	−0.1684	−0.1074	−0.0539	−0.0030	0.0478	0.1013	0.1623
$s_{1,k}$	4.9851	5.6027	5.9341	6.1545	6.3270	6.4891	6.6839
$s_{2,k}$	4.9591	5.5291	5.8174	6.0101	6.1734	6.3512	6.6100
$s_{3,k}$	4.7692	5.4420	5.7672	5.9771	6.1500	6.3350	6.6025
$s_{4,k}$	4.5586	5.2737	5.6126	5.8221	5.9828	6.1404	6.3508
$s_{5,k}$	4.4752	4.9965	5.3177	5.5592	5.7702	5.9855	6.2538

Notes: Parameters for $X_i$ were the following: $\sigma =0.1$, ${\sigma }_J=0.02$, ${\mu }_J=0.0$, ${\xi }_p=0.25$ and T = 0.5. Note that T here is considered as a parameter therefore it can be fixed for all the assets. In the experiment we have also taken the stretched grid with p = 0.01, as proposed in [14], implying the lower and upper bound for the CDF, $F_Z(z_1)$ and $F_Z(z_7)$.

Given the projection points presented in Table 3 and the underlying process $X_i(\cdot )$ we are able to replicate the marginal distribution functions of the assets observable in the market. By simulating a process $X_i(T_j)$ for every asset $i=1,\dots ,N$ with the parameters specified in Table 3 we evaluate the approximating polynomial for m = 7 with $T=1y,$ (25) $\begin{array}{r}\log S_i(T)\approx g_i^7(X_i(T))=\sum _{k=0}^{m-1}{s}_{i,k}{\ell }_k(x_k,X_i(T)),{\ell }_k(x_k,X_i(T))=\prod _{l=1,l\ne k}^m\frac{X_i(T)-x_l}{x_k-x_l},\end{array}$(25) with points $x_k$ and $s_{i,k}$ defined in Table 3 and $X_i(\cdot )$ defined in (8(8) $\begin{array}{r}{X}_i(t)=X_i(t_0)+(r-\frac{1}{2}{\sigma }_i^2-{\xi }_{p,i}\mathbb{E}[{\mathrm{e}}^{J_i}-1])t+{\sigma }_i{W}_i(t)+\sum _{k=1}^{X_{\mathcal{P},i}(t)}{J}_{k,i},\end{array}$(8) ). Ultimately, all the kernel processes $X_i$, $X_j$ need to be correlated with a given correlation coefficient. This, however, is not necessary if we are only interested in marginal distributions, thus in a fit of each asset to market implied volatility.

The calibration quality for the five assets is presented in Table 4. We report outstanding results; the error does not exceed 0.2%. This remarkable outcome is due to the core principle of the collocation method, which enables mapping a distribution of interest onto a polynomial of surrogate random variables using only a few iterations. As presented in ([14]), the method guarantees exponential convergence, even for challenging distributions.

Table 4. Table illustrates the calibration of the collocation method to market-implied volatilities.

		relative strike
asset	IV %	0.8	0.85	0.9	0.95	0.975	1	1.025	1.05	1.1	1.15	1.2
$S_1$	mkt	32.2	31.0	30.0	29.1	28.7	28.3	27.9	27.5	26.9	26.3	25.8
	coll	32.1	31.0	30.0	29.0	28.6	28.2	27.8	27.4	26.8	26.2	25.7
	err	0.1	0.0	0.0	0.1	0.1	0.1	0.1	0.1	0.1	0.1	0.1
$S_2$	mkt	29.9	28.9	28.0	27.3	27.0	26.8	26.6	26.4	26.1	26.0	25.9
	coll	29.9	28.9	28.0	27.3	27.0	26.8	26.5	26.4	26.1	25.9	25.9
	err	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.1
$S_3$	mkt	32.2	31.1	30.2	29.4	29.0	28.7	28.5	28.2	27.9	27.6	27.5
	coll	32.0	31.0	30.0	29.3	28.9	28.7	28.4	28.2	27.8	27.6	27.4
	err	0.2	0.1	0.1	0.1	0.1	0.1	0.1	0.1	0.1	0.1	0.1
$S_4$	mkt	31.9	30.6	29.4	28.4	27.9	27.4	27.0	26.6	26.0	25.5	25.1
	coll	31.7	30.4	29.2	28.2	27.7	27.3	26.9	26.5	25.9	25.4	25.0
	err	0.2	0.2	0.2	0.1	0.1	0.1	0.1	0.1	0.1	0.1	0.1
$S_5$	mkt	34.5	33.9	33.3	32.8	32.6	32.4	32.2	32.0	31.7	31.5	31.3
	coll	34.6	33.9	33.3	32.8	32.6	32.4	32.2	32.0	31.7	31.5	31.3
	err	−0.1	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0

Notes: ‘mkt’ stands for the market implied volatility obtained from the parameterization of implied volatility, ‘coll’ represents the results obtained from the collocation method and ‘err’ is the difference of implied volatilities, in $\mathrm{\%}$.

As the final step of this illustrative example we construct a basket, $B(T)=\sum _{i=1}^5{S}_i(T)$ with $S_i(t)$ given in (25(25) $\begin{array}{r}\log S_i(T)\approx g_i^7(X_i(T))=\sum _{k=0}^{m-1}{s}_{i,k}{\ell }_k(x_k,X_i(T)),{\ell }_k(x_k,X_i(T))=\prod _{l=1,l\ne k}^m\frac{X_i(T)-x_l}{x_k-x_l},\end{array}$(25) ). For this purpose, we utilize Equation (14(14) $\begin{array}{rl}B(T_j)& =\sum _{i=1}^N{S}_i(T_j)\approx \sum _{i=1}^N\exp (g_{i,j}^m(X_i(T_j)))\\ & =\sum _{i=1}^N\exp (\sum _{k=0}^{m-1}{\hat{\alpha }}_{i,j,k}{X}_i^k(T_j)),\end{array}$(14) ), where the processes $X_i$ need to be simulated with Monte Carlo. As discussed before, we consider a standard Poisson process, $X_{\mathcal{P}}(T)$, for all the underlying assets and correlated Brownian Motions, $W(t).$ These processes can be pre-simulated and then fed into the asset Equation (25(25) $\begin{array}{r}\log S_i(T)\approx g_i^7(X_i(T))=\sum _{k=0}^{m-1}{s}_{i,k}{\ell }_k(x_k,X_i(T)),{\ell }_k(x_k,X_i(T))=\prod _{l=1,l\ne k}^m\frac{X_i(T)-x_l}{x_k-x_l},\end{array}$(25) ). This procedure will produce Monte Carlo paths for the basket, $B(t)$.

The next step is to determine optimal kernel model parameters. The objective is to find the parameters that will generate basket implied volatility as close as possible to the implied volatility of the index. In this illustrative experiment, we do not consider an index; however, we can analyse the quality of the approximating formulae and the impact of the model parameters on the basket volatilities.

In the calibration procedure, we will follow Equation (24(24) $\begin{array}{r}\underset{{\sigma }_J,{\mu }_J,{\xi }_p,\rho }{min}\sum _{k=2}^3{(\mathbb{E}[B^k(T_j)]-\mathbb{E}[I^k(T_j)])}^2,\mathbb{E}[I^k(T_j)]={\int }_{\mathbb{R}}{x}^k{f}_{I(T_j)}(x)\mathrm{d}x,\end{array}$(24) ), which relies on fitting based on the variance and the skew. In Table 5 the numerical results for estimated variance and the skew for varying model parameters are shown. The results are excellent, especially for the standard deviation, ${\sigma }_B$, where the error is about 0.1.

Table 5. Table illustrates the quality of the approximation of the moments derived for a considered basket $B(T)$, T = 1, and varying set of the model parameters.

parameters					Monte Carlo		Approximation
σ	${\sigma }_J$	${\mu }_J$	${\xi }_p$	ρ	${\sigma }_B$	$\mathbb{E}[B^3]/{\sigma }_B^3$	${\sigma }_B$	$\mathbb{E}[B^3]/{\sigma }_B^3$
0.2	0.02	0.0	0.25	0.5	28.17	98.42	28.1	90.84
0.2	·	·	·	·	28.19	98.17	28.17	90.17
0.5	·	·	·	·	28.18	98.25	28.19	89.95
·	0.01	·	·	·	28.19	98.21	28.18	90.02
·	0.04	·	·	·	27.88	101.25	27.68	94.85
·	·	−0.1	·	·	29.77	84.48	29.61	78.13
·	·	−0.25	·	·	31.18	74.23	31.86	63.47
·	·	·	0.1	·	28.2	98.11	28.17	90.13
·	·	·	0.15	·	28.19	98.22	28.14	90.47
·	·	·	·	0.1	19.44	277.62	19.33	269.78
·	·	·	·	0.8	33.31	62.87	33.31	56.03

Notes: The first 5 columns of the table illustrate particular parameter configuration (dots indicate unchanged parameters). ${\sigma }_B:=\sqrt{\mathbb{V}\mathrm{ar}[B(T)]}$.

In the final part of this section, we illustrate, in Table 6, the timing results. As presented in the table, the most time-consuming part is to parameterize market discrete implied volatilities into a parametric form. Typically, such a process is performed only once, and only the relevant, optimized parameters are stored. The collocation method's calibration and the collocation points' generation is a considerably cheap operation. Finally, the Monte Carlo simulation required sampling of all the underlying processes and construction of the basket.

Table 6. Timing results in seconds for the basket consisting of 5 assets.

IV param. (sec)	Gen.$x_i$ (sec)	Gen.$y_i$ (sec)	Monte Carlo (sec)
0.5765	0.040	0.0876	0.240

Notes: The results were performed on a standard home PC where no parallelization was used. Timings are aggregated for all the assets and are reported in seconds. ‘IV param’ indicates the time needed to parameterize the implied volatilities- this process may be considered irrelevant for this framework as it is a part of data processing. ‘Gen.$x_i$’ indicates the time needed to generate the collocation points $z_i$ and $x_i$ based on the process $X_i$. ‘Gen.$y_i$’ represents the time needed to calibrate the collocation method via the inverse of the corresponding CDFs. ‘Monte Carlo’ corresponds to the complete basket simulation with five assets and 50k paths for each stock.

In the next section the impact of the model parameters on basket implied volatility will be studied.

4.2. Impact of the model parameters on basket implied volatility

This section analyses the impact of the kernel model parameters on the basket implied volatilities induced by the ELV model. As described earlier, due to the collocation method, a particular choice of model parameters does not have a material impact on the fit of individual assets- they, due to the mapping procedure, remain intact. However, it has an impact on the covariance structure of the basket. This effect can be measured in terms of implied volatilities. In the experiment we set the parameters of process $X(t)$ (all the parameter values are the same for all the assets) $\sigma =0.10$, ${\sigma }_J=0.02$, ${\mu }_J=0$, ${\xi }_p=0.25$, T = 0.5 with the correlation between all the Brownian motions, $\mathrm{d}{W}_i(t)\mathrm{d}{W}_j(t)=\rho \mathrm{d}t$, $\rho =0.5$ and measure the impact on the implied volatility of the basket, $B(T)$- the basket consists of 5 assets and it's construction is defined in Section 4.1.

Figure 4. Impact of different jump parameters on the shape of the implied volatility in Merton's jump diffusion model. First figure: impact of ${\sigma }_J$; Second figure: impact of ${\xi }_p$; Third figure: impact of ${\mu }_J.$

s

Figures 5 and  6 show that¹²: σ, ${\sigma }_J$ and ${\xi }_p$ have a moderate impact on the level of basket implied volatilities. The two remaining parameters are the most relevant: ${\mu }_J$ and ρ. To a large extent, correlation, ρ determines the level of basket implied volatilities (higher correlation lower the implied volatility) while ${\mu }_J$ in a significant way controls the implied volatility skew, as demonstrated in Figure 7. We considered negative values for ${\mu }_J$ in the experiment, implying negative shocks to the underlying process. We observe that a more negative value of ${\mu }_J$ produces more skew.

Figure 5. Impact of different jump parameters on the shape of the basket implied volatility in the ELV model. First figure: impact of ${\sigma }_J$; Second figure: impact of ${\mu }_J$.

Figure 6. Impact of different jump parameters on the shape of the basket implied volatility in the ELV model. First figure: impact of ${\xi }_p$; Second figure: impact of ρ- correlation between Brownian motions.

Figure 7. Basket implied volatilities with re-calibrated model parameters for ATM level.

Since both model parameters ${\sigma }_J$ and ${\xi }_p$ impact the basket implied volatility in a similar fashion, this means that in the calibration routine, one of them can be effectively fixed, reducing the complexity of the calibration procedure. First, the calibration process can be performed in an iterative way, where the ${\mu }_J$ (skew) is chosen; secondly, the ATM implied volatility is found such that the basked and index ATM volatility matched perfectly.

The presented results suggest that the model can only generate implied volatility skew. The basked implied volatility shapes may also be affected by the basket's composition. In Figure 4, we have shown that the kernel process allows for the generation of a smile, though limited. On the other hand, each asset's volatility contributes to the overall basket vol. This will be confirmed in the numerical experiment where we will consider a DJIA index with 30 underlying assets.

4.3. Calibration of a high dimensional basket

In this section, we perform a practical experiment where we build a bridge between a basket and an index. This section considers a DJIA (Dow Jones Industrial Average) with 30 underlying assets, as of 12/08/2021, with the spot price of 357.55. This experiment poses an extension of a case considered in Section 4.1 where a 5D case was examined. Now, however, we will also perform the model calibration, allowing us for consistent pricing of both an index and a combination of basket constituents.

The essence of the ELV model lies in the construction of a local volatility model that is constructed of basket marginal distributions that are consistent with the index. As such, the calibration procedure can be performed separately for each expiry date $T_i$, where we utilize analytic expression for the basket's moments (see Equation (24(24) $\begin{array}{r}\underset{{\sigma }_J,{\mu }_J,{\xi }_p,\rho }{min}\sum _{k=2}^3{(\mathbb{E}[B^k(T_j)]-\mathbb{E}[I^k(T_j)])}^2,\mathbb{E}[I^k(T_j)]={\int }_{\mathbb{R}}{x}^k{f}_{I(T_j)}(x)\mathrm{d}x,\end{array}$(24) ) that are matched with the moments computed from the index. The calibration is straightforward, and as indicated earlier, there are two dominant model parameters, ${\mu }_J$ and ρ, that play a central role in the skew calibration. Given the closed-form formulas for the moments, the calibration is swift and straightforward and can be performed with a simple search algorithm.

In the scope of this experiment, we undertook the calibration of 30 assets from the DJIA index, employing Equation (25(25) $\begin{array}{r}\log S_i(T)\approx g_i^7(X_i(T))=\sum _{k=0}^{m-1}{s}_{i,k}{\ell }_k(x_k,X_i(T)),{\ell }_k(x_k,X_i(T))=\prod _{l=1,l\ne k}^m\frac{X_i(T)-x_l}{x_k-x_l},\end{array}$(25) ) as the underlying methodology. The calibration quality for selected assets is presented in Figure 8. Subsequently, leveraging the unconstrained parameters inherent in the model, ${\sigma }_M$, ${\sigma }_J$, ${\mu }_J$, ${\xi }_P$, and ρ, we proceeded to calibration of the index options basket. The calibration covered a span of seven distinct maturities: 1m, 2m, 3m, 6m, 1y, 1.5y, and 2y. In Table 7, the estimated model parameters are presented, while the moments of the basket compared to the index and ‘exact’ moments obtained from Monte Carlo are tabulated in Table 8.

Figure 8. Calibration quality for four selected stocks. Red squared indicate market quote, blue line corresponds to the collocation method.

Table 7. Calibrated model parameters.

$T_i$	${\sigma }_M$	${\sigma }_J$	${\mu }_J$	${\xi }_P$	T	ρ
$1m$	0.25	0.020	−0.20	0.25	0.30	0.35
$2m$	0.15	0.020	−0.20	0.25	0.50	0.25
$3m$	0.15	0.020	−0.20	0.25	0.50	0.25
$6m$	0.15	0.020	−0.20	0.25	0.50	0.37
$1y$	0.1	0.020	−0.20	0.25	0.50	0.45
1.5y	0.1	0.020	−0.20	0.25	0.50	0.45
$2y$	0.20	0.020	−0.25	0.50	0.45	0.38

Table 8. Moment approximation quality for DJIA with 30 underlying assets, as defined in Table 5.

	Index (Market)		Basket (Monte Carlo)		Basket (Proxy)
$T_i$	${\sigma }_B$	$\mathbb{E}[B^3]/{\sigma }_B^3$	${\sigma }_B$	$\mathbb{E}[B^3]/{\sigma }_B^3$	${\sigma }_B$	$\mathbb{E}[B^3]/{\sigma }_B^3$
$1m$	19.08	6,634.4	19.0	6,714.1	18.7	6,982.6
$2m$	29.07	1,895.9	29.3	1,858.0	28.6	1,958.6
$3m$	36.47	969.8	36.0	1,009.5	35.2	1,059.5
$6m$	53.34	319.9	54.3	304.0	53.7	302.8
$1y$	78.56	107.3	81.9	95.9	82.8	85.8
1.5y	94.52	65.3	99.0	58.0	99.7	50.6
$2y$	109.07	45.4	113.5	41.3	113.1	35.5

Further insights into the calibration's effectiveness concerning an index are shown in Figures 9 and 10. We report a great alignment across all considered option expires, i.e. the ELV model (red line) closely fits the mid-market implied volatility (blue line) and is consistently within the bounds of the bid-ask spread. We report the maximum error of 0.2% in terms of the implied volatilities for constituent assets and about 1% in the case of the basket option. It is important to stress that a 1% error is within the bid-ask spread. It is critical to highlight that incorporating jumps within the ELV model's construct contributes to the more pronounced skew in implied volatility, i.e. the model without jumps corresponds to a standard representation of local volatility depicted by a purple line. The calibration highlights this assertion results in Table 7, where negative ${\mu }_J$ values were obtained for all option expiries.

Figure 9. Implied volatilities for DJIA 30 compared to ELV model with and without jumps; LHS: T = 2M, RHS: T = 3M.

Figure 10. Implied volatilities for DJIA 30 compared to ELV model with and without Jumps; LHS: T = 1y, RHS: T = 2y.

The results show the ELV model's remarkable performance in effectively capturing individual stocks' distributions and their underlying covariance structure. The model's fit to the bid-ask spreads across all available index quotations is highly desired.

4.4. Cost of basket miscalibration

In this section, we assess the impact of basket miscalibration. In the experiment, we will continue with the 30-dimensional basket calibrated in Section 4.3. We will consider a subset of assets $S_{i_1},\dots S_{i_K}$ of the basket and price an exotic derivative using two variants of the ELV framework. In the first case, we take the model with a covariance structure only driven by Brownian motion, while the alternative is the model where jumps are included. Our previous experiment has shown that including jumps facilitates a better fit to the index. Now, we quantify the potential cost of this mispricing.

In the first step of pricing, the local volatility surface for the ELV dynamics, $\overline{B}$ in (2(2) $\begin{array}{r}\mathrm{d}\overline{B}(t)=r\overline{B}(t)\mathrm{d}t+{\sigma }_{LV}(t,\overline{B}(t))\overline{B}(t)\mathrm{d}W(t),\end{array}$(2) ), needs to be constructed. This is performed based on the calibrated basket in Equation (4(4) $\begin{array}{r}{S}_i(T_j)\approx {\mathrm{e}}^{g_{i,j}(X_i(T_j))},\mathrm{and}B(T_j)\approx \sum _{i=1}^N{\omega }_i{\mathrm{e}}^{g_{i,j}(X_i(T_j))}.\end{array}$(4) ). The number of expiries, $N_T$, available for building the local volatility surface depends on the required accuracy and discretization time-step of (2(2) $\begin{array}{r}\mathrm{d}\overline{B}(t)=r\overline{B}(t)\mathrm{d}t+{\sigma }_{LV}(t,\overline{B}(t))\overline{B}(t)\mathrm{d}W(t),\end{array}$(2) ). Although, in practice, a finite and somewhat limited set of expiries, $T_j$, is available in the market, for the intermediate expiry, T, $T_j<T<T_{j+1}$ we can impose a smooth transition between the corresponding model parameters, i.e. $\begin{array}{r}\mathbf{p}(T)=\frac{T_{j+1}-T}{T_{j+1}-T_j}\mathbf{p}(T_j)+\frac{T-T_j}{T_{j+1}-T_j}\mathbf{p}(T_{j+1}),\end{array}$where $\mathbf{p}(T_j)$, $\mathbf{p}(T_{j+1})$ represent the parameters of the kernel process $X(T_j)$, $X(T_{j+1})$ observed at times $T_j$ and $T_{j+1}$, respectively.

In this experiment, we consider a basket consisting of 5 stocks, already introduced earlier in Section 4.1 where the subset of DJIA stocks has been considered, and two exotic derivatives, namely, an arithmetic Asian option, defined as: $\begin{array}{r}{V}_A(t_0)=N{\mathrm{e}}^{-rT}\mathbb{E}[{(\sum _{l=1}^{L}S(T_l)-K_A)}^{+}|\mathcal{F}(t_0)],\end{array}$with N being the notional and where $K_A$ is the strike price, and the up-and-out barrier option: $\begin{array}{r}{V}_{Bar}(t_0)=N{\mathrm{e}}^{-rT}\mathbb{E}[(B(T)-K_B{)}^{+}{1}_{\tau (H,B)>T}|\mathcal{F}(t_0)],\end{array}$where $K_B$ is the strike price, $H>B(t_0)$ is the barrier, and $\tau (x,B):=inf\{t\ge 0:B(t)>x\}$, and the notional, N. In the experiment we consider the ATM option with $K_B=123.82$ and the notional of 1. Although relatively standard in the financial world, both derivatives are sensitive to different characteristics of the underlying model.

The pricing for varying strikes and barrier levels is presented in Figure 11. The results show that although in the case of the Asian option, the skew-related miscalibration does not have a significant impact in the case of the up-and-out barrier option, the skew plays a vital role. As illustrated in Table 9 the price differences may be even ten-fold.

Figure 11. Figure illustrates prices of Asian and Barrier option prices for: LHS varying strike, K, RHS: varying the barrier level, H.

Table 9. Comparison of Asian and Barrier option pricing for ELV models with and without jumps.

Asian option price differences
K	100.05	108.46	116.87	125.28	133.69	142.10	150.51
diff	0.0392	−0.0555	−0.1439	−0.0464	0.0870	0.0602	0.0289
Barrier option price differences
H	123.8225	178.1855	232.5484	286.9113	341.2742	395.6371	450
diff	0	−1.881	−1.2491	−0.4639	−0.1986	−0.0891	−0.0735

5. Conclusions

We have proposed above the Effective Local Volatility model capable of accurate, consistent and swift pricing of basket options – both plain vanilla and of the more exotic type. The model is accurate, because – as demonstrated in Section 4.3 – it calibrates to liquid index options with only minimal error, well within the bid-ask spread. The model is consistent, because it simultaneously calibrates to individual basket constituents (via the functions $g_{i,j}(\cdot )$) and their index-implied covariance structure (via the kernel processes $X_i(\cdot )$ which we assume to be of Merton's jump diffusion type). Finally, the model works swiftly, because the typically numerically ‘expensive’ task, like sampling from stocks' marginal distributions, is handled via the collocation method, and basket calibration is performed by means of moment matching technique, leveraging the analytical expressions for basket moments which we derive in Section 3. Importantly, unlike copula-based or local-correlation models which tend to be black-boxy and intractable, the proposed approach allows users to control the basket implied volatility skew using just a few parameters in a clear and transparent way. Thus, it is straightforward to control the dependence structure of the basket and trace how implied volatility smiles transfer from individual assets to the whole basket (or possibly the index). While interesting in its own right, this feature makes it possible to price structured products including baskets on any combination of index constituents stocks: the assets simply ‘inherit’ the covariance structure implied by the index which we then feed into a one-dimensional local volatility model, allowing us to price a desired structured product in a quick and easy way.

While the model as currently proposed strikes a good balance between accuracy and tractability and works very well for basket dimensions relevant in practice, the approach can nonetheless be extended in two main directions. First, one could postulate more sophisticated dynamics for the kernel processes than the simple jump diffusion model. This might allow for even better control of the implied volatility shapes, however possibly at some cost to tractability, slowing down the calibration. Second, one could add a stochastic component to the Effective Local Volatility model which might allow for more accurate pricing of callable structures.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Notes

1 Cf. the prospectus availible in the online records of the U.S. Securities and Exchange Commission at: https://www.sec.gov/Archives/edgar/data/19617/000089109221003578/e13291-424b2.htm.

2 For example, in the Bloomberg basket options pricing template correlations are, by default, estimated over a 5 year period, whereby to eliminate noise, a given percentile of rolling 6-month cross-correlation estimates is chosen in the parameterization of the full correlation matrix.

3 We define the skew here loosely as the difference in implied volatilities between the 85 and 120$\mathrm{\%}$ ATM levels.

4 As an alternative to [27] one could consider Kou's jump-diffusion model [18] which has the additional benefit of separating the upside and downside skew. However, in this case, we opt for the simplicity and parsimony of Merton's approach.

5 Without loss of generality, we shall henceforth think of the underlying assets as stocks, however, the method developed below is obviously general and, mutatis mutandis, applies to other instruments as well.

6 The proposed framework can also be extended with a stochastic volatility process. Such an extension is trivial and will, for simplicity, be omitted.

7 The respective dynamics are given by (j = 1, 2): $\mathrm{d}{S}_j(t)=r{S}_j(t)\mathrm{d}t+v_j^{1/2}(t)S_j(t)\mathrm{d}{W}_{j,1}(t)$, $\mathrm{d}{v}_j(t)={\kappa }_j({\overline{v}}_j-v_j(t))\mathrm{d}t+{\gamma }_j{v}_j^{1/2}(t)\mathrm{d}{W}_{j,2}(t)$ with correlations $\mathrm{d}{W}_{j,1}(t)\mathrm{d}{W}_{j,2}(t)={\rho }_j\mathrm{d}t$, $\mathrm{d}{W}_{1,1}(t)\mathrm{d}{W}_{2,1}(t)={\rho }_{1,2}\mathrm{d}t$ and $\mathrm{d}{W}_{j,2}(t)\mathrm{d}{W}_{k,2}(t)=0\cdot \mathrm{d}t$. For reference, we set $S_1(t_0)=1,$ $S_2(t_0)=2.5,$ $r=0,$ ${\kappa }_1=1,$ ${\kappa }_2=0.5,$ ${\gamma }_1=1,$ ${\gamma }_2=0.6,$ ${\rho }_{S_1,v_1}=-0.5$, ${\rho }_{S_2,v_2}=-0.7$, $v_{1,0}=0.1$, $v_{2,0}=0.05$, ${\overline{v}}_1=0.1$ and ${\overline{v}}_2=0.05$.

8 The Feller's condition is a direct consequence of the so-called Fichera [11] condition for the uniqueness of solutions to elliptic and parabolic equations having diffusion coefficients vanishing on the boundary of the computational domain. It gives necessary and sufficient conditions for advection terms guaranteeing the unicity of solutions.

9 The reason why we choose a standard normal distribution in the alternative approach is twofold. First, even for a fundamental distribution as the standard normal results are highly accurate – this is also the case in e.g. [12]. By choosing a different distribution, results may be further enhanced. Secondly, as mentioned in [14], choosing the normal distribution is also motivated by the Cameron–Martin Theorem [30], which states that polynomial chaos approximations based on the normal distribution converge to any distribution.

10 The strategy proposed in this part does not require ‘re-calibration’ of ${\alpha }_{i,j,k}$ coefficients, but only neglects the coefficients of higher order.

11 (1) UnitedHealth; (2) Home Depot; (3) Goldman Sachs; (4) Microsoft Corp; (5) salesforce.com Inc.

12 Results for σ are not presented here as they resembled the impacts of ${\xi }_p$ and ${\sigma }_J$.

Appendices

Appendix 1.

Proof of Lemma 3.1

Proof.

By substitution we find: $\begin{array}{r}\mathbb{E}[{\mathrm{e}}^{a{X}_1(T)+b{X}_2(T)}]={\mathrm{e}}^{a{\mu }_1+b{\mu }_2}\mathbb{E}[\exp (a{\sigma }_1{W}_1(T)+b{\sigma }_2{W}_2(T)+\sum _{k=1}^{X_{\mathcal{P}}(T)}a{J}_{k,1}+\sum _{k=1}^{X_{\mathcal{P}}(T_j)}b{J}_{k,2})],\end{array}$with $X_1$, $X_2$ defined in (20(20) $\begin{array}{r}{X}_1(t)={\mu }_1+{\sigma }_1{W}_1(t)+\sum _{k=1}^{X_{\mathcal{P}}(t)}{J}_{k,1},X_2(t)={\mu }_2+{\sigma }_2{W}_2(t)+\sum _{k=1}^{X_{\mathcal{P}}(t)}{J}_{k,2},\end{array}$(20) ), and where $J_{k,i}\sim \mathcal{N}({\mu }_J,{\sigma }_J^2)$. Due to the correlation between Brownian motions we perform the following factorization, $\begin{array}{r}a{\sigma }_1{W}_1(T)+b{\sigma }_2{W}_2(T)\stackrel{\mathrm{d}}{=}\sqrt{a^2{\sigma }_1^2+b^2{\sigma }_2^2+2{\rho }_{1,2}ab{\sigma }_1{\sigma }_2}\tilde{W}(T)=:\hat{\sigma }\tilde{W}(T),\end{array}$with ${\hat{\sigma }}^2=a^2{\sigma }_1^2+b^2{\sigma }_2^2+2{\rho }_{1,2}ab{\sigma }_1{\sigma }_2$ and where $\tilde{W}(T)$ is an independent Brownian motion. $\begin{array}{r}\sum _{k=1}^{X_{\mathcal{P}}(T)}a{J}_{k,1}+\sum _{k=1}^{X_{\mathcal{P}}(T)}b{J}_{k,2}=\sum _{k=1}^{X_{\mathcal{P}}(T)}(a{J}_{k,1}+b{J}_{k,2})=:\sum _{k=1}^{X_{\mathcal{P}}(T)}{\hat{J}}_k,\end{array}$where $\hat{J}\equiv {\hat{J}}_k=a{J}_{k,1}+b{J}_{k,2}\sim \mathcal{N}((a+b){\mu }_J,(a^2+b^2){\sigma }_J^2)$; thus we have: $\begin{array}{r}\mathbb{E}[{\mathrm{e}}^{a{X}_1(T)+b{X}_2(T)}]=\mathbb{E}[\exp (\hat{\sigma }\tilde{W}(T)+\sum _{k=1}^{X_{\mathcal{P}}(T)}{\hat{J}}_k)]\stackrel{\perp \perp }{=}{\mathrm{e}}^{a{\mu }_1+b{\mu }_2}\mathbb{E}\left[{\mathrm{e}}^{\hat{\sigma }\tilde{W}(T)}\right]\mathbb{E}\left[{\mathrm{e}}^{\sum _{k=1}^{X_{\mathcal{P}}(T)}{\hat{J}}_k}\right].\end{array}$The first expectation simply reads $\mathbb{E}[{\mathrm{e}}^{\hat{\sigma }\tilde{W}(T)}]={\mathrm{e}}^{\frac{1}{2}{\hat{\sigma }}^{2}T},$ while the second expectation, thanks to the tower property of expectations, gives: $\begin{array}{rl}\mathbb{E}[\exp \left(\sum _{k=1}^{X_{\mathcal{P}}(T)}{\hat{J}}_k\right)]& =\sum _{n\ge 0}\mathbb{E}[\exp \left(\sum _{k=1}^n{\hat{J}}_k\right)]\frac{{\mathrm{e}}^{-{\xi }_{p}T({\xi }_{p}T{)}^n}}{n!}\\ & =\sum _{n\ge 0}\frac{{\mathrm{e}}^{-{\xi }_{p}T({\xi }_{p}T{)}^n}}{n!}{\mathbb{E}}^n\left[{\mathrm{e}}^{\hat{J}}\right]\\ & =\exp ({\xi }_{p}T\mathbb{E}[{\mathrm{e}}^{\hat{J}}-1]),\end{array}$where the last expression can be recognized as a Taylor expansion of an exponential function. Since $\begin{array}{r}\mathbb{E}[{\mathrm{e}}^{\hat{J}}]={\mathrm{e}}^{(a+b){\mu }_J+\frac{1}{2}(a^2+b^2){\sigma }_J^2},\end{array}$by substitutions the proof is complete.

Appendix 2.

Third moment derivations

For convenience we neglect the time argument $T_j$. The third moment of the basket, by definition, is defined as: $\begin{array}{r}\mathbb{E}[B^3]=\sum _{i_1=1}^N\sum _{i_2=1}^N\sum _{i_3=1}^N{\omega }_{i_1}{\omega }_{i_2}{\omega }_{i_3}\mathbb{E}\left[S_{i_1}{S}_{i_2}{S}_{i_3}\right]=:\sum _{i_1=1}^N\sum _{i_2=1}^N\sum _{i_3=1}^N{\omega }_{i_1}{\omega }_{i_2}{\omega }_{i_3}\mathbb{E}\left[S_{i_1,i_2}{S}_{i_3}\right],\end{array}$where $S_{i_1,i_2}=S_{i_1}{S}_{i_2}$. Following the same procedure as for the second moment we find: $\begin{array}{r}\mathbb{E}[B^3]=\sum _{i_1=1}^N\sum _{i_2=1}^N\sum _{i_3=1}^N{\omega }_{i_1}{\omega }_{i_2}{\omega }_{i_3}[\hat{\rho }(S_{i_1,i_2},S_{i_3}){\sigma }_{S_{i_1,i_2}}{\sigma }_{S_{i_2}}+\mathbb{E}[S_{i_1,i_2}]\mathbb{E}[S_{i_3}]].\end{array}$$\mathbb{E}[S_{i_1,i_2}]:=\mathbb{E}[S_{i_1}{S}_{i_2}]$, can be computed with (16(16) $\begin{array}{r}\mathbb{E}[S_{i_1}(T_j)S_{i_2}(T_j)]=\hat{\rho }(S_{i_1},S_{i_2}){\sigma }_{i_1}{\sigma }_{i_2}+\mathbb{E}[S_{i_1}(T_j)]\mathbb{E}[S_{i_2}(T_j)],\end{array}$(16) ), $\mathbb{E}[S_{i_3}]$ is given from the market data. Terms we need to establish are: ${\sigma }_{S_{i_1,i_2}}$ and $\hat{\rho }(S_{i_1,i_2},S_{i_3})$. For the first term we have: $\begin{array}{r}{\sigma }_{S_{i_1,i_2}}^2=\mathbb{E}[S_{i_1}{S}_{i_2}]-\mathbb{E}[S_{i_1}]\mathbb{E}[S_{i_2}],\end{array}$which again can be computed using Equation (16(16) $\begin{array}{r}\mathbb{E}[S_{i_1}(T_j)S_{i_2}(T_j)]=\hat{\rho }(S_{i_1},S_{i_2}){\sigma }_{i_1}{\sigma }_{i_2}+\mathbb{E}[S_{i_1}(T_j)]\mathbb{E}[S_{i_2}(T_j)],\end{array}$(16) ). For the correlation coefficient we have: $\begin{array}{rl}\hat{\rho }(S_{i_1,i_2},S_{i_3})& =\rho ({\mathrm{e}}^{{\hat{\alpha }}_{i_1,j,0}+{\hat{\alpha }}_{i_1,j,1}{X}_{i_1}}{\mathrm{e}}^{{\hat{\alpha }}_{i_2,j,0}+{\hat{\alpha }}_{i_1,j,1}{X}_{i_2}},{\mathrm{e}}^{{\hat{\alpha }}_{i_3,j,0}+{\hat{\alpha }}_{i_3,j,1}{X}_{i_3}})\\ & =\rho ({\mathrm{e}}^{{\hat{\alpha }}_{i_1,j,1}{X}_{i_1}+{\hat{\alpha }}_{i_2,j,1}{X}_{i_2}},{\mathrm{e}}^{{\hat{\alpha }}_{i_3,j,1}{X}_{i_3}})=:\overline{\rho }(S_{i_1,i_2},S_{i_3}),\end{array}$which yields: $\begin{array}{rl}\overline{\rho }(S_{i_1,i_2},S_{i_3})& \stackrel{\mathrm{def}}{=}\frac{\mathbb{E}[{\mathrm{e}}^{{\hat{\alpha }}_{i_1,j,1}{X}_{i_1}+{\hat{\alpha }}_{i_2,j,1}{X}_{i_2}+{\hat{\alpha }}_{i_3,j,1}{X}_{i_3}}]-\mathbb{E}[{\mathrm{e}}^{{\hat{\alpha }}_{i_1,j,1}{X}_{i_1}+{\hat{\alpha }}_{i_2,j,1}{X}_{i_2}}]\mathbb{E}[{\mathrm{e}}^{{\hat{\alpha }}_{i_3,j,1}{X}_{i_3}}]}{\sqrt{\mathbb{V}\mathrm{ar}[{\mathrm{e}}^{{\hat{\alpha }}_{i_1,j,1}{X}_{i_1}+{\hat{\alpha }}_{i_2,j,1}{X}_{i_2}}]\mathbb{V}\mathrm{ar}[{\mathrm{e}}^{{\hat{\alpha }}_{i_3,j,1}{X}_{i_3}}]}}\\ & =\frac{\mathbb{E}[{\mathrm{e}}^{{\hat{\alpha }}_{i_1,j,1}{X}_{i_1}+{\hat{\alpha }}_{i_2,j,1}{X}_{i_2}+{\hat{\alpha }}_{i_3,j,1}{X}_{i_3}}]-\mathbb{E}[{\mathrm{e}}^{{\hat{\alpha }}_{i_1,j,1}{X}_{i_1}+{\hat{\alpha }}_{i_2,j,1}{X}_{i_2}}]{\varphi }_{X_{i_3}}(-i{\hat{\alpha }}_{i_3,j,1})}{\sqrt{\mathbb{V}\mathrm{ar}[{\mathrm{e}}^{{\hat{\alpha }}_{i_1,j,1}{X}_{i_1}+{\hat{\alpha }}_{i_2,j,1}{X}_{i_2}}]({\varphi }_{X_{i_3}}(-2i{\hat{\alpha }}_{i_3,j,1})-{\varphi }_{X_{i_3}}^2(-i{\hat{\alpha }}_{i_3,j,1}))}}.\end{array}$All the expression except for the first expectation can be computed using (16(16) $\begin{array}{r}\mathbb{E}[S_{i_1}(T_j)S_{i_2}(T_j)]=\hat{\rho }(S_{i_1},S_{i_2}){\sigma }_{i_1}{\sigma }_{i_2}+\mathbb{E}[S_{i_1}(T_j)]\mathbb{E}[S_{i_2}(T_j)],\end{array}$(16) ): (A1) $\begin{array}{r}a{\sigma }_1{W}_1(T)+b{\sigma }_2{W}_2(T)+c{\sigma }_3{W}_3(T)\stackrel{\mathrm{d}}{=}\hat{\sigma }\tilde{W}(T),\end{array}$(A1) with (A2) $\begin{array}{r}{\hat{\sigma }}^2=a^2{\sigma }_1^2+b^2{\sigma }_2^2+c^2{\sigma }_3^2+2{\rho }_{1,2}ab{\sigma }_1{\sigma }_2+2{\rho }_{1,3}ac{\sigma }_1{\sigma }_3+2{\rho }_{2,3}ab{\sigma }_2{\sigma }_3,\end{array}$(A2) therefore (A3) $\begin{array}{r}\sum _{k=1}^{X_{\mathcal{P}}(T)}a{J}_{k,1}+\sum _{k=1}^{X_{\mathcal{P}}(T)}b{J}_{k,2}+\sum _{k=1}^{X_{\mathcal{P}}(T)}b{J}_{k,3}=\sum _{k=1}^{X_{\mathcal{P}}(T)}(a{J}_{k,1}+b{J}_{k,2}+b{J}_{k,3})=:\sum _{k=1}^{X_{\mathcal{P}}(T)}{\hat{J}}_k,\end{array}$(A3) where $\hat{J}\equiv {\hat{J}}_k=a{J}_{k,1}+b{J}_{k,2}+c{J}_{k,3}\sim \mathcal{N}((a+b+c){\mu }_J,(a^2+b^2+c^2){\sigma }_J^2)$; thus we have: (A4) $\begin{array}{r}\mathbb{E}[{\mathrm{e}}^{a{X}_1+b{X}_2+c{X}_3}]\stackrel{\perp \perp }{=}{\mathrm{e}}^{a{\mu }_1+b{\mu }_2+c{\mu }_3}\mathbb{E}\left[{\mathrm{e}}^{\hat{\sigma }\tilde{W}(T)}\right]\mathbb{E}\left[{\mathrm{e}}^{\sum _{k=1}^{X_{\mathcal{P}}(T)}{\hat{J}}_k}\right],\end{array}$(A4) with $\hat{\sigma }\tilde{W}(T)$ defined in (A1(A1) $\begin{array}{r}a{\sigma }_1{W}_1(T)+b{\sigma }_2{W}_2(T)+c{\sigma }_3{W}_3(T)\stackrel{\mathrm{d}}{=}\hat{\sigma }\tilde{W}(T),\end{array}$(A1) ) and ${\hat{J}}_k$ defined in (A3(A3) $\begin{array}{r}\sum _{k=1}^{X_{\mathcal{P}}(T)}a{J}_{k,1}+\sum _{k=1}^{X_{\mathcal{P}}(T)}b{J}_{k,2}+\sum _{k=1}^{X_{\mathcal{P}}(T)}b{J}_{k,3}=\sum _{k=1}^{X_{\mathcal{P}}(T)}(a{J}_{k,1}+b{J}_{k,2}+b{J}_{k,3})=:\sum _{k=1}^{X_{\mathcal{P}}(T)}{\hat{J}}_k,\end{array}$(A3) ). Using the result from proof in Appendix 1 we find: (A5) $\begin{array}{rl}{\omega }_X(a,b,c)& :=\mathbb{E}[{\mathrm{e}}^{a{X}_1+b{X}_2+c{X}_3}]\\ & =\exp (a{\mu }_1+b{\mu }_2+c{\mu }_3+\frac{1}{2}\sigma T+{\xi }_{p}T\mathbb{E}[{\mathrm{e}}^{\hat{J}}-1])\\ & =\exp [a{\mu }_1+b{\mu }_2+c{\mu }_3+\frac{1}{2}\hat{\sigma }T+{\xi }_{p}T({\mathrm{e}}^{(a+b+c){\mu }_J+\frac{1}{2}(a^2+b^2+c^2){\sigma }_J^2}-1)],\end{array}$(A5) where ${\sigma }^2$ is defined in (A2(A2) $\begin{array}{r}{\hat{\sigma }}^2=a^2{\sigma }_1^2+b^2{\sigma }_2^2+c^2{\sigma }_3^2+2{\rho }_{1,2}ab{\sigma }_1{\sigma }_2+2{\rho }_{1,3}ac{\sigma }_1{\sigma }_3+2{\rho }_{2,3}ab{\sigma }_2{\sigma }_3,\end{array}$(A2) ).

Combining these results with Lemma 3.1 we find the following expression for the correlation: $\overline{\rho }(S_{i_1,i_2},S_{i_3})$, $\begin{array}{r}\overline{\rho }(S_{i_1,i_2},S_{i_3})=\frac{{\omega }_X(c_1,c_2,c_3)-{\omega }_X(c_1,c_2){\varphi }_{X_{i_3}}(-i{c}_3)}{\sqrt{({\omega }_X(2c_1,2c_2)-{\omega }_X^2(c_1,c_2))({\varphi }_{X_{i_3}}(-2i{c}_3)-{\varphi }_{X_{i_3}}^2(-i{c}_3))}},\end{array}$where: $c_1:={\hat{\alpha }}_{i_1,j,1}$, $c_2:={\hat{\alpha }}_{i_2,j,1}$, $c_3:={\hat{\alpha }}_{i_3,j,1}$ and where ${\omega }_X(c_1,c_2,c_3)$ is defined in (A5(A5) $\begin{array}{rl}{\omega }_X(a,b,c)& :=\mathbb{E}[{\mathrm{e}}^{a{X}_1+b{X}_2+c{X}_3}]\\ & =\exp (a{\mu }_1+b{\mu }_2+c{\mu }_3+\frac{1}{2}\sigma T+{\xi }_{p}T\mathbb{E}[{\mathrm{e}}^{\hat{J}}-1])\\ & =\exp [a{\mu }_1+b{\mu }_2+c{\mu }_3+\frac{1}{2}\hat{\sigma }T+{\xi }_{p}T({\mathrm{e}}^{(a+b+c){\mu }_J+\frac{1}{2}(a^2+b^2+c^2){\sigma }_J^2}-1)],\end{array}$(A5) ) and ${\omega }_X(c_1,c_2)$ in (21(21) $\begin{array}{rl}\mathbb{E}[{\mathrm{e}}^{a{X}_1(T)+b{X}_2(T)}|\mathcal{F}(t_0)]& =\exp [a{\mu }_1+b{\mu }_2+\frac{1}{2}{\hat{\sigma }}^{2}T+{\xi }_{p}T({\mathrm{e}}^{(a+b){\mu }_J+\frac{1}{2}(a^2+b^2){\sigma }_J^2}-1)]\\ & =:{\omega }_X(a,b),\end{array}$(21) ).