# Spread Option Pricing in a Copula Affine GARCH(p,q) Model

Edoardo Berton  Lorenzo Mercuri  CREST Japan Science and Technology Agency.

###### Abstract

In this study, we construct a bivariate market model combining the copula
function with the affine GARCH(p,q) process used to describe the marginal
dynamics of the log price. We then provide a numerical procedure for pricing
European spread option contracts. To assess the accuracy of our approach we
present a comparison with the Monte Carlo simulation method.

keywords: Copula Function, Affine GARCH(p,q) process, Spread Options

##  1 Introduction

A spread option is a contract written on the price difference between two
assets, referred to as spread and defined as $S_{1,t}-S_{2,t}$ for
$t=0,...,T$. Accordingly, the payoff of a European call spread option can be
written as $\max\\{S_{1,T}-S_{2,T}-K;0\\}$. Although few spread options are
traded on large exchanges, most of the transactions take place over-the-
counter (OTC). Spread options are particularly suited for commodity markets,
where their features make them remarkably appropriate. The OTC nature of the
contract does not fit ordinary speculation purposes but rather allows to hedge
certain risks embedded in the features of commodity markets. Notably,
“processing” spreads allow the trader to reduce or gain exposure on a
commodity processing routine. For instance, a “Crack Spread” option (an option
whose underlying is the spread between crude oil and the petroleum product
refined from it) allows a refiner to cover against unfavourable fluctuations
in the prices of input or output.

The most common approach in the literature, as well as the prevailing in
practice, is the Multivariate Geometric Brownian Motion (henceforth MGBM)
approach. MGBM model is often used alongside Monte Carlo simulation and
provides the framework for further developments. The approach relies on
Pearson’s $\rho$ to collect full information about the processes’ dependence
structure. Many authors have provided their contribution to this approach by
deriving a closed form solution to the pricing problem. Notably, [19] derives
a formula to price options when the $K=0$. [16] provides a further
approximation that is currently the market standard for the contract. [5]
compute an analytical approximation by reducing the price range to a bounded
interval that is particularly tight for certain parameters values. Finally,
[1] introduce a lower bound for the price that has been proven to be more
accurate than the approximation in [16].

The vivacity of the market for spread options has recently drawn much
attention. Several works have provided new theoretical insights on the topic
of spread options. [4] and [27] present new methodologies for the computation
of the price of spread options, while [21] provides an assessment of spread
option price sensitivities through Malliavin calculus. Further enrichments to
the spread option framework are introduced by [26] and [25]. The former
proposes a joint model for day-ahead electricity prices in interconnected
markets and observes that the model yields analytical prices. The latter
develops a multi-factor model and shows how to compute the joint
characteristic function of a pair of futures prices.

An alternative method to price spread options by means of Copula functions is
provided in [6] and [12]. These works propose a Copula-based dependence
structure and different processes for underlying assets. Previous literature
exhibits two major limitations: (i) past pricing methods rely heavily on
Gaussian framework, even though different processes may better reflect the
underlying’s behaviour, and (ii) most of the literature (although not all)
imposes linear correlation as a measure of dependence.

In this study, we develop a framework for spread option pricing that relies on
the application of affine GARCH(p,q) models. The class of GARCH models is
introduced in [2] as a useful adaptation of the ARCH model presented in [11].
GARCH models assume that the volatility process at time $t$ is influenced not
only by previous error terms but also by its lagged values, favouring
volatility clustering dynamics. Another advantage of GARCH modelling is that
it does not require returns to be independent. Indeed, despite being
considered Gaussian distributed conditional to past values, returns’
unconditional distribution exhibits leptokurtosis.

From the seminal work in [10], several GARCH models for option pricing have
been proposed in the financial literature. A major breakthrough occurred with
the GARCH model with Gaussian innovations developed in [14] where the
specification of the variance process yields a recursive procedure for the
computation of the characteristic function of the log price. Therefore the
price of a European option can be obtained by means of the Fourier Transform
as in [13]. To distinguish this model from to the original GARCH, in
literature the model in [14] is known as ”Affine” GARCH model with Gaussian
innovation. Affine Garch model with non-guassian innovations are the model
proposed in [7] with Inverse Gaussian innovations and the model in [20] with
Tempered Stable Innovations.

We embed the dependence structure of the processes within a Copula function.
The Copula-based setup provides a flexible method of capturing the dependence
between the marginals without loss of fit on other crucial features of
financial data (heavier tails, non linear dependence and asymmetrical
dependence). Considering an appropriate Copula family, we derive a formula for
a spread option that involves the first partial derivatives of the Copula
function. Numerical results of this formula are benchmarked against an
adequate Monte Carlo simulation to assess the accuracy and the performance of
the introduced method.

This paper is organized as follows. Sections 2 and 3 review the main results
about the affine GARCH model with normal innovations and the main properties
of the copula functions respectively. In Section 4 we propose a numerical
method for pricing the spread options based on the evaluation of a single
integral. Section 5 assesses the numerical precision of the formula proposed
in Section 4 assuming that the dependence structure between assets is
described by a Plackett copula and the univariate distribution of log-returns
can be generated from two Heston-Nandi GARCH models. Section 6 concludes the
paper.

##  2 Heston-Nandi GARCH model with Gaussian Innovations

In this section we review the affine GARCH model with normal innovations
proposed in [14], henceforth referred to as HN-GARCH.

This specification retains the features of the GARCH introduced in [2] and
overcomes the drawback related to symmetrical effects of shocks. In addition,
it is able to capture simultaneously the stochastic nature of volatility and
its negative correlation with spot returns, enabling quick adjustments in
variance dictated by changes in market levels.

Continuous-time stochastic volatility models that have been previously
introduced, see [8],[13], are difficult to implement and computationally
intensive, whereas the HN-GARCH process is much easier to apply to available
data, and volatility is readily observable from the history of asset prices.
Besides, the GARCH process defined in [14] allows option pricing in a semi-
analytical closed-form formula, leveraging to a certain extent on numerical
integration. Lastly, it can be proven that, as the observation interval
shrinks, the results are numerically close to the continuous-time stochastic
volatility model of [13].

The model has two basic assumptions. The first pertains the specification of
the GARCH that models the volatility of the spot asset log price as well as
the process for the spot log price itself.

######  Assumption 1.

The spot asset price, $S_{t}$, follows the process below over time steps of
length $\Delta$:

$log(S_{t})=log(S_{t-\Delta})+r_{f}+\lambda h_{t}+\sqrt{h_{t}}z_{t}$ (1a)

where:

$h_{t}=\omega+\sum_{j=1}^{p}\beta_{j}h_{t-j\Delta}+\sum_{j=1}^{q}\alpha_{j}(z_{t-j\Delta}-\gamma\sqrt{h_{t-j\Delta}})^{2}$
(1b)

Where $r_{f}$ is the continuously compounded risk-free interest rate over the
time interval $\Delta$ and $z_{t}$ is a standard normal noise. $h_{t}$ is the
conditional variance of the log returns between time $t-\Delta$ and time $t$
and is known from the information set at time $t-\Delta$.

We can see in (1a) that the conditional variance $h_{t}$ is included in the
mean as a risk premium. Therefore the expected spot return is proportional to
the variance by a parameter $\lambda$111Since volatility equals
$\sqrt{h_{t}}$, the return premium per unit of risk is proportional to the
square root of $h_{t}$ as well exactly as in [9].. Successive analysis and
implementation will focus on the ($p=1,q=1$) specification of the HN-GARCH,
deemed the most parsimonious yet effective in modeling the time series. The
process therefore includes one lag of both the variance and the error term,
and can be written as

$h_{t}=\omega+\beta h_{t-1}+\alpha(z_{t-1}-\gamma\sqrt{h_{t-1}})^{2}$ (2)

The expected value of variance at period $t+1$ is given by

$E[h_{t+1}]=\omega+\alpha+(\beta+\alpha\gamma^{2})h_{t}$ (3)

Imposing the constraint $\beta+\alpha\gamma^{2}<1$, the first order process is
mean reverting, with finite mean and variance. Further, the process will tend
to a long-time unconditional variance level of

$V_{L}=\frac{\omega+\alpha}{1-\beta-\alpha\gamma^{2}}$ (4)

The process for $h_{t+\Delta}$ is fully observable at time $t$ as a function
of the spot price

$h_{t+\Delta}=\omega+\beta h_{t}+\alpha\frac{(R_{t}-r_{f}-\lambda h_{t}-\gamma
h_{t})^{2}}{h_{t}}$ (5)

where $R_{t}$ is the log return at time $t$ defined as
$log(S_{t})-log(S_{t-\Delta})$. The $\alpha$ parameter determines the kurtosis
of the distribution, and imposes deterministic time-varying variance if set
equal to zero. $\gamma$ accounts for asymmetric impact of shocks; a large
negative shock increases the variance more than a large positive shock,
consistently with financial theory.

The covariance between the variance process and the spot price is defined by
the following equation

$Cov(h_{t+\Delta},log(S_{t}))=-2\alpha\gamma h_{t}$ (6)

Given the non-negativity constraint imposed on the $\alpha$ parameter, a
positive value of $\gamma$ would imply for the spot price and variance to
exhibit negative correlation.

We report the following result provided in [14].

######  Proposition 1.

The value of a call option with one period to expiry obeys the Black-Scholes-
Rubinstein formula.

This result implies that the conditional risk neutral distribution is still
normal and, as in [10], the risk neutral dynamics for the log price and for
the variance can be rewritten as follows:

$log(S_{t})=log(S_{t-\Delta})+r_{f}-\frac{1}{2}h_{t}+\sqrt{h_{t}}z^{\ast}_{t}$
(7a) $h_{t}=\omega+\beta
h_{t-1}+\alpha(z^{\ast}_{t-1}-\gamma^{\ast}\sqrt{h_{t-1}})^{2}$ (7b)

where

$z^{\ast}_{t}=z_{t}+(\lambda+\frac{1}{2})\sqrt{h_{t}}$
$\gamma^{\ast}=\gamma+\lambda+\frac{1}{2}$ $\lambda^{\ast}=-\frac{1}{2}.$

Figure 1: Simulated distributions of Heston-Nandi GARCH prices over different
periods. Parameters: $S_{0}=50.52\$$, $\sigma^{2}=0.036\%$, $\alpha=7.5e-6$,
$\beta=0.91$, $\gamma=91.98$, $\omega=5.4e-30$.

###  2.1 Heston-Nandi Closed Form Option Pricing

The main advantage of the introduction of this model is to allow closed-form
option pricing. That is, compute option price based on the analytical
martingale distribution of the underlying at expiry, rather than through
numerical procedures such as Monte Carlo methods.

Note that the conditional generating function of the spot price at time $t$ is
computed mathematically as

$G(u)=\mathsf{E_{t}}[S_{T}^{u}]$ (8)

Considering that for a general process X the moment generating function is
defined as

$MGF_{X}(u)=\mathsf{E}[e^{uX}]$ (9)

we conclude that the conditional generating function of the spot price $S_{T}$
in (8) is also the moment generating function of the logarithm of $S_{T}$.

In order to obtain the probability density function and probability
distribution of the log price under the risk neutral measure $\mathbb{Q}$ at
expiry, [14] obtains a log linear form for the moment generating function of
the log price, that we report in the following proposition for the HN-
GARCH(1,1) model.

######  Proposition 2.

The moment generating function of the log price has the following form

$MGF_{\ln(S_{T})}(u)=S_{t}^{u}\exp{(A_{t}+B_{t}h_{t+\Delta})}$ (10)

Where $S_{t}$ is the spot price of the asset, and

$A_{t}=A_{t+\Delta}+ur_{f}+B_{t+\Delta}\omega-\frac{1}{2}\log{(1-2\alpha
B_{t+\Delta})}$ (11a) $B_{t}=u(\lambda+\gamma)-\frac{1}{2}\gamma^{2}+\beta
B_{t+\Delta}+\frac{1}{2}\frac{(u-\gamma)^{2}}{1-2\alpha B_{t+\Delta}}$ (11b)

Coefficients $A_{t}$ and $B_{t}$ can be derived recursively from the terminal
conditions $A_{T}=B_{T}=0$.

From (10), it is possible to obtain the characteristic function by
substituting $u$ with $iu$. Therefore the characteristic function of the
logarithm of the spot price reads:

$\phi_{X}(u)=S_{t}^{iu}\exp{(A_{t}+B_{t}h_{t+\Delta})}$ (12)

The characteristic function, under the martingale measure, can be obtained by
setting $\gamma=\gamma^{*}$ and $\lambda=-\frac{1}{2}$.

Given the risk neutral characteristic function, $\phi_{X}^{*}(u)$ of the
logarithm of the spot price, the probability density function $f_{X}(x)$ and
the cumulative distribution function $F_{X}(x)$ can be obtained through
inverse Fourier transform as follows:

$f_{X}(x)=\frac{1}{\pi}\int_{0}^{\infty}\mathbb{R}\left[e^{-iux}\phi_{X}(u)\right]du$
(13a)
$F_{X}(x)=\frac{1}{2}-\frac{1}{\pi}\int_{0}^{\infty}\mathbb{R}\left[\frac{e^{-iux}\phi_{X}(u)}{iu}\right]du$
(13b)

Where $\mathbb{R}[z]=\frac{z+\bar{z}}{2}$ represents the real part of a
complex number $z$.

Figure 2: Characteristic function $\phi_{X}(u)$ and derived probability
density function and distribution function over 90 days. Parameters:
$S_{0}=80$, $\sigma^{2}=0.5\%$, $\alpha=7e-5$, $\beta=0.56$, $\gamma=0.45$,
$\omega=1.6e-5$.

Figure 2 represents the characteristic function (Figure 2a) and the derived
probability density (Figure 2b) and distribution function (Figure 2b). The
Inverse Fourier Transform can be used to price a European call option
following the approach developed first in [13].

##  3 Copula Function and Dependence Structure

In this section we briefly report the theoretical concept of bivariate Copula,
outlining the main theorems alongside some of the principal properties. We
then describe the most relevant copula families and their applications.
Further, we introduce the Archimedean Copulas and give some examples of the
most significant specifications. Finally, we review the notion of concordance
measure and define the two most common measures: Kendall’s $\uptau$ and
Spearman’s $\rho$.

######  Definition 1.

A function $C:\left[0,1\right]^{2}\rightarrow\left[0,1\right]$ is called a
bivariate Copula if

$C\left(u,v\right)=\mathsf{Pr}\left(U\leq u,V\leq v\right)$ (14)

where $U,V$ are two independent Uniform random variables.

In order to be a valid Copula, a function must satisfy certain properties that
follow from probability theory. Specifically, a Bivariate Copula must comply
with the following three properties:

  1. 1.

Normalized marginals: $C(u,1)=C(1,u)=u$. From the uniform marginals property,
it holds that $\mathsf{Pr}(U\leq u,V\leq 1)=\mathsf{Pr}(U\leq u)=u$ and vice
versa.

  2. 2.

Groundedness: $C(u,v)$ is non-decreasing in each element. If at least one
component is equal to zero, then $C(u,v)=0$.

  3. 3.

2-increasingness: for each $a,b,c,d\in\mathbf{I}$, with $a\leq b$ and $c\leq
d$, $\mathsf{Pr}(U\in\left[a,b\right],V\in\left[c,d\right])$ must be non
negative. Hence

$0\leq C(b,d)-C(a,d)-C(b,c)+C(a,c)\leq 1$ (15)

The above properties can be extended to any $d$-dimensional Copula, although
for the sake of simplicity we address only the bivariate specification.

Let us now recall the fundamental theorem of Copula theory. Sklar’s Theorem
provides the theoretical foundations for the application of Copulas and shows
how joint distributions can be set up by a number of marginal distributions
and a copula function $C$ that specifies the dependence structure among them.

######  Theorem 1.

A function F: $\mathbb{R}^{2}\rightarrow[0,1]$ is the distribution function of
some bivariate random vector $(X_{1},X_{2})$ if and only if there exist a
Copula $C:[0,1]^{2}\rightarrow[0,1]$ and univariate distributions
$F_{1},F_{2}:\mathbb{R}\rightarrow[0,1]$ such that

$C(F_{1}(x_{1}),F_{2}(x_{2}))=F(x_{2},x_{2})\;\;\;\;\;\;\;\;x_{1},x_{2}\in\mathbb{R}$
(16)

If the distribution function $F_{j}$ of all components $X_{j}$ with $j=1,2$ is
continuous, then the correspondence between C and F is uniquely determined.

An interesting implication of Sklar’s theorem is that, conditional on the
marginal distributions $F_{1}$ and $F_{2}$ being absolutely continuous, the
corresponding unique Copula distribution can be rewritten as

$C(u,v)=F(F_{1}^{-1}(u),F_{2}^{-1}(v))$ (17)

The latter illustrates how Copulas can be extracted from multivariate joint
cumulative distribution functions with continuous margins. Noting that
$C(u,v)$ corresponds to the probability that $X_{1}$ and $X_{2}$ lie below the
$u$-th and $v$-th quantile respectively, it appears clear how Copulas express
dependence on a quantile scale.

It is most often meaningful to define how the Copula function is distributed
conditional on the marginals.

######  Definition 2.

The conditional distribution function of the Copula $C(U,V)$ given $V=v$,
which we denote $c_{v}(u)$, is given by:

$c_{v}(u)=\lim_{\Delta v\rightarrow 0}\frac{C(u,v+\Delta v)-C(u,v)}{\Delta
v}=\frac{\partial C(u,v)}{\partial v}$ (18)

The conditional distribution $c_{v}(u)$ is a distribution in
$\left[0,1\right]$ and is Uniform if and only if the Copula corresponds to the
Independence Copula $\Pi$. The conditional distribution function is of
significant importance when simulating a random sample that is jointly
distributed according to the Copula $C$.

The second order mixed partial derivative of a Copula function provides the
joint density function, denoted $c(u,v)$. Formally, in the bivariate case

$c(u,v)=\frac{\partial^{2}C(u,v)}{\partial u\partial v}$ (19)

Not all Copulas posses, however, a closed form density function. For instance,
the Countermonotonicity Copula and Comonotonicity Copula are not absolutely
continuous and do not admit the computation of the density function. Yet, most
parametric specifications provide analytical densities to ease further
computations.

In the following section we mention the Fréchet–Hoeffding bounds and their
significant dependence structure implications, as well as the most relevant
families of Copulas.

######  Definition 3.

Consider any Copula $C(u,v)$, then it always holds that

$\max\\{u+v-1,0\\}\leq C(u,v)\leq\min\\{u,v\\}$ (20)

Unlike numeric dependence measures, a Copula is a mathematical object that
collects full information about the dependence structure without projecting it
onto a single number. The inequality in (20) defines the so called
Fréchet–Hoeffding bounds that represent the two extreme dependence structures
a Copula can map. Analogously to numeric dependence measures, these bounds may
be interpreted as the cases of extreme dependence.

A remarkable result of (20) is the family of Fundamental Copulas. These
Copulas arise as the functions that describe perfect positive dependence,
perfect negative dependence and independence.

| bivariate definition | dependence structure
---|---|---
Comonotonicity | $M(u,v)=\min\\{u,v\\}$ | extreme positive
Copula | dependence
Countermonotonicity | $W(u,v)=\max\\{u+v-1,0\\}$ | extreme negative
Copula | dependence
Independence | $\Pi(u,v)=uv$ | independence
Copula
Table 1: Fundamental Copulas in the bivariate case and the dependence
structure they represent.

Unlike Comonotonicity and Independence Copula, the Countermonotonicity Copula
does not exist for dimensions $d>2$. It is easy to see that, in the bivariate
case, the Comonotonicity Copula is the joint distribution of a random vector
$(U,V)$ where $V=T(U)$ is a strictly increasing transformation of $U$.
Diametrically, the Countermonotonicity Copula is the Copula that collects
information regarding the dependence structure of the random vector $(U,1-U)$.
Finally, it follows from Sklar’s Theorem that the components $(U,V)$ of a
random vector are independent if and only if their Copula is the Independence
Copula.

Copula functions can be explicitly defined or extracted from an existing joint
distribution by virtue of Sklar’s Theorem results. The family of Implicit
Copulas stems from the latter method and includes the Gaussian Copula and the
$t$ Copula. Copulas belonging to this family reproduce the dependence
structure of the existing corresponding joint distribution, although it is
often the case that these Copulas cannot be extracted in a simple closed form.
Table 2 reports the most common examples of Implicit Copulas.

| bivariate definition
---|---
Gaussian Copula | $C_{\rho}^{G}(u,v)=\mathbf{\Phi}_{\rho}(\Phi^{-1}(u),\Phi^{-1}(v))$
$t$ Copula | $C_{\nu,\rho}^{t}(u,v)=\textbf{t}_{\nu,\rho}(t_{\nu}^{-1}(u),t_{\nu}^{-1}(v))$
Table 2: Gaussian Copula and $t$ Copula in the bivariate case. $\rho$ is the
correlation of the two random variables, $\Phi(\cdot)$ represents the standard
normal distribution function and $\mathbf{\Phi}$ the joint distribution
function. $\mathbf{t}_{\nu,\rho}$ indicates the multivariate $t$ distribution
with $\nu$ degrees of freedom whereas $t_{\nu}$ denotes the standard
univariate $t$ distribution with $\nu$ degrees of freedom.

In contrast with Implicit Copulas, explicitly defined Copulas arise from the
need of a simple closed-form algebraic expression to compute joint
probability. Among these Explicit Copulas, Archimedean Copulas enjoy
considerable popularity in a number of practical applications in several
disciplines. This collection of Copulas has in fact a convenient algebraic
expression and is capable of capturing asymmetry.

######  Definition 4.

A Copula $C_{\varphi}$ is an Archimedean Copula if it possesses the functional
form

$C_{\varphi}(u,v)=\varphi(\varphi^{-1}(u),\varphi^{-1}(v))$ (21)

For a suitable non-increasing function $\varphi:[0,\infty)\rightarrow[0,1]$
with $\varphi(0)=1$ and $\lim_{x\rightarrow\infty}\varphi(x)=0$, called the
“Archimedean Generator”.

[17] postulates that $C_{\varphi}$ is a valid distribution if the archimedean
generator $\varphi$ is $d$-monotone. The generator $\varphi$ is said to be
d-monotone if it has derivatives of all orders in $(a,b)$ and if
$(-1)^{k}f^{(k)}(x)\geq 0$ for any $x\in(a,b)$ and for $k=0,1,...,d$. Each
completely monotone generator function $\varphi$ can be applied to define a
Copula of arbitrary dimension $d\geq 2$. Table 3 shows the most common
Archimedean Copula specifications with their generator function.

Although Archimedean Copulas are typically parametrized in terms of a function
$\varphi$, in practice one chooses a parametric family of Laplace transforms
governed by some real parameter $\theta$, so that $\varphi=\varphi_{\theta}$.
This allows the Copula to be completely parametrized by the single parameter
$\theta$, so that $C_{\theta}$ is equivalent to $C_{\varphi_{\theta}}$. Hence,
it is the parameter $\theta$ that provides the range of different dependence
structures. An evident drawback of this specification, considering non-
negativity of $\varphi$, is that this class of Copulas cannot model negative
dependence structures.

name | $\varphi_{\theta}(x)$ | $\varphi^{-1}_{\theta}(y)$ | $\theta\in$
---|---|---|---
Clayton Copula | $(1+x)^{-1/\theta}$ | $y^{-\theta}$ | $(0,\infty)$
Gumbel Copula | $\exp\\{-x^{1/\theta}\\}$ | $(-\log(y))^{\theta}$ | $[1,\infty)$
Ali-Mikhail-Haq Copula | $\frac{1-\theta}{e^{x}-\theta}$ | $\log(\frac{1-\theta}{y}+\theta)$ | $[0,1)$
Table 3: Archimedean Copulas with the respective generator function $\varphi$
and $\theta$ parameter range.

Copula functions allow for a shift in dependence measurement, in the following
section we stylize some concepts about rank correlation. Most multivariate
distributions rely on a dependence structure described by matrices of pairwise
linear correlation coefficients. Pearson’s $\rho$ represents in fact the most
common measure of correlation within a pair of random variables in theoretical
statistics as well as in financial applications. This simple yet fundamental
measure has some drawbacks when it comes to modelling multivariate
distributions. First of all, it is only able to capture pairwise correlation.
When dealing with a multivariate random vector with $d>2$ in fact, a
correlation matrix of dimension $d\times d$ needs to be computed, with at
least $d(d-1)/2$ correlations to be estimated. Secondly, the relationship
between the random components might be non-linear and thus be distorted by
Pearson’s $\rho$. Lastly, financial applications of $\rho$ show that this
measure may result in an incorrect estimation of risk.

In order to account for the broader intuition underpinning correlation, [24]
introduces the notion of Concordance Measure. A Concordance Measure is a
mapping I from the set of bivariate distributions with continuous marginals to
the interval [-1,1] satisfying the following set of conditions:

  * •

Symmetry: $I(X_{1},X_{2})=I(X_{2},X_{1})$;

  * •

Coherence: let $(X_{1},X_{2})$ and $(Y_{1},Y_{2})$ be two bivariate random
vectors having Copulas $C_{X}(u,v)$ and $C_{Y}(u,v)$ respectively, and these
are ordered as $C_{X}(u,v)\leq C_{Y}(u,v)$ for every $(u,v)\in[0,1]^{2}$. Then
$I(X_{1},X_{2})\leq I(Y_{1},Y_{2})$;

  * •

Independence: if $X_{1}$ and $X_{2}$ are independent, $I(X_{1},X_{2})=0$. The
opposite is not necessarily true;

  * •

Change in sign: $I(-X_{1},X_{2})=-I(X_{1},X_{2})$. Intuitively, a change in
sign inverts the ordering of one random variable, hence the concordance
measure inverts itself.

  * •

Convergence: if $\\{(X_{1}^{(n)},X_{2}^{(n)})\\}_{n\in N}$ is a sequence of
random vectors with associated sequence of joint distribution functions
$\\{F_{n}\\}_{n\in N}$ and $F_{n}$ converges point-wise to $F$, the joint
distribution of a random vector $(X_{1},X_{2})$, then

$\lim_{n\rightarrow\infty}I(X_{1}^{(n)},X_{2}^{(n)})=I(X_{1},X_{2})$

Following, we introduce the most common concordance measures, but let us first
define “concordant” a pair of observations $(x_{1},x_{2})$ and $(y_{1},y_{2})$
such that $(x_{1}-y_{1})(x_{2}-y_{2})>0$. Conversely, the same pair is said to
be “discordant” if $(x_{1}-y_{1})(x_{2}-y_{2})<0$.

######  Definition 5.

(Kendall’s $\uptau$) Let $(X_{1},X_{2})$ be a random vector with copula C as
joint distribution function. Then Kendall’s $\uptau$ is defined as

$\uptau=4\int_{0}^{1}\int_{0}^{1}C(u,v)\frac{\partial^{2}C(u,v)}{\partial
u\partial v}dudv-1=4E[C(U,V)]-1$ (22)

Kendall’s $\uptau$ is defined for a bivariate random vector $(U,V)$ with
continuous marginals, so that the joint distribution is uniquely defined by
the copula $C$, and hence the concordance measure does not reflect any
information about the single marginals. The intuition behind (22) can be
interpreted as the probability of observing a concordant pair minus the
probability of observing a discordant pair, formalized as follows

$\uptau=\mathsf{Pr}((X_{1}-Y_{1})(X_{2}-Y_{2})>0)-\mathsf{Pr}((X_{1}-Y_{1})(X_{2}-Y_{2})<0)$
(23)

For iid samples, $(X_{1}^{(1)},X_{2}^{(1)})...(X_{1}^{(n)},X_{2}^{(n)})$, the
empirical version of (22) is defined as

$\hat{\uptau}=\frac{\text{n. of concordant pairs}-\text{n. of discordant
pairs}}{\text{n. of total pairs}}$ $=\frac{\sum_{1\leq i<j\leq
n}\text{sign}[(X_{1}^{(j)}-X_{1}^{(i)})(X_{2}^{(j)}-X_{2}^{(i)})]}{n(n-1)/2}$
(24)

######  Definition 6.

(Spearman’s $\rho$) Given the random vector $(X_{1},X_{2})$ with continuous
marginal distributions $F_{i}$ with $i=1,2$, let the vector $(U,V)$ be the
transformed vector $(U,V)=(F_{1}(X_{1}),F_{2}(X_{2}))$. Note that the vector
$(U,V)$ has joint distribution function C. Spearman’s $\rho$ is defined as the
Pearson’s correlation coefficient of $(U,V)$, formally

$\rho_{S}=\text{Corr}(U,V)=\text{Corr}(F_{1}(X_{1}),F_{2}(X_{2}))$ (25)

Again, this dependence measure does not depend on the marginal laws
$F_{1},F_{2}$ as the influence has been removed by transforming the random
variables into Uniforms. Mathematically, (25) is equivalent [see [22] for more
details] to

$\rho_{s}=12-\int_{0}^{1}\int_{0}^{1}C(u,v)dudv-3$ (26)

Properties of Pearson’s $\rho$ are naturally inherited by Spearman’s $\rho$,
such as symmetry and values of 0 and 1 in case of independence and perfect
positive dependence respectively. Analogously to Kendall’s $\uptau$,
$\rho_{S}$ is increasing in the point-wise ordering of the Copulas, such that,
if $C_{1}(u,v)\leq C_{2}(u,v)$ for all $(u,v)\in[0,1]^{2}$, then
$\rho_{S,1}\leq\rho_{S,2}$. Let $\hat{\rho_{S}}$ denote the empirical
counterpart of Spearman’s $\rho$. $\hat{\rho_{S}}$ is somewhat related to the
empirical version of the correlation coefficient, but levels of the random
variables are substituted by ranks. Formally, consider the iid samples
$(X_{1}^{(1)},X_{2}^{(1)}),...,(X_{1}^{(n)},X_{2}^{(n)})$, the sample version
can be computed as

$\hat{\rho_{S}}=\frac{\sum_{i=1}^{n}(\text{rank}[X_{1}^{(i)}]-\frac{n+1}{2})(\text{rank}[X_{2}^{(i)}]-\frac{n+1}{2})}{\sqrt{\sum_{i=1}^{n}(\text{rank}[X_{1}^{(i)}]-\frac{n+1}{2})^{2}}\sqrt{\sum_{i=1}^{n}(\text{rank}[X_{2}^{(i)}]-\frac{n+1}{2})^{2}}}$
(27)

where rank$[X_{j}^{(i)}]$ represents the rank statistic of $X_{j}^{(i)}$ with
$j=1,2$.

To close this section, note that many families of Copulas can be parametrized
in terms of a single parameter $\theta$, which typically collects the
information about the dependence structure between the marginals. It is often
the case that one can express the parameter $\theta$ as a function of either
$\uptau,\rho_{S}$ or both.

In the next part, we introduce the Copula specification we rely on for
computations.

###  3.1 Plackett Copula

In the continuance of the study, the Plackett Copula is used to model the
dependence structure of two asset prices in order to price a spread option. It
seems then natural in this section to highlight the theoretical framework of
such Copula and to review its features.

It is useful to firstly illustrate the concepts of “contingency table” and
“odds ratio”, inasmuch as the Plackett Copula heavily leverages and extends
these concepts to continuous margins.

| up | down |
---|---|---|---
up | $a$ | $b$ | $a+b$
down | $c$ | $d$ | $c+d$
| $a+c$ | $b+d$ |
Table 4: 2$\times$2 contingency table

Let us suppose, for instance, that the contingency table in Table 4 shows the
number of days in a year where two given stocks’ closing prices are up/down
relative to the opening prices. For simplicity, let A denote the row-stock and
B the column-stock. It is equivalent to express $a,b,c,d$ in Table 4 as the
absolute number of days or as the percentage of the total days $n$, where
$n=a+b+c+d$.

Intuitively, $a+b$ and $c+d$ represent the unconditional probabilities for
stock A to close at a higher price or at a lower price respectively. The
“odds” for this stock price to increase are $(a+b)/(c+d)$, it can be
interpreted as the probability of a favourable outcome in terms of the number
of unfavourable outcomes. Analogously, the odds of stock B’s price to increase
are computed as $(a+c)/(b+d)$. The odds can also be computed on the
conditional probabilities of one stock’s development, given the other stock’s
performance. The conditional probability of Stock A closing with a higher
price given Stock B result is “up” is $a/(a+c)$ whereas the opposite
conditional probability will be computed as $c/(a+c)$. One could then compute
the odds as the ratio of one probability to the other.

To assess the dependence between row variable and column variable one should
calculate the “cross-product ratio” or “odds ratio” as

$\theta=\displaystyle\frac{\displaystyle\frac{a}{a+c}/\displaystyle\frac{c}{a+c}}{\displaystyle\frac{b}{b+d}/\displaystyle\frac{d}{b+d}}=\frac{ad}{bc}$
(28)

The resulting parameter $\theta$ captures the dependence structure of the
stylized distribution. Note that, when $\theta>1$ probability tends to
concentrate in the main diagonal (up-up and down-down) which intuitively
corresponds to positive dependence; conversely, a value of $\theta<1$ would
correspond to negative dependence. Lastly, $\theta=1$ implies independence,
suggesting that each observed entry is equal to its expected value.

[23] introduces a bivariate distribution that provides a pragmatic extension
of the concept of contingency table dependence, allowing for continuous
margins.

Let $X$ and $Y$ be continuous random variables with joint distribution $H$ and
marginal distributions $F$ and $G$ respectively. Let us now rewrite “up” and
“down” as $X\leq x$ and $X>x$, and analogously for $Y$ random variable. Table
4 can be rearranged as illustrated in Table 5.

| $Y\leq y$ | $Y>y$ |
---|---|---|---
$X\leq x$ | $H(x,y)$ | $F(x)-H(x,y)$ | $F(x)$
$X>x$ | $G(y)-H(x,y)$ | $1-F(x)-G(y)+H(x,y)$ | $1-F(x)$
| $G(y)$ | $1-G(y)$ |
Table 5: 2$\times$2 contingency table in the continuous margin domain.

The cross-product ratio is then rewritten as

$\theta=\frac{H(x,y)[1-F(x)-G(y)+H(x,y)]}{[F(x)-H(x,y)][G(y)-H(x,y)]}$ (29)

By means of Theorem 1, if the margins are continuous then their Copula is
uniquely determined, so one can substitute $H(x,y)$ with $C(F(x),G(y))$ and
$F(x),G(y)$ with $u,v$ respectively. Accordingly, $\theta$ will become

$\theta=\frac{C(u,v)[1-u-v+C(u,v)]}{[u-C(u,v)][v-C(u,v)]}$ (30)

Holding $\theta$ constant, solving for $C(u,v)$ yields the following

$C(u,v)=\frac{[1+(\theta-1)(u+v)]\pm\sqrt{[1+(\theta-1)(u+v)]^{2}-4uv\theta(\theta-1)}}{2(\theta-1)}$
(31)

Note that [18] proves that (31) with the “+” sign before the radical is never
a Copula, whereas with “-” sign it always is. Furthermore, one needs to
encircle the parameter $\theta$ for (31) to comply with theoretical
requirements. Consequently the analytical expression of the so-called
“Plackett Copula” is

$\hfilneg
C_{\theta}(u,v)=\frac{[1+(\theta-1)(u+v)]-\sqrt{[1+(\theta-1)(u+v)]^{2}-4uv\theta(\theta-1)}}{2(\theta-1)}\hskip
10000.0pt$ (32a) for $\theta>0$ and $\neq 1$. $C_{\theta}(u,v)=uv$ (32b) for
$\theta=1$.

The latter are 2-increasing and absolutely continuous. It can be shown by
proving that

$\frac{\partial^{2}C_{\theta}(u,v)}{\partial u\partial v}\geq
0\;\;\;\;\;\;\;\;u,v\in\mathbf{I}^{2}$

and

$C_{\theta}(u,v)=\int_{0}^{u}\int_{0}^{v}\frac{\partial^{2}C_{\theta}(s,t)}{\partial
s\partial t}dtds$

Moreover, Plackett Copulas form a comprehensive family of Copulas, ranging
continuously from the lower Fréchet–Hoeffding bound in the limit case
$\lim_{\theta\rightarrow 0}C_{\theta}(u,v)=W(u,v)$, to the upper
Fréchet–Hoeffding bound in the limit case
$\lim_{\theta\rightarrow\infty}C_{\theta}(u,v)=M(u,v)$. It is trivial to note
that when $\theta=1$, $C_{1}(u,v)=uv=\Pi(u,v)$. Another advantage of Plackett
Copula is that it possesses-closed form analytical probability density
function

$c(u,v)=\frac{\theta\left[1+(u-2uv+v)(\theta-1)\right]}{\\{\left[1+(\theta-1)(u+v)\right]^{2}-4uv\theta(\theta-1)\\}^{3/2}}$
(33)

The Plackett Copula lends itself particularly well to modelling dependence
structure of financial phenomena, and it is well spread in practice. As only
one parameter $\theta$ captures the entire dependence dynamic, it is notably
easy to estimate via Maximum Likelihood from past data. One could also
estimate $\theta$ by computing the observed cross-product ratio in the dataset
using some arbitrary quadrant separation.

Figure 3: Scatter plot of Plackett Copula in the following cases: Negative
dependence ($\theta<1$), Independence ($\theta=1$) and Positive dependence
($\theta>1$).

The scatter plots in Figure 3 show the pattern in the joint distribution, when
$\theta$ ranges from $\theta=0.04$ to $\theta=30$. The two lateral plots
exhibit clear negative and positive dependence respectively, whilst the middle
one represents independence. Furthermore, one can express Spearman’s $\rho$ as
a function of $\theta$ as follows

$\rho_{S}=\left\\{\begin{aligned}
&\frac{\theta+1}{\theta-1}-\frac{2\theta}{(\theta-1)^{2}}\log{\theta}\;\;\;\;\;\;\;\;\text{if
}\theta\neq 1\\\
&0\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\text{if
}\theta=1\end{aligned}\right.$ (34)

Figure 4: Relationship between Plackett’s $\theta$ and Spearman’s $\rho$.

###  3.2 Plackett’s $\theta$

The relationship between the marginal distributions under the Plackett Copula
is entirely captured by the non-negative $\theta$ parameter. A simple yet
attractive approach to $\theta$ estimation is the observed cross-product
ratio, previously described in (28), $\theta^{\ast}=ad/bc$ where $a,b,d$ and
$c$ are the observed frequencies in the four quadrants determined by lines in
$\mathbb{R}^{2}$ parallel to the axes through a point $(p,q)$. An optimum
choice for $(p,q)$ is presented in [18] to be the sample median vector, which
in turn minimizes the asymptotic variance of the estimator $\theta^{\ast}$.
Considering the following contingency table, where $p$ is the median of
variable $X$ and $q$ is the median for variable $Y$, let $a,b,c$ and $d$ be
the relative frequencies in the quadrants.

By the definition of median, it holds that $a+c=b+d=0.5$. And so it follows
that $c=(0.5-a)$ and $b=(0.5-d)$. Again, from the definition of median it must
hold also that $a+b=c+d=0.5$. We equate previous results and obtain that
$a=d$. Then, taking the cross-product ratio, we can substitute

$\theta^{\ast}=\frac{ad}{bc}=\frac{ad}{(0.5-a)(0.5-d)}=\frac{a^{2}}{(0.5-a)^{2}}$
(35)

where $a$ represents the observed frequency of observations in which neither
of the variables exceeds its median value.

##  4 Spread Option Pricing with a Copula Function

In this section we propose a method for pricing a spread option based on a
numerical integration of one-dimensional integral. We assume that the joint
distribution of log-returns is modelled through a Copula function. For sake of
brevity we discuss only the case of a European spread call option. The formula
proposed in this section is easy to apply to any model where the univariate
cumulative distribution function has an analytical (or semi-analytical)
formula. Under a martingale measure $\mathbb{Q}$, the price of a spread call
option with strike price $K$ can be obtained as:

$P_{t_{0}}=e^{-r\left(T-t_{0}\right)}\mathsf{E}^{\mathbb{Q}}\left[\max\left\\{S_{1,T}-S_{2,T}-K;0\right\\}\left|\mathcal{F}_{t_{0}}\right.\right].$
(36)

Denoting with $X_{i,t_{0},T}$ with $i=1,2$ the log-return for each asset over
the period $\left[t_{0},T\right]$, the joint density function
$f_{X_{1,t_{0},T},X_{2,t_{0},T}\left|\mathcal{F}_{t_{0}}\right.}\left(x_{1},x_{2}\right)$
of log-returns given the information available at time $t_{0}$ can be written
using a Copula $C\left(u,v\right)$ function as:

$f_{X_{1,t_{0},T},X_{2,t_{0},T}\left|\mathcal{F}_{t_{0}}\right.}\left(x_{1},x_{2}\right)=c\left(F_{X_{1,t_{0},T}\left|\mathcal{F}_{t_{0}},\right.}\left(x_{1}\right),F_{X_{2,t_{0},T}\left|\mathcal{F}_{t_{0}},\right.}\left(x_{2}\right)\right)f_{X_{1,t_{0},T}\left|\mathcal{F}_{t_{0}}\right.}\left(x_{1}\right)f_{X_{2,t_{0},T}\left|\mathcal{F}_{t_{0}}\right.}\left(x_{2}\right)$

where the function $c\left(u,v\right)$ is the mixed second order derivative of
the Copula function $C\left(u,v\right)$, i.e.
$c\left(u,v\right):=\frac{\partial^{2}C\left(u,v\right)}{\partial u\partial
v}$; $F_{X_{i,t_{0},T}\left|\mathcal{F}_{t_{0}},\right.}\left(x_{i}\right)$ is
for asset $i$ the conditional univariate cumulative distribution function of
the log-returns and
$f_{X_{i,t_{0},T}\left|\mathcal{F}_{t_{0}}\right.}\left(x_{i}\right)$ is the
associated conditional density function. For any twice differentiable Copula
function, the formula in (36) can be obtained using the following result:

######  Theorem 2.

Let $\left(\Omega,\mathcal{F},\left(\mathcal{F}_{t}\right)_{t\geq
0},\mathbb{P}\right)$ be a filtered probability space. If the Copula function
is twice differentiable (i.e. $C\left(u,v\right)\in\mathcal{C}^{2,2}$), the
non-arbitrage price for the spread call option is:

$P_{t_{0}}=e^{-r\left(T-t_{0}\right)}\left[I_{t_{0},1}^{T}-I_{t_{0},2}^{T}-I_{t_{0},3}^{T}\right],$
(37)

where the quantities $I_{t_{0},1}^{T}$, $I_{t_{0},2}^{T}$ and
$I_{t_{0},3}^{T}$ are defined as:

$I_{t_{0},1}^{T}=\int_{d_{3}}^{1}S_{1,t_{0}}\exp\left(F^{-1}_{X_{1,t_{0},T}\left|\mathcal{F}_{t_{0}}\right.}\left(u\right)\right)\partial_{u}C\left(u,d_{2}\left(u\right)\right)\mbox{d}u$
$I_{t_{0},2}^{T}=\mathsf{E}^{\mathbb{Q}}\left[S_{2,T}\left|\mathcal{F}_{t_{0}}\right.\right]-\int_{0}^{1}S_{2,t_{0}}\exp\left[F^{-1}_{X_{2,t_{0},T}\left|\mathcal{F}_{t_{0}}\right.}\left(v\right)\right]\partial_{v}C\left(d_{1}\left(v\right),v\right)\mbox{d}v$
$I_{t_{0},3}^{T}=K\left[1-\int_{0}^{1}\partial_{v}C\left(d_{1}\left(v\right),v\right)\mbox{d}v\right]$

and the terms $d_{1}\left(v\right)$, $d_{2}\left(u\right)$ and $d_{3}$ can be
determined using the following three formulas

$d_{1}\left(v\right)=F_{X_{1,t_{0},T}\left|\mathcal{F}_{t_{0}}\right.}\left(\ln\left[\frac{S_{2,t_{0}}\exp\left[F^{-1}_{X_{2,t_{0},T}\left|\mathcal{F}_{t_{0}}\right.}\left(v\right)\right]+K}{S_{1,t_{0}}}\right]\right),$
$d_{2}\left(u\right)=F_{X_{2,t_{0},T}\left|\mathcal{F}_{t_{0}}\right.}\left(\ln\left(\frac{S_{1,t_{0}}\exp\left(F^{-1}_{X_{1,t_{0},T}\left|\mathcal{F}_{t_{0}}\right.}\left(u\right)\right)-K}{S_{2,t_{0}}}\right)\right)$

and

$d_{3}=F_{X_{1,t_{0},T}\left|\mathcal{F}_{t_{0}}\right.}\left(\ln\left(\frac{K}{S_{1,t_{0}}}\right)\right)$

where the univariate conditional distributions
$F_{X_{i,t_{0},T}\left|\mathcal{F}_{t_{0}}\right.}\left(\cdot\right)$ with
$i=1,2$ are determined under the martingale measure $\mathbb{Q}$;
$F^{-1}_{X_{i,t_{0},T}}\left(\cdot\right)$ with $i=1,2$ denote the
corresponding quantile functions.

###### Proof.

We start from the expectation in (36) that can be rewritten as:

$I_{t_{0}}^{T}=\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}\max\left\\{S_{1,t_{0}}e^{x_{1,t_{0},T}}-S_{2,t_{0}}e^{x_{2,t_{0},T}}-K;0\right\\}dF_{X_{1,t_{0},T},X_{2,t_{0},T}\left|\mathcal{F}_{t_{0}}\right.}\left(x_{1},x_{2}\right).$

where
$F_{X_{1,t_{0},T},X_{2,t_{0},T}\left|\mathcal{F}_{t_{0}}\right.}\left(x_{1},x_{2}\right)$
is the conditional joint cumulative distribution function under the martingale
measure $\mathbb{Q}$. Exploiting the representation in (16) and the twice
differentiability of the Copula function $C\left(u,v\right)$ we obtain the
following formula for the differential of the conditional joint distribution:

$dF_{X_{1,t_{0},T},X_{2,t_{0},T}\left|\mathcal{F}_{t_{0}}\right.}\left(x_{1},x_{2}\right)=c\left(F_{X_{1,t_{0},T}\left|\mathcal{F}_{t_{0}}\right.}\left(x_{1}\right),F_{X_{2,t_{0},T}\left|\mathcal{F}_{t_{0}}\right.}\left(x_{2}\right)\right)\mbox{d}F_{X_{1,t_{0},T}\left|\mathcal{F}_{t_{0}},\right.}\left(x_{1}\right)\mbox{d}F_{X_{2,t_{0},T}\left|\mathcal{F}_{t_{0}}\right.}\left(x_{2}\right)$
(38)

Function $c\left(\cdot,\cdot\right)$ is defined in (19). Substituting (38) in
the integral $I_{t_{0}}^{T}$, we obtain the following result:

$I_{t_{0}}^{T}=\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}\left(S_{1,t_{0}}e^{x_{1}}-S_{2,t_{0}}e^{x_{2}}-K\right)^{+}c\left(F_{X_{1,t_{0},T}\left|\mathcal{F}_{t_{0}}\right.}\left(x_{1}\right),F_{X_{2,t_{0},T}\left|\mathcal{F}_{t_{0}}\right.}\left(x_{2}\right)\right)\mbox{d}F_{X_{1,t_{0},T}\left|\mathcal{F}_{t_{0}},\right.}\left(x_{1}\right)\mbox{d}F_{X_{2,t_{0},T}\left|\mathcal{F}_{t_{0}}\right.}\left(x_{2}\right),$

where $\left(x\right)^{+}=x\mathbbm{1}_{\left\\{x>0\right\\}}$ and
$\mathbbm{1}_{\left\\{x>0\right\\}}$ is the indicator function that assumes
value 1 if $x>0$, 0 otherwise.
Defining
$u:=F_{X_{1,t_{0},T}\left|\mathcal{F}_{t_{0}},\right.}\left(x_{1}\right)$ and
$v:=F_{X_{2,t_{0},T}\left|\mathcal{F}_{t_{0}},\right.}\left(x_{2}\right)$, we
obtain:

$I_{t_{0}}^{T}=\int_{0}^{1}\int_{0}^{1}\left(S_{1}^{t_{0}}\exp\left[F^{-1}_{X_{1,t_{0},T}\left|\mathcal{F}_{t_{0}}\right.}\left(u\right)\right]-S_{2}^{t_{0}}\exp\left[F^{-1}_{X_{2,t_{0},T}\left|\mathcal{F}_{t_{0}}\right.}\left(v\right)\right]-K\right)^{+}c\left(u,v\right)\mbox{d}u\mbox{d}v.$
(39)

The integral $I_{t_{0}}^{T}$ can be decomposed into three terms:

$I_{t_{0}}^{T}=I_{t_{0},1}^{T}-I_{t_{0},2}^{T}-I_{t_{0},3}^{T}$ (40)

defined as follows:

$I_{t_{0},1}^{T}:=\int_{0}^{1}\int_{0}^{1}S_{1,t_{0}}\exp\left[F^{-1}_{X_{1,t_{0},T}\left|\mathcal{F}_{t_{0}}\right.}\left(u\right)\right]\mathbbm{1}_{\left\\{S_{1,t_{0}}\exp\left[F^{-1}_{X_{1,t_{0},T}\left|\mathcal{F}_{t_{0}}\right.}\left(u\right)\right]>S_{2,t_{0}}\exp\left[F^{-1}_{X_{2,t_{0},T}\left|\mathcal{F}_{t_{0}}\right.}\left(v\right)\right]+K\right\\}}c\left(u,v\right)\mbox{d}u\mbox{d}v$
(41)
$I_{t_{0},2}^{T}:=\int_{0}^{1}\int_{0}^{1}S_{2,t_{0}}\exp\left[F^{-1}_{X_{2,t_{0},T}\left|\mathcal{F}_{t_{0}}\right.}\left(v\right)\right]\mathbbm{1}_{\left\\{S_{1,t_{0}}\exp\left[F^{-1}_{X_{1,t_{0},T}\left|\mathcal{F}_{t_{0}}\right.}\left(u\right)\right]>S_{2,t_{0}}\exp\left[F^{-1}_{X_{2,t_{0},T}\left|\mathcal{F}_{t_{0}}\right.}\left(v\right)\right]+K\right\\}}c\left(u,v\right)\mbox{d}u\mbox{d}v$
(42)

and

$I_{t_{0},3}^{T}:=\int_{0}^{1}\int_{0}^{1}K\mathbbm{1}_{\left\\{S_{1,t_{0}}\exp\left[F^{-1}_{X_{1,t_{0},T}\left|\mathcal{F}_{t_{0}}\right.}\left(u\right)\right]>S_{2,t_{0}}\exp\left[F^{-1}_{X_{2,t_{0},T}\left|\mathcal{F}_{t_{0}}\right.}\left(v\right)\right]+K\right\\}}c\left(u,v\right)\mbox{d}u\mbox{d}v.$
(43)

We introduce the following function:

$d_{1}\left(v\right)=F_{X_{1,t_{0},T}\left|\mathcal{F}_{t_{0}},\right.}\left(\ln\left[\frac{S_{2,t_{0}}\exp\left[F^{-1}_{X_{2,t_{0},T}\left|\mathcal{F}_{t_{0}}\right.}\left(v\right)\right]+K}{S_{1,t_{0}}}\right]\right),$
(44)

we obtain this identity

$\mathbbm{1}_{\left\\{S_{1,t_{0}}\exp\left[F^{-1}_{X_{1,t_{0},T}\left|\mathcal{F}_{t_{0}}\right.}\left(u\right)\right]>S_{2,t_{0}}\exp\left[F^{-1}_{X_{2,t_{0},T}\left|\mathcal{F}_{t_{0}}\right.}\left(v\right)\right]+K\right\\}}=\mathbbm{1}_{u>d_{1}\left(v\right)}$

that we use to rewrite $I_{t_{0}}^{T}\left(2\right)$ and
$I_{t_{0}}^{T}\left(3\right)$ in the following way:

$\displaystyle I_{t_{0},2}^{T}$ $\displaystyle=$
$\displaystyle\int_{0}^{1}S_{2,t_{0}}\exp\left[F^{-1}_{X_{2,t_{0},T}\left|\mathcal{F}_{t_{0}}\right.}\left(v\right)\right]\left[\int_{d_{1}\left(v\right)}^{1}c\left(u,v\right)\mbox{d}u\right]\mbox{d}v$
(45) $\displaystyle=$
$\displaystyle\int_{0}^{1}S_{2,t_{0}}e^{F^{-1}_{X_{2,t_{0},T}\left|\mathcal{F}_{t_{0}}\right.}\left(v\right)}\left[1-\partial_{v}C\left(d_{1}\left(v\right),v\right)\right]\mbox{d}v$
$\displaystyle=$
$\displaystyle\mathsf{E}^{\mathbb{Q}}\left[S_{2,T}\left|\mathcal{F}_{t_{0}}\right.\right]-\int_{0}^{1}S_{2,t_{0}}\exp\left[F^{-1}_{X_{2,t_{0},T}\left|\mathcal{F}_{t_{0}}\right.}\left(v\right)\right]\partial_{v}C\left(d_{1}\left(v\right),v\right)\mbox{d}v$
$\displaystyle I_{t_{0},3}^{T}$ $\displaystyle=$ $\displaystyle
K\int_{0}^{1}\left[\int_{0}^{1}\mathbbm{1}_{u>d_{1}\left(v\right)}c\left(u,v\right)\mbox{d}u\right]\mbox{d}v$
(46) $\displaystyle=$ $\displaystyle
K\int_{0}^{1}\left[\int_{d_{1}\left(v\right)}^{1}c\left(u,v\right)\mbox{d}u\right]\mbox{d}v$
$\displaystyle=$ $\displaystyle
K\int_{0}^{1}\left[\partial_{v}C\left(1,v\right)-\partial_{v}C\left(d_{1}\left(v\right),v\right)\right]\mbox{d}v$
$\displaystyle=$ $\displaystyle
K\int_{0}^{1}\left[1-\partial_{v}C\left(d_{1}\left(v\right),v\right)\right]\mbox{d}v$
$\displaystyle=$ $\displaystyle
K\left[1-\int_{0}^{1}\partial_{v}C\left(d_{1}\left(v\right),v\right)\mbox{d}v\right]$

Now we compute the integral $I_{t_{0},1}^{T}$. The following identity holds
true:

$\mathbbm{1}_{\left\\{S_{1,t_{0}}\exp\left[F^{-1}_{X_{1,t_{0},T}\left|\mathcal{F}_{t_{0}}\right.}\left(u\right)\right]>S_{2,t_{0}}\exp\left[F^{-1}_{X_{2,t_{0},T}\left|\mathcal{F}_{t_{0}}\right.}\left(v\right)\right]+K\right\\}}=\mathbbm{1}_{\left\\{d_{2}\left(u\right)>v\right\\}}\mathbbm{1}_{\left\\{1>u>d_{3}\right\\}}$

where

$d_{2}\left(u\right):=F_{X_{2,t_{0},T}\left|\mathcal{F}_{t_{0}}\right.}\left(\ln\left(\frac{S_{1,t_{0}}\exp\left(F^{-1}_{X_{1,t_{0},T}\left|\mathcal{F}_{t_{0}}\right.}\left(u\right)\right)-K}{S_{2,t_{0}}}\right)\right)$
(47)

and

$d_{3}=F_{X_{1,t_{0},T}\left|\mathcal{F}_{t_{0}}\right.}\left(\ln\left(\frac{K}{S_{1,t_{0}}}\right)\right).$
(48)

The integral $I_{t_{0}}^{T}\left(1\right)$ becomes

$\displaystyle I_{t_{0}}^{T}\left(1\right)$ $\displaystyle=$
$\displaystyle\int_{d_{3}}^{1}S_{1,t_{0}}\exp\left(F^{-1}_{X_{1,t_{0},T}\left|\mathcal{F}_{t_{0}}\right.}\left(u\right)\right)\left[\int_{0}^{d_{2}\left(u\right)}c\left(u,v\right)\mbox{d}v\right]\mbox{d}u$
(49) $\displaystyle=$
$\displaystyle\int_{d_{3}}^{1}S_{1,t_{0}}\exp\left(F^{-1}_{X_{1,t_{0},T}\left|\mathcal{F}_{t_{0}}\right.}\left(u\right)\right)\partial_{u}C\left(u,d_{2}\left(u\right)\right)\mbox{d}u$

The result is obtained considering (49), (45), (46), (47) and (48). ∎

It is worth noting that the result in (37) can be applied to any Copula
function discussed in Section 3 and requires the evaluation of one-dimensional
integrals with bounded support that can be easily computed using several
integral techniques available in any programming language.

##  5 Numerical Analysis

In this section, we conduct two numerical analyses. First, we analyze the
computation efficiency and the accuracy of the result in (37). In the second
analysis, we compare our formula with a numerical integration method applied
for the evaluation of the double integral in (36). Moreover, we determine the
confidence interval at the $95\%$ level for the Monte-Carlo prices [see [3]
for more details] to assess whether the prices obtained applying our formula
belong to this interval.

We consider, as an example, the case of the Plackett Copula that has an
analytical expression for partial derivatives of any order. In order to
reproduce the heteroskedasticity observed in the financial time series, we
require that the univariate log prices, under martingale measure, follow two
HN Affine GARCH(1,1) models as described in Section 2. From this latter
assumption, the conditional univariate distribution of log-returns can be
determined using the Inverse Fourier Transform.

The dataset, that we use for the estimation of the parameters, is composed of
the futures prices of the contract nearest to expiry for Brent crude oil and
West Texas Intermediate (WTI) crude oil. Both contracts are quoted in U.S.
Dollars. Three years of futures contract prices is gathered on a daily basis
from March 1${}^{\text{st}}$, 2017 to February 29${}^{\text{th}}$, 2020. The
dataset consists of 774 closing price records for Brent and 792 closing price
records for WTI.

The Heston-Nandi GARCH models are estimated through Maximum Likelihood
Estimation on the series of log returns of both contracts. Table 6 shows the
estimated parameters of the two GARCH models and their respective standard
errors. The annualized volatility of the series, computed as
$\sqrt{252(\omega+\alpha)/(1-\beta-\alpha\gamma^{2})}$, is 0.2972 for Brent
futures returns and 0.2985 for WTI futures returns.Plackett’s parameter
$\theta$ is estimated using the result in (35) ($\theta^{\star}=50.52$).

$\omega$ $\alpha$ $\beta$ $\gamma$ $\lambda$ $\beta+\alpha\gamma^{2}$ Brent
9.124e-33 7.081e-6 0.914 96.505 -0.418 0.979 (183.187) (0.625) (4.415)
(270.264) (12.989) WTI 2.845e-4 7.155e-6 0.175 0.161 -0.522 0.175 (1.774)
(2.897) (12.587) (2.961) (1.322)

Table 6:  Estimates of the parameters of the GARCH models through MLE
(standard errors of the estimates).

For the usage of the formula in (37), the numerical integrals in (49), (45)
and (46) are computed applying the midpoint integration method that, in our
case, requires a partition of the support $[0,1]$. For each natural number
$N$, we divide the support $[0,1]$ into $N$ sub-intervals of the same length
$\Delta=\frac{1}{N}$. Thus, we approximate the integrals in (49), (45) and
(46) with a Riemann summation of the following form:

$\sum_{i=1}^{N}f\left(i-1+\frac{1}{2N}\right)\frac{1}{N}\rightarrow\int_{0}^{1}f(x)dx\text{
a.s. }N\rightarrow+\infty,$

where $f\left(x\right)$ is an integrable function in $[0,1]$. Table 7 analyzes
the accuracy and the computation time for the formula in (37) as $N$ varies.

K | $N=100$ | $N=500$ | $N=1000$ | $N=5000$ | $N=10000$
---|---|---|---|---|---
0.0 | 6.081 | 6.142 | 6.147 | 6.149 | 6.149
| (1.12) | (1.12) | (1.12) | (1.12) | (1.12)
2.5 | 3.964 | 4.014 | 4.017 | 4.019 | 4.019
| (1.11) | (1.11) | (1.12) | (1.12) | (1.12)
5 | 2.224 | 2.262 | 2.265 | 2.266 | 2.266
| (1.13) | (1.14) | (1.14) | (1.20) | (1.20)
7.5 | 1.073 | 1.1071 | 1.109 | 1.111 | 1.111
| (1.13) | (1.13) | (1.13) | (1.13) | (1.13)
10 | 0.493 | 0.528 | 0.539 | 0.533 | 0.533
| (1.12) | (1.12) | (1.12) | (1.13) | (1.13)
Table 7:  Comparison of European spread call option prices with time to
maturity 3 months. The prices are computed using the formula in (37) where the
integrals are evaluated by means of the midpoint integration approach. We
report the computational times in parenthesis and $N$ denotes the number of
sub-intervals in the partition of the support $[0,1]$.

We observe that computational time increases slowly as $N$ increases; all
prices seem to reach a stable behaviour for $N\geq 5000$.

We describe alternative approaches for pricing a spread European Call option
that we compare with our formula in (37). The first approach is based on the
numerical evaluation of the double integral in (39). The partition is now
associated to the set $\left[0,1\right]^{2}$ of the Cartesian plane. In this
case, the double integral can be approximated by a double Riemann summation of
the following form:

$\sum_{i=1}^{N}\sum_{j=1}^{N}f\left(i-1+\frac{1}{2N},j-1+\frac{1}{2N}\right)\frac{1}{N^{2}}\rightarrow\int_{0}^{1}\int_{0}^{1}f(x,y)dxdy\text{
a.s. }N\rightarrow+\infty,$

where $f(x,y)$ is a integrable function defined on the support
$\left[0,1\right]^{2}$. The second alternative approach is based on the Monte
Carlo simulation applying the following algorithm:

  1. 1.

We generate two independent random numbers $\hat{u}$, $\hat{s}$ from an
Uniform distribution with support $[0,1]$.

  2. 2.

We determine the simulated log-return $\hat{x}_{1}$ for the first asset as:

$\hat{u}=F_{X_{1,t_{0},T}|\mathcal{F}_{t_{0}}}(x_{1}).$

  3. 3.

We determine a random number $\hat{v}$ that satisfies the following equation:

$\hat{s}=c_{u}\left(v\right).$ (50)

$c_{u}\left(\cdot\right)$ is the partial derivative of the Copula function
with respect to the variable $u$ and corresponds to the conditional
distribution of $V$ given the event $U=u$. In the Plackett Copula, equation
(50) has an analytical solution as explicited in [15].

  4. 4.

We determine the simulated log-return $\hat{x}_{2}$ for the second asset
solving with respect to $x_{2}$ the following equation:

$\hat{u}=F_{X_{2,t_{0},T}|\mathcal{F}_{t_{0}}}(x_{2}).$

  5. 5.

We repeat $M$ times steps 1 - 4.

The formula in (36) is finally approximated as follows:

$\hat{P}_{t_{0}}=e^{-r\left(T-t_{0}\right)}\frac{1}{M}\sum_{j=1}^{M}\max\left\\{S_{1,t_{0},T}e^{\hat{x}_{1,j}}-S_{2,t_{0},T}e^{\hat{x}_{2,j}}-K;0\right\\}.$
(51)

The results are reported in Table 8 where all prices, obtained by applying
(37), belong to the Monte Carlo confidence interval determined with $10^{5}$
simulations. Prices obtained with (37) coincide with those obtained using the
double integration method. As expected, the main advantage of our approach is
the reduction of the computational time.

K | Single | Double | $CI^{-}_{95\%}$ | $CI^{+}_{95\%}$
---|---|---|---|---
0.0 | 6.149 | 6.149 | 6.127 | 6.182
| (1.12) | (34.66) | (51.22)
2.5 | 4.019 | 4.019 | 3.982 | 4.031
| (1.12) | (36.54) | (101.95)
5 | 2.266 | 2.266 | 2.251 | 2.293
| (1.20) | (34.32) | (103.62)
7.5 | 1.111 | 1.111 | 1.107 | 1.140
| (1.13) | (34.24) | (45.96)
10 | 0.533 | 0.532 | 0.526 | 0.551
| (1.13) | (34.07) | (44.94)
Table 8:  Comparison of 90 days European spread call option prices computed
using Formula (37) (Single), double integration method (Double) and Monte
Carlo 95% confidence interval. To determine the Monte Carlo upper and lower
bound we simulate $10^{5}$ realization of the spread payoff.

##  6 Conclusion

In this paper we proposed a new formula for pricing a spread option based on a
generic twice differentiable Copula function. The accuracy of the newly
introduced formula is investigated using both Monte Carlo simulations and a
numerical method for approximating multiple integrals. We showed that the
pricing of spread options based on a copula function is more efficient than
the alternative methods considered. The main advantage refers to the fact that
the pricing formula depends on a single integral that can be easily computed,
at least numerically, in any standard programming language.

## Acknowledgment

This work was partly supported by JST CREST Grant Number JPMJCR2115

## References

  * [1] P. Bjerksund and G. Stensland.  Closed form spread option valuation.  Department of Finance and Management Science, Norwegian School of Economics and Business Administration, Discussion Papers, 14, 12 2006.
  * [2] T. Bollerslev.  Generalized autoregressive conditional heteroskedasticity.  Journal of Econometrics, 31(3):307–327, 1986.
  * [3] P. P. Boyle.  Options: A monte carlo approach.  Journal of Financial Economics, 4(3):323–338, 1977.
  * [4] R. Caldana and G. Fusai.  A general closed-form spread option pricing formula.  Journal of Banking & Finance, 37(12):4893–4906, 2013.
  * [5] R. Carmona and V. Durrleman.  Pricing and hedging spread options.  SIAM Review, 45(4):627–685, 2003.
  * [6] S. Chiou and R. Tsay.  A copula-based approach to option pricing and risk assesment.  Journal of Data Science, 6:273–301, 01 2008.
  * [7] P. Christoffersen, S. Heston, and K. Jacobs.  Option valuation with conditional skewness.  Journal of Econometrics, 131(1-2):253–284, 2006.
  * [8] J. C. Cox.  Notes on option pricing i: Constant elasticity of diffusions.  Unpublished draft, Stanford University, 1975.
  * [9] J. C. Cox, J. E. Ingersoll, and S. A. Ross.  A theory of the term structure of interest rates.  Econometrica, 53(2):385–407, 1985.
  * [10] J. C. Duan.  The garch option pricing model.  Mathematical Finance, 5(1):13–32, 1995.
  * [11] R. F. Engle.  Autoregressive conditional heteroscedasticity with estimates of the variance of united kingdom inflation.  Econometrica, 50(4):987–1007, 1982.
  * [12] H. Herath, P. Kumar, and A. Amershi.  Crack spread option pricing with copulas.  Journal of Economics and Finance, 37:1–22, 01 2011.
  * [13] S. Heston.  A closed-form solution for options with stochastic volatility with applications to bond and currency options.  Review of Financial Studies, 6:327–343, 1993.
  * [14] S. L. Heston and S. Nandi.  A Closed-Form GARCH Option Valuation Model.  The Review of Financial Studies, 13(3):585–625, 06 2000.
  * [15] M. E. Johnson.  Multivariate Statistical Simulation.  Springer Berlin Heidelberg, 2011.
  * [16] E. Kirk.  Corrleation in energy markets.  Managing Energy Price Risk, pages 71–78, 1995.
  * [17] S. V. Malov.  On Finite-Dimensional Archimedean Copulas, pages 19–35.  2001\.
  * [18] K. V. Mardia.  Measures of multivariate skewness and kurtosis with applications.  Biometrika, 57(3):519–530, 12 1970.
  * [19] W. Margrabe.  The value of an option to exchange one asset for another.  The Journal of Finance, 33(1):177–186, 1978.
  * [20] L. Mercuri.  Option pricing in a garch model with tempered stable innovations.  Finance Research Letters, 5(3):172–182, 2008.
  * [21] F. J. Mhlanga and S. M. Kgomo.  On the sensitivity analysis of spread options using malliavin calculus.  2021\.
  * [22] R. B. Nelsen.  An introduction to copulas. Springer, 2006.
  * [23] R. L. Plackett.  A class of bivariate distributions.  Journal of the American Statistical Association, 60(310):516–522, 1965.
  * [24] M. Scarsini.  On measures of concordance.  Stochastica, 8(3):201–218, 1984.
  * [25] L. Schneider and B. Tavin.  Seasonal volatility in agricultural markets: modelling and empirical investigations.  Annals of Operations Research, 2021.
  * [26] T. Sønderby Christensen and F. E. Benth.  Modelling the joint behaviour of electricity prices in interconnected markets.  Quantitative Finance, 20(9):1441–1456, 2020.
  * [27] J. Van Belle, S. Vanduffel, and J. Yao.  Closed-form approximations for spread options in lévy markets.  Applied Stochastic Models in Business and Industry, 35(3):732–746, 2019.