Co-movements, option pricing and risk management: an application to WTI versus Brent spread options

Co-moments of asset returns play a major role in financial contagion during crises. We study the properties of a particular specification of the generalized bivariate normal distribution which allows for co-volatility and co-skewness. With this probability distribution, formulae for single-name and exchange options can be evaluated quickly since they are based on one-dimensional integrals. We provide a very precise approximation formula for spread option prices and derive the corresponding greeks. We perform a day-to-day re-estimation of the probability distribution on a dataset of WTI vs Brent spread options, showing the ability of this specification to capture the salient empirical features observed in the market. Finally, we show the impact of co-movements on portfolio risk management.


Introduction
High-order co-moments play an important role in the literature on financial contagion, indeed, they are used to identify channels of contagion among asset returns. 1 Fry et al. (2010) extend correlation-based tests of contagion [see, for instance, Forbes and Rigobon (2002)] and develop statistical tests based on co-moments using a generalized normal distribution. These tests are generalized in  and . Fry-McKibbin et al. (2014) perform an extensive set of numerical experiments to show how option prices and hedging strategies are affected by co-moments. While the papers above focus on tests based on single-channel contagion, more recently Hsiao and Morley (2022) formulate tests based on multiple channels, which establish contagion by identifying changes in co-moments. The authors apply their newly-developed statistical tests to four events in which financial crises originated in one country have spread through the world. 1 See the survey of Dungey et al. (2005) for a review of methodologies on contagion.
The literature offers other approaches that allow to capture financial contagion in the context of option pricing, one of which consists of the use of mutually self-exciting jump processes (Aït-Sahalia et al., 2014;Aït-Sahalia & Hurd, 2015;Kokholm, 2016). In the context of single-name options, Melick and Thomas (1997) offer a method and a new model based on mixture of lognormals to estimate the option-based risk neutral density.
Co-moments have also a key role in the pricing of spread options. This is well documented in Fry-McKibbin et al. (2014), where the authors show the effects of high-order co-moments on spread option prices. The literature that studies spread options offers a variety of models and numerical techniques (see Carmona and Durrleman 2003, for a review). Margrabe (1978) provides the formula to price exchange options, a special case in which the exercise price is null. Apart from this case, no closed formula exists, and one must resort to numerical techniques. In this context, two approaches stand out. The first approach is the two-dimensional Fourier transform method of Hurd and Zhou (2010), which is exact, but requires a bivariate numerical integration. The second is the one-dimensional approximation of Caldana and Fusai (2013), which generalizes the lower bound of Bjerksund and Stensland (2014) (improving the approximation formula of Kirk, 1995), valid in the two-dimensional lognormal setup, for any model for which the joint characteristic function is available in closed form. This method requires only a one-dimensional Fourier inversion and provides a very accurate approximation of the true price. From the modeling side, the literature offers a variety of approaches. Huang and Kou (2006) and Cheang and Chiarella (2011) are based on jump diffusion processes. Subordinated Levy processes are used in Ballotta and Bonfiglioli (2016), while Dempster and Hong (2002) and Schneider and Tavin (2018) use multi-factor models based on stochastic volatility.
In this paper, we study a particular specification of the generalized bivariate normal distribution used by Fry-McKibbin et al. (2014) to analyze the impact of high-order co-moments on the risk associated with options written on contagious underlyings. The specification allows for co-volatility and co-skweness, and turns out to be particularly appropriate for the representation of non-linear dependence among asset returns by means of high-order co-moments. Under the specification considered in this paper, closed-form expressions of the conditional densities and the quasi-closed expressions for the marginal densities can be derived in a simple way. The relevant co-moments can also be computed by means of a one-dimensional numerical integral. This specification admits closed form formulae for single-name and exchange option in terms of one-dimensional integrals, which can be evaluated quickly. We also provide an approximated formula for spread option prices. The formula gives a lower bound for the true price. Nevertheless, extensive numerical experiments show that the bounds provided are very precise.
We use WTI versus Brent spread option prices to perform a day-to-day re-estimation of the joint (risk-neutral) probability distribution. Our dataset consists of European Call spread options with several maturities and strike prices, which allow the owner to buy WTI futures in exchange of Brent futures plus some extra cash, namely the strike price. The period under investigation covers a little more than a year of option prices including the months in which the effects of the worldwide pandemic named Covid-19 have been dramatic for the energy sector. During those months, the futures on WTI have experienced a tremendous crash, due to the lack of oil demand and the limited storage facilities in the United States. In the same period, the Brent experienced a less pronounced but still strong negative impact. The use of this dataset is particularly valuable for our purposes for at least two reasons. First, an estimation exercise on a dataset which covers a period of crisis is useful in that it allows to test the quality of fit of the considered probability distribution in periods in which classic distributions usually fail. In other words, this dataset gives us the opportunity to test our specification during a period of stressed market conditions. Second, the dataset allows us to uncover the transmission channels of contagion in terms of risk-neutral high-order comoments during that period. Then, we use the parameters obtained to highlight the contagion effects, in terms of high-order co-moments, that occurred in the last energy crisis due to the worldwide pandemic.
The remainder of the paper is organized as follows: in Sect. 2, we provide a specification of the generalized bivariate normal distribution which allows quasi-closed formulae for option prices. In Sect. 3, we apply the model to option pricing. The relevant formulae for the Greeks are also provided in Appendix B. In Sect. 4, we calibrate the model to WTI versus Brent spread option data. Section 5 assesses the impact of co-moments in risk management. Finally, in Sect. 6, we conclude.

A tractable bivariate Generalized Normal model
We consider a world with two real-valued state variables denoted by Y 1 , Y 2 . Later in this paper we will interpret the state variables as the fundamental uncertainty driving returns of two asset in the context of option pricing and risk management. The standard assumption in option pricing is the normality of the fundamental uncertainty. The huge option pricing literature extends the normality assumption in many directions. In this paper we focus on the use of a particular specification of the generalized normal distribution considered by Fry et al. (2010) and Fry-McKibbin et al. (2014), who extend the work of Lye and Martin (1993) to the multidimensional case. In particular, we specify the density of the bivariate random variable (Y 1 , Y 2 ) as a special case of the generalized normal distribution proposed by Fry-McKibbin et al. (2014): θ 2 y 2 2 + θ 3 y 1 y 2 + θ 4 y 1 y 2 2 + θ 5 y 2 1 y 2 − θ 6 y 2 1 y 2 2 − η , (1) where θ 1 = θ 2 = 1 1−ρ 2 , θ 3 = ρ 1−ρ 2 , ρ ∈ (−1, 1). We restrict our attention to the case θ 6 > 0 to guarantee that the two random variables admit finite moments of all orders. Once the parameters θ 1 , . . . , θ 6 are specified, η is determined by the relation From (1) we recover the bivariate normal density by setting θ 4 = θ 5 = θ 6 = 0, where the parameters θ 1 , θ 2 control the variance of Y 1 , Y 2 , respectively and θ 3 drives the linear dependence between the states variables. The role of the remaining parameters is to allow for a more sophisticated dependence structure. In particular, θ 4 and θ 5 drive the co-skewness, that is the dependence of the bivariate random variables (Y 1 , Y 2 2 ) and (Y 2 1 , Y 2 ). In addition, θ 6 models the interaction of the couple (Y 2 1 , Y 2 2 ), known as co-volatility. Figure 1 shows the effects of different values of these parameters. The parameters driving co-volatility distort symmetrically the density function, assigning more weight to events where one of the two components of the random variable is close to zero. The parameters driving co-skewness modify asymmetrically the previous effect, according to their sign and magnitude.
The probability density function in (1) allows for greater flexibility with respect to the normal distribution, while keeping analytic tractability. In fact, straightforward manipulation  (1) for different parametrization of the triple (θ 4 , θ 5 , θ 6 ). a (0, 0, 0); b (0, 0, 2); c (1.5, 0, 2); d (−1.5, 0, 2); e (1.5, 1.5, 2); f (0, 1.5, 2); In all panels we set ρ = 0 of the joint density function allows to derive conditional densities in closed form. Additionally, marginal densities can also be expressed in closed-form up to the solution of a one dimensional integral (which can be computed instantaneously) for the normalizing constant. These facts, which are summarized in the next proposition (see the proof in 1), turn out to be crucial for fast computation of option prices and model calibration. In terms of notation, here and throughout the rest of the paper φ(·; μ, σ 2 ) denotes the density function of a normal random variable with expected value μ and variance σ 2 whereas φ(·) denotes the standard normal density function, i.e. φ(·; 0, 1). Finally Φ(·) denotes the cumulative distribution function of a standard normal random variable.

The marginal density of Y
5. The parameter η can be computed by imposing that the integral over R of one of the two marginal densities above is unity, thus solving a one-dimensional integral.
Proposition 1 not only provides a representation of the density function in terms of a onedimensional integral, but also suggests an algorithm for the simulation from such bivariate density.
1. Use the inverse CDF method to simulate a random variate, y * 2 , from the marginal density of Y 2 , f Y 2 (·). 2. Use y * 2 and the conditional density of Y 1 |Y 2 and obtain y * 1 by simulating from the normal distribution with mean μ 1|2 (y * 2 ) and variance σ 2 1|2 (y * 2 ). Note that the inverse CDF needed in step 1. can be efficiently approximated by evaluating the non-normalized density at a properly selected grid of points (see Tsay, 2010, p.623). Therefore, the first step of the above algorithm only requires evaluating the kernel of Y 2 at a grid of points without having to compute the normalization constant η.

Co-moments
In this section, we give formulae for the relevant co-moments (co-skewness and co-volatility). A first way to derive the formula for co-skewness (identified by looking at changes in the interaction between expected returns of the first asset and volatility of the second asset) is the following: Alternatively, the same measure of co-skewness can be obtained via Replicating the arguments above, it is straightforward to derive the formula for coskewness identified by looking at changes in the interaction between expected returns of the second asset and volatility of the first asset: Finally, co-volatility can be calculated as Note that the formula for co-moments can also be derived using the results derived by Fry et al. (2010). To recall their result, let us write the log of the bivariate density as e h(y 1 ,y 2 ) dy 1 dy 2 and let us denote by θ the vector of parameters. Then Our representation of the relevant co-moments is not as elegant as that provided in (2). Nevertheless, it provides a practical and quick way to compute the relevant quantities. The representation of Fry et al. (2010) is indeed in implicit form, since the expression for η is not available.

Generalized normal distribution and option prices
In this section, we propose a lower bound for European spread options which pays at the maturity τ the amount c τ = (S 1,τ − S 2,τ − K ) + , where S 1,τ , S 2,τ are the value at maturity of the two underlying assets, K is the strike price and x + denotes the positive part of x. In the special case K = 0, the lower bound turns out to be exact, thus generalizing Margrabe (1978)'s formula for European exchange options. In the case K > 0, the lower bound is very close to the real price, thus providing a fast and reliable way to compute spread option prices under our specification for the generalized normal distribution. We specify the risk-neutral probability distribution of the underlying assets along the lines of Fry-McKibbin et al. (2014). We consider a financial market with two risky assets and denote by r the instantaneous risk-free rate. Our standing assumption is about the distribution of the two underlying assets at maturity τ , which we specify directly under the risk-neutral probability measure: (3), the distribution of the random vector (Y 1 , Y 2 ) is described by (1), the parameters σ i is the volatility of asset i, r = r if the underlying assets are stocks 0 if the underlying assets are futures is the moment generating function of Y i . The value of a European spread option 2 written on the two assets defined in (3) is given by Apart from the special case K = 0, no closed form solutions exist for such contracts. However, the literature offers a fast way to compute spread option prices by means of convenient approximations. In this context, Kirk (1995) proposes a first approximation, subsequently refined by Bjerksund and Stensland (2014). Caldana and Fusai (2013) generalize Bjerksund and Stensland (2014) to any model where the joint characteristic function can be computed explicitly. The approximation proposed in Bjerksund and Stensland (2014) produces option prices which are very close to the true price. It consists in substituting the payoff at maturity with the following approximated payoff: where With our distributional assumptions, this leads to option prices that are very rough (sometimes too rough) approximations of the true price. However, a slight modification of the approximation above, which consists in setting c = S 2,t er (τ −t) instead of c = E t S b 2,τ , proves to be effective. In the next proposition, we compute the approximated option prices.
Proposition 2 Under the assumptions of Proposition 1, the quantitỹ Single-name option prices can be obtained quickly via one-dimensional integrals. For instance, the value of a European call option written on asset S i is denotes the time-t conditional expectation under the risk neutral probability measure. Simple calculations show that the above option price can be expressed as reduces tõ We observe that (5) can be considered as an approximation for spread option prices only for the case K ≥ 0. The complementary case can be easily dealt with, because the call option can be thought of a put on the opposite price spread, i.e. S 2,τ − S 1,τ , and put values are obtained using the put-call parity (see the discussion at the beginning of Sect. 4).
Hedging the risk associated with trading in financial options involves the use of the socalled Greeks. These are the derivatives (in the mathematical sense) of the option price with respect to some quantities of interest. Two widely used hedging strategies for spread options are delta hedging and delta-gamma hedging (see, for instance, Venkatramanan and Alexander, 2011). The former involves taking positions on the underlying stocks, proportionally to the derivative of the option price with respect to each underlying, the so-called deltas of the option. The latter implies having in the hedging portfolio not only some shares of the stocks, but also some options on the two stocks, in a way proportional to the gammas of the option (the second derivatives of the option prices with respect to the underlying stocks). Both hedging strategies require the frequent rebalancing of the hedging portfolio, calling for a fast way to compute the necessary hedge ratios. In Appendix B, we provide the deltas and the gammas for the option price formula (5).
The formula in (5), which is exact in the special case of European exchange options, reduces the two-dimensional option pricing problem to the solution of one-dimensional integral, which can be easily and quickly computed with commonly available routines. This makes fast calibration of spread options possible. Although we are not the first who are able to calibrate such options (see, for instance, Schneider and Tavin, 2018) the approximation provided in (5) allows, for the first time, the use of high order co-moments in the calibration. However, before proceeding with the applications of the models, we perform a series of numerical experiments to evaluate the quality of the approximation provided by (5).

Quality of the approximation
To show the quality of the approximation detailed in Proposition 2, we run a set of numerical evaluations by comparing the approximated option prices given in (5) with the true prices. We Table 1 For column Θ i , each element of the table is defined as max (ρ,θ 4 wherec t is the lower bound derived in proposition 2, and c t is the true price obtained by means of a numerical routine for two-dimensional integration The errors are expressed in percentage refer to the true option price as the one computed by a routine for two-dimensional numerical integration, using the original representation for the density given in formula (1). We consider WTI versus Brent hypothetical spread options allowing the owner to exchange a futures on the WTI with a futures written on Brent at a given maturity. We set S 1,t = 51.26 and S 2,t = 55.4. These are the New York Mercantile Exchange (NYMEX) closing quotes, as observed on March 20, 2020, of the futures written on WTI and Brent, respectively with delivery on September 1, 2020. At the same date, we observe several spread options written on the two underlying futures, and several strike prices, ranging from K = 0 to K = 7. We use these parameters, in conjunction with a risk-free interest rate r = 0.007 and σ 1 = 0.25, σ 2 = 0.2. To isolate the effects of parameters governing co-movements on the approximation quality, we consider 6 different sets of parameters' values. In particular, we fix the set Λ = {−0.9, −0.5, −0.1, 0.1, 0.5, 0.9} and define: The values in the sets Θ 1 -Θ 6 are chosen to span the relevant parameter space as much as possible. In practice, we selected ρ between − 0.9 and 0.9, θ 4 , θ 5 between -2.8 and 2.8 and θ 6 between 0.25 and 8. The choice of this ranges for the parameters related to the co-movements has been made based on the estimation results presented in the next section. In fact, for our day-to-day re-estimation we fitted values of θ 4 , θ 5 ranging roughly from − 2 to 2 and fitted values of θ 6 ranging rougly from 0.3 to 6. We compare the true option prices, obtained by means of two-dimensional quadrature, with the lower bounds obtained by applying a one-dimensional quadrature to (5). We present the results in Table 1. The approximation provided by the lower bound deteriorates as the option goes deep out of the money. This effect is clearly expected and often encountered in such kinds of approximation [see, for instance, the extensive numerical experiments in Caldana and Fusai (2013)]. Nevertheless, the approximation quality is satisfactory and unaffected by the parameters governing co-movements.

Applications to option pricing
The approximation provided in Proposition 2 allows for fast estimation of the risk-neutral joint probability distribution. In this section, we evaluate the quality of fit when estimation is performed on observed spread option prices. We also provide empirical evidence of the fact that the fundamental uncertainty implied from spread option prices can be far from being normally distributed.
Our dataset consists of WTI-Brent spread option closing prices along with WTI and Brent future quotes at NYMEX, downloaded from Bloomberg. This option pays at maturity  Table 2, we report the number of options available for different levels of moneyness and maturities. The time period spans from June 3, 2019, to April 30, 2020, for a total of 231 trading dates. For each trading date, we have a set of prices of spread options with different strike prices and maturities and the corresponding value of the underlying futures. We estimate the probability distribution for each of the 231 dates available, by selecting the parameters σ 1 , σ 2 , ρ, θ 4 , θ 5 , θ 6 that minimize the sum of squared difference between the observed option prices and the theoretical option prices from our specification. Then, we evaluate the quality of fit of the model by both in terms of Mean Absolute Error (MAE), that is the averaged absolute value of the difference between observed and theoretical option prices and in terms of Mean Relative Error (MRE), defined as the averaged absolute value of the relative error. In Fig. 2, we plot the calibration error in terms of MAE for each of the 231 days of calibrated option prices. The figure clearly shows the ability of the considered probability distribution to replicate observed option prices. The deterioration of the quality of fit observed in the last part of the plot reflects the turbulence observed in the oil market due to the start of the Covid-19 pandemic in february 2020. In that period, the futures on WTI experienced a dramatic crash due to the lack of storage facilities. Despite this fact, the loss in quality of fit seems to be quite moderate if compared with the tremendous turmoil observed in  the market during those months. This bring us to the conclusion that our specification is able to reproduce observed option prices even in periods of unprecedented turbulence. In Table  3, we report the MREs, classified by maturities and levels of moneyness. The table reveals a satisfactory fit in terms of relative error. Higher levels of moneyness seem to deteriorate the goodness of fit of the model, although the deterioration seems to be restrained. The impact of different levels of maturity on the quality of the fit appears to be random, as no clear patterns are identifiable. Overall, our analysis leads us to conclude in favor of the particular specification of the generalized bivariate normal distribution proposed in this paper, in terms of its ability to replicate observed option prices, even in presence of high turmoil in the market.
To illustrate how the WTI-Brent's co-movements are affected by high-order co-moments, we select 12 trading dates across the period under investigation and show the implied density functions of the couple (Y 1 , Y 2 ) estimated during those dates. The couple drives the fundamental uncertainty in the market. 3 In Fig. 3, we plot the first set of estimated density functions. These plots refer to the period of time between September and November, 2019. Despite the fact that during those months the oil market was experiencing a time of relative calm, the density functions implied by our specification of the generalized bivariate normal distribution with high-order co-moments are far from following a bivariate normal distribution. In particular, all six panels show signs of abundant (risk neutral) co-skewness. Co-volatility is also abundant, especially in the last four panels of Fig. 3. In Fig. 4, we plot the second set of implied density functions. The second row of the figure is particularly interesting, since it refers to the months in which the oil market was experiencing the tremendous crash provoked by the fall of the demand subsequent to the world-wide pandemic. The implied density functions show patterns in which co-skewness and co-volatility are prominent, showing that high-order (risk-neutral) co-moments are important transmission channels of contagion in the WTI-Brent market during the last turbulent period.

Generalized normal distribution and risk management
In this section we provide a further applications aimed at evaluating the impact of comovements on risk management, in terms of risk measures like Value-at-Risk (VaR) and Expected Shortfall (ES).
In what follows, we consider two stocks S 1 , S 2 , whose future behavior after τ days is modeled, under the historical probability measure, by: where the density of the bivariate random variable (Y 1 , Y 2 ) is given in (1). 4 We set the time horizon to τ = 1, denote by X 1 , X 2 the logarithmic returns of the two stocks, and focus on risk measures for the portfolio with return described by the random variable where ω ∈ [0, 1] is the weight of the first asset. The main goal of the section is to assess whether the presence of co-movements increases the two risk measures compared to the normal case. In particular, in both cases, we compute for the portfolio the 5%-VaR and 5%-ES and report the relative percentage differences V a R 1 where the risk measures with superscript '1' and 'N' are the ones from model (1) and the normal model, respectively. We investigate the effect of coskewness on VaR/ES and therefore make the θ 4 parameter vary in the interval (−0.95, 0.95) and fix the remaining parameters as follows: ρ = θ 5 = 0, and θ 6 = 0.5. For the normal distribution we assume zero correlation and hence this case is also obtained from (1) setting ρ = θ 5 = θ 5 = θ 6 = 0. In addition, we set μ 1 = μ 2 = 0 and σ 1 = σ 2 = 1. This allows us to isolate the effects of co-movements, in line with the objectives of our paper. We consider an equally-weighted portfolio (ω = 0.5), one with a weight ω = 0.25 and one with ω = 0.75. 5 Results, reported on Fig. 5, are obtained using the Monte Carlo method to estimate the two risk measures, based on the average of 10,000 estimates (each obtained from a sample of length 1,000). For both ω = 0.5 and ω = 0.25 the graphs are u-shaped. In particular, when ω = 0.25, the two risk measures are severely underestimated by the normal model when the value of the parameter associated to co-skewness is larger (in absolute value) than one half. Furthermore, in Panel (b) the line corresponding to ES lies consistently above the one for VaR, meaning that in this case the normal model underestimates ES more severely, at least for values of θ 4 near plus or minus one. Interestingly, when ω = 0.75, the relative percentage differences are decreasing with θ 4 . The two risk measures calculated for this portfolio are larger for the model allowing for co-movements only for values of θ 4 very close to minus one.
We also make a comparison in terms of the following conditional Expected Shortfall The above quantity has been recently investigated in the context of the bivariate Kumaraswamy distribution by Ghosh and Marques (2021) and is particularly useful to measure the overflow of the two risks X 1 and X 2 at the same time.
We calculate relative percentage differences between model (1) and the normal one. For the former we make the θ 4 parameter vary in the interval (−0.95, 0.95) and set θ 5 = 0.75 and θ 6 = 0.5. Also, to isolate the effect of co-movements, again we set μ 1 = μ 2 = 0 and σ 1 = σ 2 = 1. Relative percentage differences are given in Fig. 6. For both models, we consider three different values for the ρ parameters, i.e. ρ ∈ {−0.5, 0, 0.5}. It is clear that for ρ = 0.5 the conditional Expected Shortfalls from the normal model are always larger than the measures based on the bivariate density (1). However, for ρ = 0 and θ 4 < 0 and for ρ = −0.5 and θ 4 < 0.5 we observe that the considered risk measure is larger for model (1).
To complete this section, we use real data to evaluate the impact of co-movements on risk measures like VaR and ES. The aim is to compare risk measures for the portfolio return (7) and (X 1 , X 2 ) described by density (1) and those based on a normal model. We consider a pair of assets, namely the S&P 500 index (Ticker SPX) and gold (Ticker XAU). The rationale for selecting these two assets is that, contrary to the previous case, they allow diversification in a portfolio. We download 1,231 daily observations for log-returns of the two assets and assume that (X 1 , X 2 ) are the log-returns from SPX and XAU, respectively and are modelled either as the one-day-ahead log-returns from Eq. (6) or via the bivariate normal distribution. We estimate the parameters of the model with the maximum likelihood (ML) method using daily data and a rolling window of 1000 observations. The set of estimated parameters covers the between risk measures for the portfolio return (7) based on log-returns from equation (6) and those from a bivariate normal distribution same period of the previous analysis, i.e. June 3, 2019 to April 30, 2020. We use the estimated parameters to derive portfolio VaR and ES and do the same in the normal case. In Figure 7 we plot the time series of relative percentage differences between risk measures calculated based on (1) and based on a bivariate normal distribution. A few remarks are in order. Firstly, for all the days in the sample and for all the values of ω we consider, risk measures that take into account co-movements are larger than the normal-based measures. Indeed, VaR and ES based on (1) and on the estimated parameters are about 11% to 26% larger that the same Fig. 6 Relative percentage differences C E S 1 5%,5% (X 1 , X 2 )/C E S N 5%,5% (X 1 , X 2 ) × 100 between conditional ES based the log-returns from Eq. (6) and those from a bivariate normal distribution

(red lines) and
E S 1 5% (R p )/E S N 5% (R p ) − 1 × 100 (black lines) between risk measures for the portfolio return R p = ωX 1 + (1 − ω)X 2 . The risk factors (X 1 , X 2 ) denote log-returns from SPX and XAU, respectively and are modelled either as the one-day-ahead log-returns from Eq. (6) or via the bivariate normal distribution

Conclusions
In this paper, we studied a specification of the generalized bivariate normal distribution which combines great flexibility with simple and fast computations for option pricing. The considered probability distribution allows to take into accounts high-order co-moments, namely co-skewness and co-volatility, when modeling financial uncertainty. It is a special case of Fry-McKibbin et al. (2014), but has the important advantage of being tractable. Indeed, by conditioning arguments, most important quantities of interest can be expressed in terms of one dimensional integral, which in turn can be evaluated quickly.
We derive a closed form formula for exchange option prices, thus generalizing Margrabe (1978), and an approximated formula for spread option prices and greeks which turns out to be very accurate. We exploited the gain in speed [with respect to the more general model of Fry-McKibbin et al. (2014)], due to the tractability of the density function, to perform 231 calibration exercises, one for each trading dates in the available dataset, on WTI-Brent spread options traded at NYMEX. Our exercise shows that: i) The probability distribution provides a good fit for spread option data also in turbulent periods and ii) The implied density functions display patterns in which co-volatility and co-skewness play an important role.
From the marginal density of Y 1 another possible representation for the constant η is found: The condition θ 2 4 − 2θ 2 θ 6 < 0 is needed for the conditional variance σ 2 2|1 (y 1 ) to be strictly positive ∀y 1 .
Putting the two expectations together gives the desired result.

B Greeks
In this section we report the Greeks for spread options.