1 Introduction↩︎

A financial security is a financial contract whose value at maturity is determined by the price process of the underlying asset. A derivative is a financial security whose value/price is derived from one or more underlying assets. There are four types of derivatives: options, forwards, futures and swaps. A swap is a financial derivative in which two counter parties agree to trade future cash flows, with the size of the cash flow decided at the start of the deal.
Understanding the movement of a market is critical in the financial sector in order to correctly hedge and speculate on the underlying asset. Volatility and variance swaps are forward contracts on realized volatility and variance of the underlying stock respectively. In a stock market volatility and variance of stock are good indicators for many reasons such as the future fluctuation of the stock price. Investors or traders have an insight on future fluctuation of the stock price. Therefore, volatility and variance swaps provide investors with tools to hedge against or speculate on fluctuations in stock price volatility.

Most of the literature has focused on pricing volatility and variance swaps for a single asset, with several studies conducted in this area under various stochastic volatility models. The classical [1] formula assumed the volatility of a stock is constant. However, such an assumptions for financial model is not consistent for pricing financial derivatives in a stock market. To address this limitation, the stochastic volatility framework was introduced and has been further refined by financial researchers to reduce model risk and improve pricing accuracy. In [2], a new probabilistic approach using the Heston model was proposed to study volatility, variance, covariance and correlation swaps for financial markets. The impact of discrete sampling on the valuation of options on realized variance in the Heston model is investigated in [3], which establishes analytical methodology for pricing and hedging options on realized variance in the Heston model supplemented with jumps in asset returns and variance. [4] studied variance and volatility swap, and introduced the pricing strategies for some popular models (like Black-Scholes model, Merton jump diffusion model, Heston model, 3/2 model, and GARCH models). They have estimated the various parameters in the models using option-based or price-based approach, and concluded that Black-Scholes model, Merton jump diffusion, and GARCH models didn’t perform well for pricing variance and volatility swap. In some diffusion based models the volatility is driven by a Brownian motion, which can be correlated with the underlying asset; such models account the different "stylized facts". Non-Gaussian processes of Ornstein-Uhlenbeck (OU) type offer the possibility of capturing important distribution deviations from Gaussianity and for flexible modeling of dependence structure [5], [6]. The non-Gaussian OU stochastic volatility model was used by [7] to study volatility and variance swaps. [8] developed analytical solution for pricing volatility and variance swaps for Barndorff-Nielsen and Shephard (BNS) process driven financial markets. [9] modeled variance and volatility swap using the superposition of BNS type model. Analytical formulas for the arbitrage free prices for the weighted variance and weighted volatility swap is conducted by [10] under the frame work of the BNS type stochastic volatility model. All the aforementioned authors focused on pricing volatility and variance swaps for a single underlying asset. However, in today’s complex financial transaction’s investors often diversify their portfolios accross multiple assets to minimize the risk of the unobserved financial crisis. This motivates the development of multivariate variance swap, which allow contracts to be based on the joint behavior of two or more assets. [11] and [12] proposed methods for pricing generalized variance swap in financial markets with Markov-modulated volatilities and in the multi-asset setting using the Barndorff–Nielsen and Shephard (BNS) model. Their approach utilized the trace and maximum Eigenvalue of the asset return covariance matrix as measures of portfolio variability. The authors concluded that the maximum eigenvalue provides a more informative measure than the trace based swap, as it not only captures the variances but also incorporates the covariances among asset returns. In our study, we adopt the same multi-asset stochastic volatility framework, although various studies such as [13], [14] propose alternative formulations based on Lévy-driven asset dynamics.

The Generalized Variance of n-dimensional random vector variable \(X\) is defined as the determinant of its variance-covariance matrix introduced by [15]. [16] proposed the generalized variance metric to measure the composite multivariate degree of inequality (dispersion) of any multidimensional cloud. This research proposal focuses on formulating the multiasset stochastic volatility dynamics under both the Heston and Barndorff–Nielsen–Shephard (BNS) models. We define the generalized variance method based on the covariance structure of asset returns and propose a method to price the multivariate variance swap accordingly. Although the framework is applicable to any \(n\) asset portfolio, we illustrate our approach using a subset of 9 selected stocks, grouped into three portfolios of 3 assets each, to highlight how pricing varies across different asset combinations.

The rest of the paper is organized as follows. In Section 2, we provide a brief overview of the variance swap, the multivariate variance swap, and the portfolio return covariance matrix. In sections 3 and 4, we derive the multivariate variance swap for the multi-asset dynamics of Heston model and BNS model respectively. In Section 5, we provide the numerical results of the model fitting and parameter estimation for multivariate variance swap. Figures, tables, and correlation matrix are discussed under this section. Finally, a brief conclusion is provided in Section [sec6], and some additional derivations are included in the Appendix.

2 Variance swap and the portfolio return covariance matrix↩︎

The payoff at expiration \(T\) is \[N (\sigma_R^2 - K_{\text{var}}),\] where \(N\) is the notional amount, \(\sigma_R^2\) is the realized variance, and \(K_{\text{var}}\) is the variance strike price. Calculating realized variance in discrete and continuous time have the following formulas:

• In discrete time, \[\sigma^{2}_{R_{d}}=\frac{n}{T(n-1)}\sum_{i=1}^{n} \left( R_i-\bar{R}\right)^2,\] where \(R_i=\log\big(\frac{S_{i+1}}{S_i}\big)\), \(\Bar{R}\) is mean of the log returns, \(n\) is the number of log return observations, \(S_i\) is the price of the underlying asset at time \(t_i\), and \(\frac{n}{T}\) is the annualization factor if the maturity \(T\) is assumed in years.

• In continuous time, \[\sigma^{2}_{R_{c}}=\lim_{n\to\infty} \frac{n}{T(n-1)}\sum_{i=1}^{n} \left(\log\left( \frac{S_{i+1}}{S_i}\right)\right)^2 =\frac{1}{T}\int_{0}^{T}\sigma^{2}_{t} dt,\]
  where \(\sigma^{2}_{t}\) is the variance over an infinitesimal time period \(t\).

To price the multivariate variance swap, we need to construct the covariance matrix of the returns (i.e portofolio return covariance matrix). Assume that \(r_1,r_2,r_3,...,r_n\) be the dynamics of the log returns of the given stochastic volatility model. Hence, the portfolio return covariance matrix is given by \[\Sigma= \begin{pmatrix} \text{Var}\left(r_1\right) & \text{Cov}\left(r_1,r_2\right) &\cdots & \text{Cov}\left(r_1,r_n\right) \\ \text{Cov}\left(r_2,r_1\right) & \text{Var}\left(r_2\right) &\cdots& \text{Cov}\left(r_2,r_n\right) \\ \vdots & \vdots & \ddots & \vdots\\ \text{Cov}\left(r_n,r_1\right) & \text{Cov}\left(r_n,r_2\right) &\cdots & \text{Var}\left(r_n\right) \end{pmatrix}.\] The instantaneous variance of the multi-asset dynamics is given by the determinant of the covariance matrix of the log-returns dynamics. Hence, to price a multivariate variance swap using the generalized variance method, we compute the determinant \(\Sigma\) (i.e., \(|\Sigma|\)) to obtain the instantaneous variance of the portfolio return. The price \(P_{var}\) of a multivariate variance swap with strike price \(K_{\text{var}}\) is \[P_{var} = \mathop{\mathrm{\mathbb{E}}}\left[ e^{-rT} \left( \sigma_R^2 - K_{\text{var}} \right) \right],\] where \[\sigma_R^2 = \frac{1}{T} \int_0^T |\Sigma| \, dt.\]

3 Multivariate variance swap for Heston model↩︎

In this section, we introduce the multivariate Heston model and derive the corresponding expression for the variance swap using the generalized variance approach. We assume that the price process of \(n\) assets, denoted by \(S_t = (S_t^{1}, S_t^{2}, \ldots, S_t^{n}),\) follows a multivariate extension of the Heston stochastic volatility model. This framework allows us to capture both the individual dynamics of each asset and the correlation structure among their volatilities. The formal formulation of the multivariate Heston model and the corresponding multivariate variance swap is presented below.

3.1 Multivariate Heston model↩︎

Let \((\Omega, \mathcal{F}, \{\mathcal{F}_t\}, P)\) be a probability space with filtration \(\mathcal{F}_t\), \(t \in [0,T]\). [2] assumes that the underlying asset \(S_t\) operates in a risk-neutral world, with the following dynamics: \[\begin{align} \frac{dS_t}{S_t} &= r_t dt + \sigma_t dB_t, \\ d\sigma_t^2 &= k (\theta^2 - \sigma_t^2) dt + \gamma \sigma_t dW_t, \quad \sigma_0^2 > 0, \end{align} \label{eq:hes1}\tag{1}\] where \(r_t\) is a deterministic interest rate, \(\sigma_0\) and \(\theta\) are the initial and long-term volatilities, \(k > 0\) is the reversion speed, \(\gamma > 0\) is the volatility of volatility, and \(B_t\) and \(W_t\) are independent standard Wiener processes.

The variance \(\sigma_t^2\) follows a Cox-Ingersoll-Ross process, as shown in the second part of equation 1 . The change-of-time method is discussed in [17] and [2] applied the change-of-time method to solve this process, and its solution for \(\sigma_t^2\) is given by \[\sigma_t^2 = e^{-kt} \left( \sigma_0^2 - \theta^2 + \tilde{W}(\phi^{-1}_t) \right) + \theta^2. \label{eq:hes2}\tag{2}\] [2] state that \(\tilde{W}(t)\) is an \(\mathcal{F}_t\)-measurable one-dimensional Wiener process, and \(\phi^{-1}_t\) is the inverse function of \(\phi_t\), where \[\phi_t=\gamma^{-1}\int^t_0\{e^{\kappa\phi_s}(\sigma^2_0-\theta^2_0+\tilde{W}(t)+\theta^2e^{2\kappa\phi_s})\}^{-1}ds\] such that \(\tilde{W}(\phi^{-1}_t)\) is a random process with \[\mathop{\mathrm{\mathbb{E}}}[\tilde{W}(\phi^{-1}_t)] = 0.\] Thus from equation@eq:eq:hes2 we obtain \[\mathop{\mathrm{\mathbb{E}}}[\sigma_t^2] = e^{-kt} (\sigma_0^2 - \theta^2) + \theta^2. \label{eq:hes4}\tag{3}\]
For \(n\) stock price processes \(S^i_t\), where \(i = 1, 2, 3, \dots, n\), the multiple stock following Heston model is defined as \[\begin{align} \frac{dS^i_t}{S^i_t} &= \mu^i_t dt + \sigma^i_t dB^i_t, \\ d(\sigma^i)^2_t &= k_i (\theta_i^2 - (\sigma^i)^2_t) dt + \gamma_i \sigma^i_t dW^i_t, \quad (\sigma^i)^2_0 > 0, \end{align} \label{eq:hes5}\tag{4}\] where, for \(i=1,2,3,...,n\), \(\mu^i_t\) is a deterministic function, and \(k_i, \theta_i, \gamma_i > 0\) are the reversion speed, long-term volatility, and volatility of volatility, respectively. The correlation of the stock prices is \([B^l_t, B^m_t] = c_{lm}(t) \, dt,\) where \(c_{lm}(t)\) is a deterministic function of time t for \(l,m \in \{1,2,3,...,n\}\). At fixed time \(t\), \(c_{lm}(t)\) can be calculated using stock price data and be regarded as a constant \(c_{lm}\). The wiener processes \(W^i_t\) in the squared volatility process are independent of each other for \(i=1,2,3,...,n\), and \(B^i_t\) are independent of \(W^i_t\).

3.2 Multivariate variance swap using generalized variance method for Heston model↩︎

For simplicity, we begin by working with three stocks and then generalize the results to \(n\) stocks, assuming that the stocks follow the Heston model. The \(3\times3\) portfolio return covariance matrix for the Heston model equation 4 is \[\begin{align} \label{Omega951} \Sigma_1 = \begin{pmatrix} (\sigma^1)_t^2 & c_{12}\sigma^1_t\sigma^2_t & c_{13}\sigma^1_t\sigma^3_t\\ c_{21}\sigma^2_t\sigma^1_t & (\sigma^2)_t^2 & c_{23}\sigma^2_t\sigma^3_t\\ c_{31}\sigma^3_t\sigma^1_t & c_{32}\sigma^3_t\sigma^2_t & (\sigma^3)_t^2 \end{pmatrix}. \end{align}\tag{5}\] Now, we need to compute \(\sigma^2_R\), where \(\sigma^2_R\) denotes the realized variance computed from the determinant of the portfolio covariance matrix \(\Sigma_1\). It is defined as the average of the instantaneous variance in the period [0,T], which is given by \[\sigma_R^2 = \frac{1}{T}\int^T_0|\Sigma_1|dt. \label{eq:hes20}\tag{6}\]
In Appendix 7, we obtain the determinant of 5 as below \[\begin{align} |\Sigma_1| = |C| \prod_{i=1}^{3} (\sigma^i_t)^2, \label{eq:hes21} \end{align}\tag{7}\] where \(C = (c_{lm})_{1\leq l,m \leq 3}\) is the correlation matrix of stock prices calculated by the stock price data. By the independence of \((\sigma^i_t)^2\), \(i = 1, 2, 3\), we have \[\begin{align} \mathop{\mathrm{\mathbb{E}}}[| \Sigma_1 |] = |C| \prod_{i=1}^{3} \mathop{\mathrm{\mathbb{E}}}[(\sigma^i_t)^2]. \label{eq:hes22} \end{align}\tag{8}\] We notice that the 3-dimensional case of 3 is \[\begin{align} \mathop{\mathrm{\mathbb{E}}}[(\sigma^i_t)^2] = e^{-k_i t} \big((\sigma^i_0)^2 - \theta_i^2\big) + \theta_i^2, \quad i = 1, 2, 3. \end{align}\] Using the above equation, we obtain \[\begin{align} & \prod_{i=1}^{3} \mathop{\mathrm{\mathbb{E}}}[(\sigma^i_t)^2] \nonumber\\ = \; & e^{-(k_3+k_2+k_1) t} \big((\sigma^3_0)^2 - \theta_3^2\big) \big((\sigma^2_0)^2 - \theta_2^2\big) \big((\sigma^1_0)^2 - \theta_1^2\big) + e^{-(k_3+k_2) t} \big((\sigma^3_0)^2 - \theta_3^2\big) \big((\sigma^2_0)^2 - \theta_2^2\big) \theta_1^2 \notag\\ & + e^{-(k_3+k_1) t} \big((\sigma^3_0)^2 - \theta_3^2\big) \big((\sigma^1_0)^2 - \theta_1^2\big) \theta_2^2 + e^{-(k_2+k_1) t} \big((\sigma^2_0)^2 - \theta_2^2\big) \big((\sigma^1_0)^2 - \theta_1^2\big) \theta_3^2 \notag\\ & + e^{-k_3 t} \big((\sigma^3_0)^2 - \theta_3^2\big) \theta_2^2 \theta_1^2 + e^{-k_2 t} \big((\sigma^2_0)^2 - \theta_2^2\big) \theta_3^2 \theta_1^2 + e^{-k_1 t} \big((\sigma^1_0)^2 - \theta_1^2\big) \theta_3^2 \theta_2^2 + \theta_3^2 \theta_2^2 \theta_1^2. \label{eq:hes23} \end{align}\tag{9}\]

4 Multivariate swap for Barndorff-Nielsen and Shephard model↩︎

In this section, we introduce the multi-asset Barndorff-Nielsen and Shephard (BNS) model and derive its variance swap using the generalized variance method. The BNS framework is particularly well suited for capturing discontinuities in asset prices through its incorporation of jumps in the volatility process. By extending the model to a multi-asset setting, we account for both individual stock dynamics and shared jump risk across assets. The generalized variance approach, based on the determinant of the return covariance matrix, provides a tractable means of quantifying instantaneous portfolio variance under the BNS dynamics. The formal structure and derivation of the variance swap under this model are outlined below.

4.1 Multivariate BNS model↩︎

Consider a financial market with a risk-free asset yielding a constant return rate \(r\) and two stocks traded up to a fixed exercise date \(T\). Barndorff-Nielsen and Shephard [6] [5] modeled the stock price process \(S = \{S_t\}_{t \geq 0}\) on a filtered probability space \((\Omega, \mathcal{F}, \{\mathcal{F}_t\}_{0 \leq t \leq T}, P)\), carrying a standard Brownian motion \(W_t\) and an independent, positive, non-decreasing Lévy process \(Z_{\lambda t}\). The dynamics are \[\begin{align} S_t & = S_0 e^{X_t}, \nonumber\\ dX_t & = (\mu + \beta \sigma_t^2) dt + \sigma_t dW_t + \rho dZ_{\lambda t}, \tag{11}\\ d\sigma_t^2 & = -\lambda \sigma_t^2 dt + dZ_{\lambda t}, \quad \sigma_0^2 > 0, \tag{12} \end{align}\] where \(\mu, \beta, \rho, \lambda \in \mathbb{R}\), \(\lambda > 0\), and \(\rho \leq 0\) represents the leverage effect. The process \(Z_t\) is a subordinator (a Lévy process with non-decreasing paths and no Gaussian component), referred to as the background driving Lévy process (BDLP). The filtration \(\mathcal{F}_t\) is the augmentation of the filtration generated by \((W, Z)\).

Non-Gaussian Ornstein-Uhlenbeck processes, driven by subordinators, can model properties such as heavy-tailed log-returns, aggregational Gaussianity, and volatility clustering [6]. Following [18], the BDLP \(Z\) satisfies

• \(Z\) has no deterministic drift, and its Lévy measure has density \(w(x)\). The cumulant transform is: \[\kappa(\theta) = \log \mathop{\mathrm{\mathbb{E}}}[e^{\theta Z_1}] = \int_{\mathbb{R}_+} (e^{\theta x} - 1) w(x) \, dx,\] where it exists.

• Let \(\hat{\theta} = \sup\{\theta \in \mathbb{R} : \kappa(\theta) < +\infty\}\), then \(\hat{\theta} > 0\).

• \(\lim_{\theta \to \hat{\theta}} \kappa(\theta) = +\infty\).

Under an equivalent martingale measure (EMM), as in [18] and [19], the dynamics 11 and 12 become: \[\begin{align} dX_t & = b_t dt + \sigma_t dW_t + \rho dZ_{\lambda t}, \tag{13}\\ d\sigma_t^2 & = -\lambda \sigma_t^2 dt + dZ_{\lambda t}, \quad \sigma_0^2 > 0, \tag{14} \end{align}\] where \[b_t = r - \lambda \kappa(\rho) - \frac{1}{2} \sigma_t^2,\] and \(W_t\) and \(Z_t\) are a Brownian motion and Lévy process, respectively, under the EMM. The solution to 14 is: \[\sigma_t^2 = e^{-\lambda t} \sigma_0^2 + \int_0^t e^{-\lambda (t-s)} dZ_{\lambda s}, \label{eq:hes13}\tag{15}\] where \(\sigma_t^2\) is strictly positive and bounded below by \(e^{-\lambda t} \sigma_0^2\). The instantaneous variance of the log-return from 13 is \((\sigma_t^2 + \rho^2 \lambda \text{Var}[Z_1]) dt\).

For \(n\) stocks with log-returns \(dX^i_t\) for \(i\in\{1,2,3,...,n\}\), the multi-stock follows BNS model is given by \[\begin{align} dX^i_t & = b^i_t dt + (\sigma^i)_t dB^i_t + \rho_i dZ^*_{\lambda t} \tag{16}\\ d(\sigma^i)^2_t & = -\lambda (\sigma^i)^2_t dt + dZ^i_{\lambda t}, \quad (\sigma^i)^2_0> 0 \tag{17} \end{align}\] where \(B^i_t\) are Wiener processes with \([B^l_t, B^m_t] = c_{lm}(t) \, dt,\) where \(c_{lm}(t)\) is a deterministic function of time t for \(l,m \in \{1,2,3,...,n\}\), and \(Z^i_t\) are independent Lévy processes. At fixed time \(t\), \(c_{lm}(t)\) are the same values as in Heston model.

4.2 Multivariate variance swap using generalized variance method for BNS model↩︎

We consider three stocks assuming that the stocks follow the BNS model. The \(3 \times 3\) portfolio return covariance matrix for the BNS model equation 16 is \[\begin{align} \label{Omega952} \Sigma_2 = \begin{pmatrix} (\sigma^1)_t^2 + \rho_1^2 \lambda \text{Var}[Z_1^*] & c_{12} \sigma^1_t\sigma^2_t + \rho_1 \rho_2 \lambda\text{Var}[Z_1^*] & c_{13} \sigma^1_t\sigma^3_t + \rho_1 \rho_3 \lambda\text{Var}[Z_1^*]\\ c_{21} \sigma^2_t\sigma^1_t + \rho_2 \rho_1 \lambda\text{Var}[Z_1^*] & (\sigma^2)_t^2 + \rho_2^2 \lambda \text{Var}[Z_1^*] & c_{23} \sigma^2_t\sigma^3_t + \rho_2 \rho_3 \lambda\text{Var}[Z_1^*]\\ c_{31} \sigma^3_t\sigma^1_t + \rho_3 \rho_1 \lambda\text{Var}[Z_1^*] & c_{32} \sigma^3_t\sigma^2_t + \rho_3 \rho_2 \lambda\text{Var}[Z_1^*] & (\sigma^3)_t^2 + \rho_3^2 \lambda \text{Var}[Z_1^*] \end{pmatrix}. \end{align}\tag{18}\] Similar to equation 6 , we also need to compute the realized variance \(\sigma^2_R\) defined as the average of the instantaneous variance in the period \([0,T]\), which is given by \[\begin{align} \label{sigma94295R} \sigma_R^2 = \frac{1}{T}\int^T_0|\Sigma_2|dt. \end{align}\tag{19}\] In Appendix 7, we obtain the determinant of 18 as below \[\begin{align} \label{Sigma95295det} & | \Sigma_2 | = \nonumber\\ & |C| \prod_{i=1}^{3}(\sigma^i_t)^2 + \lambda \text{Var}[Z_1^*] |C| \Big( \delta_{11}\rho_1^2(\sigma^3_t)^2 (\sigma^2_t)^2 + \delta_{22}\rho_2^2(\sigma^3_t)^2 (\sigma^1_t)^2 + \delta_{33}\rho_3^2 (\sigma^2_t)^2 (\sigma^1_t)^2 \nonumber\\ & \qquad \qquad \qquad \qquad \qquad \qquad + 2 \delta_{21}\rho_2\rho_1(\sigma^3_t)^2 \sigma^2_t \sigma^1_t + 2 \delta_{31}\rho_3\rho_1 \sigma^3_t (\sigma^2_t)^2 \sigma^1_t + 2 \delta_{32}\rho_3\rho_2 \sigma^3_t \sigma^2_t (\sigma^1_t)^2 \Big), \end{align}\tag{20}\] where \((\delta_{ij})_{1\leq i,j \leq 3} = C^{-1}\) can be also calculated by the stock price data. Based on our construction in equation 17 where the variance processes are driven by independent processes, we have \[\begin{align} \label{E95Sigma952} & \mathop{\mathrm{\mathbb{E}}}[| \Sigma_2 |] = \nonumber\\ & |C| \prod_{i=1}^{3}\mathop{\mathrm{\mathbb{E}}}[(\sigma^i_t)^2] + \lambda \text{Var}[Z_1^*] |C| \Big( \delta_{11}\rho_1^2 \mathop{\mathrm{\mathbb{E}}}[(\sigma^3_t)^2] \mathop{\mathrm{\mathbb{E}}}[(\sigma^2_t)^2] + \delta_{22}\rho_2^2 \mathop{\mathrm{\mathbb{E}}}[(\sigma^3_t)^2] \mathop{\mathrm{\mathbb{E}}}[(\sigma^1_t)^2] \nonumber\\ & \qquad \qquad \qquad \qquad \qquad \qquad \quad + \delta_{33}\rho_3^2 \mathop{\mathrm{\mathbb{E}}}[(\sigma^2_t)^2] \mathop{\mathrm{\mathbb{E}}}[(\sigma^1_t)^2] + 2 \delta_{21}\rho_2\rho_1 \mathop{\mathrm{\mathbb{E}}}[(\sigma^3_t)^2] \mathop{\mathrm{\mathbb{E}}}[\sigma^2_t] \mathop{\mathrm{\mathbb{E}}}[\sigma^1_t] \nonumber\\ & \qquad \qquad \qquad \qquad \qquad \qquad \quad + 2 \delta_{31}\rho_3\rho_1 \mathop{\mathrm{\mathbb{E}}}[\sigma^3_t] \mathop{\mathrm{\mathbb{E}}}[(\sigma^2_t)^2] \mathop{\mathrm{\mathbb{E}}}[\sigma^1_t] + 2 \delta_{32}\rho_3\rho_2 \mathop{\mathrm{\mathbb{E}}}[\sigma^3_t] \mathop{\mathrm{\mathbb{E}}}[\sigma^2_t] \mathop{\mathrm{\mathbb{E}}}[(\sigma^1_t)^2] \Big). \end{align}\tag{21}\] From equation 15 , we compute the expectation of the variance process: \[\begin{align} \label{E} \mathop{\mathrm{\mathbb{E}}}[\sigma_t^2] & = e^{-\lambda t} \mathop{\mathrm{\mathbb{E}}}[\sigma_0^2] + \int_0^t e^{-\lambda (t-s)} \mathop{\mathrm{\mathbb{E}}}[dZ_{\lambda s}] \nonumber\\ & = e^{-\lambda t} \sigma_0^2 + \int_0^t e^{-\lambda (t-s)} \;\kappa_1 \lambda ds \nonumber\\ & = e^{-\lambda t}(\sigma_0^2 - \kappa_1) + \kappa_1, \end{align}\tag{22}\] where \(\kappa_1\) is the 1st cumulant of the subordinator driving the variance process, i.e. the mean of jump size, which can be calculated by the stock price data. In the 3-dimensional case, we have \[\begin{align} \mathop{\mathrm{\mathbb{E}}}[(\sigma_t^i)^2] = e^{-\lambda t}\big((\sigma_0^i)^2 - \kappa_1^i\big) + \kappa_1^i, \quad i = 1,2,3. \end{align}\] Using the above equation, we obtain the followings: \[\begin{align} \label{E0} & \prod_{i=1}^{3} \mathbb{E}[(\sigma_t^{i})^2] \nonumber\\ =\;& e^{-3\lambda t} \big((\sigma_0^{1})^2 - \kappa_1^{1}\big)\big((\sigma_0^{2})^2 - \kappa_1^{2}\big)\big((\sigma_0^{3})^2 - \kappa_1^{3}\big) \nonumber\\ & + e^{-2\lambda t} \Big[ \kappa_1^{1} \big((\sigma_0^{2})^2 - \kappa_1^{2}\big)\big((\sigma_0^{3})^2 - \kappa_1^{3}\big) + \kappa_1^{2} \big((\sigma_0^{1})^2 - \kappa_1^{1}\big)\big((\sigma_0^{3})^2 - \kappa_1^{3}\big) \nonumber\\ &\qquad \qquad + \kappa_1^{3} \big((\sigma_0^{1})^2 - \kappa_1^{1}\big)\big((\sigma_0^{2})^2 - \kappa_1^{2}\big) \Big] \nonumber\\ & + e^{-\lambda t} \Big[ \kappa_1^{1} \kappa_1^{2} \big((\sigma_0^{3})^2 - \kappa_1^{3}\big) + \kappa_1^{1} \kappa_1^{3} \big((\sigma_0^{2})^2 - \kappa_1^{2}\big) + \kappa_1^{2} \kappa_1^{3} \big((\sigma_0^{1})^2 - \kappa_1^{1}\big) \Big] + \kappa_1^{1} \kappa_1^{2} \kappa_1^{3}, \end{align}\tag{23}\]

\[\begin{align} \label{E1} & \mathbb{E}\big[(\sigma_t^{3})^2\big] \mathbb{E}\big[(\sigma_t^{2})^2\big] \nonumber\\ =\;& e^{-2\lambda t} \big((\sigma_0^{3})^2 - \kappa_1^{3}\big)\big((\sigma_0^{2})^2 - \kappa_1^{2}\big) + e^{-\lambda t} \Big[ \kappa_1^{3}\big((\sigma_0^{2})^2 - \kappa_1^{2}\big) + \kappa_1^{2}\big((\sigma_0^{3})^2 - \kappa_1^{3}\big) \Big] + \kappa_1^{3} \kappa_1^{2}, \end{align}\tag{24}\]

\[\begin{align} \label{E2} & \mathbb{E} \big[(\sigma_t^{3})^2\big] \mathbb{E}\big[(\sigma_t^{1})^2\big] \nonumber\\ =\;& e^{-2\lambda t} \big((\sigma_0^{3})^2 - \kappa_1^{3}\big)\big((\sigma_0^{1})^2 - \kappa_1^{1}\big) + e^{-\lambda t} \Big[ \kappa_1^{3}\big((\sigma_0^{1})^2 - \kappa_1^{1}\big) + \kappa_1^{1}\big((\sigma_0^{3})^2 - \kappa_1^{3}\big) \Big] + \kappa_1^{3} \kappa_1^{1}, \end{align}\tag{25}\]

\[\begin{align} \label{E3} & \mathbb{E}\big[(\sigma_t^{2})^2\big] \mathbb{E}\big[(\sigma_t^{1})^2\big] \nonumber\\ =\;& e^{-2\lambda t} \big((\sigma_0^{2})^2 - \kappa_1^{2}\big)\big((\sigma_0^{1})^2 - \kappa_1^{1}\big) + e^{-\lambda t} \Big[ \kappa_1^{2}\big((\sigma_0^{1})^2 - \kappa_1^{1}\big) + \kappa_1^{1}\big((\sigma_0^{2})^2 - \kappa_1^{2}\big) \Big] + \kappa_1^{2} \kappa_1^{1}. \end{align}\tag{26}\]

5 Model Fitting and Parameter Estimation↩︎

This section presents numerical results, model calibration, and parameter testing procedures. For the purpose of model calibration, nine stocks were selected using the “quantmod" package [21] in R [22], covering the period from January 1, 2021, to January 1, 2024. The mean, variance, and kurtosis of the 9 assets can be also found in appendix 8. These stocks were then randomly reshuffled into three groups, each consisting of three stocks, to evaluate the robustness and precision of the proposed model in different combinations of assets. The log-returns of the daily closing prices were computed for each stock. Subsequently, the covariance matrix of the log-returns was estimated at 10-day time intervals, and the determinant of each covariance matrix was calculated to capture the joint variability and dependency structure among the selected assets over time. The correlation matrix of Coca-Cola, Apple, and Tesla, the correlation matrix for Google, Microsoft, and Meta, and the correlation matrix for J.P. Morgan, Nvidia, and Amazon are in Figure 2.

Figure 1 presents histograms of the log-return distributions of the three stocks grouped together for analysis. The histograms reveal that all underlying assets display similar distributional patterns, suggesting that their return behaviors share comparable statistical characteristics. This similarity in log-return distributions implies that the assets may respond to market factors in a consistent manner, which provides a reasonable basis for conducting a joint variance swap analysis.

Figure 2 illustrates the correlation matrix among the three underlying assets considered in this analysis. The correlation heatmap indicates that the assets are generally moderately to weakly correlated with one another, except for META and Google, which exhibit a relatively strong positive correlation with a coefficient of 0.72. This suggests that META and Google tend to move in similar directions more frequently compared to the other asset pairs.

Figure 3 illustrates the cumulative returns of the closing prices for the nine selected assets, grouped as follows: the left panel presents Coca-Cola, Apple, and Tesla; the middle panel displays Google, Microsoft, and Meta; and the right panel includes JPMorgan, Nvidia, and Amazon. The figure clearly shows that the cumulative return distributions of these assets vary substantially. Some assets, such as Tesla and Nvidia, exhibit high volatility, whereas others, like J.P. Morgan, demonstrate relatively stable performance. The rationale for selecting this diverse set of covariates is to evaluate the robustness and generalization capability of the model’s predictive accuracy across assets with different volatility and market behaviors.

Figure 4 presents the realized variance of the covariance matrix along with the fitted results obtained from the Heston and BNS models. A nonlinear least squares (NLS) estimator was employed in R for parameter calibration and model validation. The estimation results indicate that all parameters in both the Heston and BNS models were statistically significant, with \(p\)-values on the order of \(2.6 \times10^{-16}\), confirming the overall robustness and suitability of the fitted models.

Based on these findings, both the Heston and BNS models demonstrated strong performance in capturing and predicting the determinant of the realized variance. However, the BNS model outperformed the Heston model in several key aspect it achieved higher predictive accuracy, reduced the mean squared error, and exhibited faster convergence. Quantitatively, the BNS model was able to explain approximately 40% more of the variation in the data that the Heston model fails to capture. Hence, on average, the BNS model provided about a 40% improvement in performance over the standard Heston model.

Multiple error metrics are employed to evaluate the performance and accuracy of our models, including the Absolute Percentage Error (APE), Average Absolute Error (AAE), Average Relative Percentage Error (ARPE), and Root Mean Square Error (RMSE). These metrics provide a comprehensive assessment of model goodness-of-fit by capturing different aspects of prediction error both in magnitude and variability. The formal definitions of these error measures, as summarized in Tables 1 and 2, are expressed as follows:

\[\begin{align} \text{APE} &= \frac{1}{\bar{\sigma}_R^2}\sum_{i=1}^n \frac{|\sigma_R^2(t_i)-\hat{\sigma}_R^2(t_i)|}{n}, \end{align}\] \[\begin{align} \text{AAE} &= \frac{1}{n} \sum_{i=1}^n |\sigma_R^2(t_i)-\hat{\sigma}_R^2(t_i)|, \end{align}\] \[\begin{align} \text{ARPE} &= \frac{1}{n} \sum_{i=1}^n \frac{|\sigma_R^2(t_i)-\hat{\sigma}_R^2(t_i)|}{\sigma_R^2(t_i)} , \end{align}\] \[\begin{align} \text{RMSE} &= \sqrt{\frac{1}{n} \sum_{i=1}^n \Big(\sigma_R^2(t_i)-\hat{\sigma}_R^2(t_i)\Big)^2}. \end{align}\]

Tables 1 and 2 present the numerical error calculations for the realized variance of the multi-asset determinant covariate matrices. The reported error metrics provide a quantitative comparison between the models, highlighting their relative predictive performance. As shown in these tables, the BNS model consistently demonstrates lower error values across all measures, indicating a superior fit and improved accuracy compared to the standard Heston model. This result suggests that the BNS framework captures the underlying market dynamics more effectively, particularly in modeling volatility and cross-asset relationships.

6 Conclusion↩︎

In this paper, we introduced a generalized variance method to compute the instantaneous variance of the log-returns for multi-stock dynamics under both the Heston and Barndorff-Nielsen–Shephard (BNS) models. This method is based on the computation of the determinant of the portfolio return covariance matrix, providing a more tractable approach for estimating the realized variance of selected stocks in our numerical analysis.

We selected nine stocks and randomly reshuffled them into three groups, each containing three distinct stocks. The fitted models from both the Heston and BNS frameworks showed a strong agreement with the empirical realized variance over the selected time horizons. Notably, the BNS model demonstrated an superior fit, capturing approximately 40% more of the variance in the data that the Heston model failed to explain. Meanwhile, the BNS model demonstrated lower errors indicating an improved accuracy compared to the Heston model. This aligns with our expectations, as the BNS model is capable of capturing jumps more effectively in the underlying asset dynamics.

7 Appendix: Calculations of \(|\Sigma_1|\) and \(|\Sigma_2|\)↩︎

If the 3-dimensional covariance matrix for Heston model is given by 5 , it is trivial that \((\sigma^i)_t^2 = \sigma^i_t \sigma^i_t\), \(i = 1, 2, 3\). Then \[\begin{align} \label{DCD} \Sigma_1 = DCD, \end{align}\tag{31}\] where \(C = (c_{lm})_{1\leq l,m \leq 3}\) is the correlation matrix of stock prices which can be calculated using stock price data, and \(D = \text{diag}(\sigma^1_t, \sigma^2_t, \sigma^3_t)\). Hence, \(|\Sigma_1| = |C||D|^2\) gives equation 7 .

If the 3-dimensional covariance matrix for BNS model is given by 18 , then \[\begin{align} \Sigma_2 = \Sigma_1 + \lambda \text{Var}[Z_1^*] \rho \rho^\top, \end{align}\] where \(\rho=(\rho_1,\rho_2,\rho_3)^\top\). By the matrix determinant lemma, \[\begin{align} \label{Sigma952} | \Sigma_2 | = | \Sigma_1 + (\lambda \text{Var}[Z_1^*] \rho) \rho^\top | = | \Sigma_1 |\;(1 + \lambda \text{Var}[Z_1^*] \rho^\top \Sigma_1^{-1} \rho). \end{align}\tag{32}\] We notice that, at fixed time \(t\), the correlation matrix \(C(t)\) of stock prices can be calculated using stock price data and be regarded as a constant matrix \(C\). Simultaneously, \(C^{-1}\) can be also calculated using the same stock price data. Therefore, we denote \(C^{-1} = (\delta_{ij})_{1\leq i,j \leq 3}\). By equation 31 , we have \[\begin{align} \label{Sigma-1} \Sigma_1^{-1} = D^{-1}C^{-1}D^{-1} = \begin{pmatrix} \frac{\delta_{11}}{(\sigma^1_t)^2} & \frac{\delta_{12}}{\sigma^1_t\sigma^2_t} & \frac{\delta_{13}}{\sigma^1_t\sigma^3_t}\\ \frac{\delta_{21}}{\sigma^2_t\sigma^1_t} & \frac{\delta_{22}}{(\sigma^2_t)^2} & \frac{\delta_{23}}{\sigma^2_t\sigma^3_t}\\ \frac{\delta_{31}}{\sigma^3_t\sigma^1_t} & \frac{\delta_{32}}{\sigma^3_t\sigma^2_t} & \frac{\delta_{33}}{(\sigma^3_t)^2} \end{pmatrix}. \end{align}\tag{33}\] Finally, we obtain \(| \Sigma_2 |\) in equation 20 using 7 , 32 , and 33 .

8 Appendix: Summary Statistics↩︎

References↩︎