# Multivariate Option Pricing with Pair-Copulas

Authors: Anna Barban; Luca Di Persio
Journal: Journal of Probability
 (2014)

---

## Abstract

We propose a copula-based approach to solve the option pricing problem in the risk-neutral setting and with respect to a structured derivative written on several underlying assets. Our analysis generalizes similar results already present in the literature but limited to the trivariate case. The main difficulty of such a generalization consists in selecting the appropriate vine structure which turns to be of D-vine type, contrary to what happens in the trivariate setting where the canonical vine is sufficient. We first define the general procedure for multivariate options and then we will give a concrete example for the case of an option written on four indexes of stocks, namely, the S&P 500 Index, the Nasdaq 100 Index, the Nasdaq Composite Index, and the Nyse Composite Index. Moreover, we calibrate the proposed model, also providing a comparison analysis between real prices and simulated data to show the goodness of obtained estimates. We underline that our pair-copula decomposition method produces excellent numerical results, without restrictive assumptions on the assets dynamics or on their dependence structure, so that our copula-based approach can be used to model heterogeneous dependence structure existing between market assets of interest in a rigorous and effective way.


---

## Body

## 1. Introduction

In what follows we will consider a European option written on 4 assets. We will assume that the risk-neutral setting holds true, according to the framework defined in [1]; that is, the price of each of the four considered underlying assets only depends on its history; moreover, it is independent form the others past behaviour. Previous condition allows us to write the price of the above mentioned option as the following discounted value:
(1)
p
t
=
e
-
r
f
T
-
t
E
Q
G
(
S
1
(
T
)
,
…
,
S
4
(
T
)
)
∣
F
t
,
where T
>
0 is the maturity time of the option, r
f is a real positive parameter usually representing the interest rate given by a bank for our cash deposit and it is assumed to be constant over the whole option’s life, G
(
S
1
(
T
)
,
…
,
S
4
(
T
)
) is the payoff of the option written on the four assets S
1
,
…
,
S
4, whose prices, at any time t
∈
0
,
T, are S
i
(
t
) for i
=
1
,
…
,
4. We would like to underline that (1) is given directly under the risk-neutral (martingale) measure Q and it represents the (fair or no-arbitrage) price of the option with payoff G determined by the so-called martingale approach, see, for example, [2, Section  5.4].Under suitable assumptions, the price determined in (1) can be rewritten according to the expected value definition as follows:
(2)
p
t
=
e
-
r
f
T
-
t
∬
0
+
∞
G
(
s
1
,
…
,
s
4
)
f
Q
(
s
1
,
…
,
s
4
)
d
s
1
⋯
d
s
4
,
where f
Q
(
·
) is the joint density probability function of the underlying assets with respect to the risk-neutral probability measure Q.Our aim is to apply a pair-copula decomposition approach to reproduce the joint density of the 4 underlying assets as correlated random variables. Latter goal requires first to find a method able to (statistically) describe the behaviour of each underlying asset. In particular, we assume that the underlying assets returns evolve as generalized autoregressive conditional heteroskedastic processes with parametersp
=
1 and q
=
1, namely, a GARCH
(
p
,
q
) process. Latter choice is motivated by the fact that the GARCH processes properly describe the evolution of variables that do not have constant volatility over time; see, for example, [3]. Therefore, GARCH processes turn out to be very effective models for heteroskedastic processes, as in the case of financial time series which exhibit structured interrelations.In particular, in what follows we consider 4 assets, namely,S
1
,
…
,
S
4 and a discrete set of times whose elements represent a day between day 0, beginning of transactions, and day T which can be taken as the end of our daily observations of the asset’s prices. The quantity S
i
(
t
), for t
∈
0
,
…
,
T and i
∈
1
,
…
,
4, stands for theclosing price of the underlying i at the trading day t. The one-day log-return of the ith asset is defined as
(3)
r
i
,
t
+
1
:
=
log
⁡
S
i
(
t
+
1
)
S
i
(
t
)
,
i
,
t
∈
1
,
…
,
4
×
0
,
…
,
T
-
1
.
Note that the objective 4-variate distribution of the log-returns is specified conditional to all return information available at time t, that is conditional to the sigma algebra F
t
:
=
σ
(
(
r
1
,
s
,
…
,
r
4
,
s
)
:
s
≤
t
).In order to derive the joint risk-neutral return process from the objective one, we assume that the objective marginal distributions of the assets returns evolve as aGARCH
(
1,1
) process with Gaussian innovations; see, for example, [4]; namely,
(4)
r
i
,
t
+
1
=
μ
i
+
η
i
,
t
+
1
,
σ
i
,
t
+
1
2
=
γ
i
+
α
i
η
i
,
t
+
1
2
+
β
i
σ
i
,
t
2
,
η
i
,
t
+
1
∣
F
t
~
P
N
0
,
σ
i
,
t
2
,
where, for every i
=
1
,
…
,
4, μ
i is the expected daily log-return of the asset S
i, while η
i
,
t
+
1 is the Gaussian innovation under the objective, orreal world, probability P linked to the return r
i
,
t
+
1. In particular, η
i
,
t
+
1, conditioned to all returns information available at time t, that is conditioned to F
t, has mean zero and variance σ
i
,
t
2 which evolves as a GARCH
(
1,1
) processes with parameters γ
i
,
α
i
,
β
i
>
0 such that α
i
+
β
i
<
1.Since the marginal distributions, at a given timet, are specified conditional to a common set of information, that is, with respect to the σ-algebra F
t, then we are allowed to exploit copula theory techniques to derive the joint conditional distribution. In what follows, inspired by [5, 6], see also [1, 7], we assume that the objective and the risk-neutral copulas are the same.The main idea behind our option pricing model is that we can use a convenient map to transform each marginal objective distribution to its risk-neutral counterpart, as in [4]. Then we define a proper 4-dimensional copula function, obtained by the pair-copula construction method (see, for example, [8]) instead of deriving the joint risk-neutral distribution directly. Finally, we determine the fair price of the option by taking the discounted expected value of the option’s payoff under the risk-neutral measure.In particular, assuming that the risk-neutral probabilityQ satisfies alocal risk-neutral valuation relationship (LRNVR) (see [4, Def.2.1, Th.2.2]), then the law of the returns under Q is given by
(5)
r
i
,
t
+
1
=
r
f
-
1
2
σ
i
,
t
2
+
η
i
,
t
+
1
*
,
σ
i
,
t
+
1
2
=
γ
i
+
α
i
r
i
,
t
+
1
-
μ
i
2
+
β
i
σ
i
,
t
2
,
η
i
,
t
+
1
*
∣
F
t
~
Q
N
0
,
σ
i
,
t
2
,
where η
i
,
t
+
1
* are the Gaussian innovations under Q. Exploiting the transform defined in (5) it is possible to model the marginal behaviour of the marginal log-returns under the measure Q.


### 1.1. Pair-Copula Decomposition

In this section, we show how to model the interdependence of the underlying assetsS
1
,
…
,
S
4, under the risk-neutral measure Q, by a copula-based approach. We would like to underline that in the d dimensional case for d
>
2 and differently from the bivariate one, the main difficulty is that of finding the right d-dimensional copula which properly reproduces the dependence structure existing between the single pairs of components. We solve such a problem by a multivariate copula construction method, namely, the pair-copula decomposition approach; see, for example, [8].Let us start introducing some basic concepts for the pair-copula representation in view of the analysis of the three-dimensional case; see [9].Definition 1 (tree).
A treeT
=
N
,
E is an acyclic graph, where N is its set of nodes and E is its set of edges (unordered pairs of nodes).Definition 2 (regular vine).
A regular vine tree ond variables is a structure of connected trees T
1
,
…
,
T
d
-
1, with nodes N
i and edges E
i, i
=
1
,
…
,
d
-
1, such that the following conditions are satisfied: (1)
T
1 is a tree with a set of nodes N
1
=
1
,
…
,
d and a set of edges denoted by E
1;
(2)
fori
=
2
,
…
,
d
-
1 the tree T
i has node set N
i
=
E
i
-
1 and a set of edges denoted by E
i;
(3)
two edges in treeT
i are joined in tree T
i
+
1 if they share a common node in tree T
i.Here we focus on two special cases of regular vines, namely, theD-vine and thecanonical vine; see, for example, [8] for details.Before going into details, we would like to underline that while from the pair-copula decomposition point of view the three-dimensional case can be easily treated (for example, the D-vine and the canonical vine specification coincide), difficulties arise when the number of correlated dimensions increases.In three dimensions the construction method proceeds as follows: letf denote the joint density probability function of three random variables X
1
,
X
2
,
X
3, and then f can be decomposed as follows:
(6)
f
x
1
,
x
2
,
x
3
=
f
1
x
1
f
2
x
2
f
3
x
3
c
12
F
1
x
1
,
F
2
x
2
×
c
23
F
2
x
2
,
F
3
x
3
c
13
∣
2
×
F
1
∣
2
x
1
∣
x
2
,
F
3
∣
2
x
3
∣
x
2
,
where, for every i
,
j
=
1
,
…
,
3, F
i stands for the marginal probability distribution function of X
i, f
i is the corresponding marginal density function, and c
i
j is the associated bivariate copula density, while c
13
∣
2 represents the bivariate copula density of X
1
,
X
3 conditional to the second component.In order to reproduce a density function through a pair-copula decomposition in dimensiond
>
3, we have first to select an appropriate decomposition structure. In particular, from a financial point of view, previous statement means that we first have to analyze the economic relationship between the variables of interest in order to find the variable that governs the interactions in the data set. If such a variable exists, then a canonical vine specification should be the best choice; otherwise, we will use the D-vine decomposition structure; see, for example, [8].Once the vine structure has been chosen, we are left with the selection of the order of the variables, the families of pair-copulas, and the specification of their parameters.The order of the variables is set according to the following rule: first we compute the dependence measure Kendall’s tau for each pair of variables and we order the variables in the first tree on the basis of their dependence structure; this means that the two variables with the strongest dependence are put in the first two nodes of the first tree and so on. The next trees are determined as a consequence; see [10]. Then we select the best-fitting pair-copula family; namely, we have to select the pair-copula specification that guarantees the best fitting with observed data. Such a problem can lead to different solutions depending on different setting. In our specific case, we are justified to consider only the Gaussian, the t-Student, the Clayton and the Gumbel copulas since they are the most used types in financial applications; see, for example, [1, 5–7]. The particular copula is selected exploiting standard information criteria such as Akaike information criterion (AIC) and the Bayesian information criterion (BIC) respectively; see, for example, [11–13] and references therein, for each edge in each tree. In particular, the choice is made in order to minimize the AIC and the BIC coefficient, respectively. Finally, the selected copula’s parameters are estimated using the maximum likelihood criterion for each pair.Once we have obtained the results for the first tree, we need to construct a sample for the conditional bivariate distribution, in order to find the pair-copula associated with the copula densityc
13
∣
2 in (6). Let us consider the following definition of bivariate conditional distribution function in terms of copula:
(7)
h
u
,
v
,
Θ
:
=
F
u
∣
v
=
∂
C
u
,
v
u
,
v
,
Θ
∂
v
,
where u and v are uniform random variables and Θ is the set of parameters characterizing the copula of the joint distribution of u and v.We would like to underline that the functionh plays a key role in the pair-copula decomposition approach in a d-dimensional setting when d
>
3, since it allows to reproduce the conditional behaviour of a random vector in terms of a bivariate function. Indeed, by an iterative algorithm, we can rewrite the conditional distribution using a proper choice of h and of the copula written on its marginals.Once the pair-copula method has been theoretically implemented, thus obtaining the joint distribution of interest, we are left with the need to calibrate the model in order to obtain satisfactory numerical results. The calibration procedure will be the main goal of the next section.


### 1.2. Calibration of the Model

The pair-copula decomposition and theh-function described in Section 1.1 can be calibrated starting from the definition of the Gaussian innovations η
i
,
t of the model, namely, a GARCH
(
1,1
), we have chosen for the marginal distributions of the underlying assets.In particular, if we are in ad-dimensional setting, d
≥
2, we define the standardized innovations as follows:
(8)
Z
1
,
s
,
Z
2
,
s
,
…
,
Z
d
,
s
s
:
=
η
1
,
s
σ
1
,
s
,
η
2
,
s
σ
2
,
s
,
…
,
η
d
,
s
σ
d
,
s
,
∀
(
i
,
s
)
∈
1
,
…
,
d
×
0
,
T
,
where σ
i
,
s is the standard deviation of the underlying i at time s. By the well-known results for the GARCH model (see, for example, [14, Part V, Sec.16]) the Z
i
,
s are i
.
i
.
d, standard Gaussian random variables, even if the stochastic processes Z
1
,
…
,
Z
d are not independent. By the Sklar’s Theorem we have that the joint distribution F (under the objective measure P) of the underlying innovations can be written in terms of its marginals; see, for example, [15, Sec.1]. In particular, since Z
1
,
…
,
Z
d are continuous stochastic processes, there exists an unique copula C such that
(9)
F
z
1
,
…
,
z
d
=
C
F
1
z
1
,
…
,
F
d
z
d
,
for all z
i
∈
R, i
=
1
,
…
,
d. Moreover, we assume that the copula C in (9) is a parametric copula of parameter θ
P. Latter assumption is justified since the dependence structure existing between the variables z
i is usually given by a copula function that depends on a particular parameter. Given the standardized innovations, we want to infer their dependence structure under P using the pair-copula decomposition method. Then, we will use the dependence structure, namely, the copula function that better reproduces the joint behaviour of the underlyings, to price an option written on four indexes which will be considered as underlyings.In order to use the pair-copula decomposition, and in particular theh-function defined in (7), we compute the corresponding estimated standardized innovations in the interval (
0,1
), applying the cumulative distribution function Φ of the standard normal distribution N
(
0,1
). Thus, the corresponding uniform variables U
i, for i
=
1
,
…
,
4 are defined as follows:
(10)
U
1
,
s
,
…
,
U
4
,
s
:
=
Φ
η
^
1
,
s
σ
^
1
,
s
,
…
,
Φ
η
^
4
,
s
σ
^
4
,
s
.
Hereafter, we will work using the random vector (
U
1
,
s
,
…
,
U
4
,
s
), with s
∈
0
,
T, where T
>
0 is the maturity date of ourinvestment.The next step is that of finding the copula on the joint distribution of the random vectorU
1
,
…
,
U
4. In particular, we have to choose the bivariate copulas that best fit our data. Such a copula structure will be then exploited to reproduce the 4-variate copula density and thus the joint distribution density, for the random vector U
1
,
…
,
U
4 under the objective measure P.


### 1.3. Option Pricing

In this section, we will consider (2) to solve the associated option pricing problem exploiting results obtained through previous sections. To reproduce the joint density function f
Q we use, as for the density under the objective probability P, a multivariate copula constructed through the pair-copula approach. In particular, we assume that the copula under Q belongs to the same family as the one determined under P, possibly being characterized by different parameters.The option pricing problem will be solved in two steps. First we simulate the underlying assetsS
1
,
…
,
S
4 under the risk-neutral measure Q; then, we simulate the option price and we calibrate the model to determine the set Θ
Q of the copulas’ parameters under Q. The set Θ
Q is composed of all the parameters of the pair-copulas used to create the proper vine specification. For example, in the 3-dimensional case, we have Θ
Q
=
(
θ
12
Q
,
θ
13
Q
,
θ
23
∣
1
Q
), for such a particular vine structure.In order to simulate the assetsS
1
,
…
,
S
4 under Q we consider the risk-neutral transformation of the returns defined in (5); see [4] for details. This turns out to be a recursive method to estimate both the returns and the volatility terms under the risk-neutral measure. Moreover, in order to maintain the dependence structure obtained from the market data, we simulate, using the vine specification previously defined, the variables U
i, i
=
1
,
…
,
4, according to the algorithm described in [10].Concerning the option price estimation, let us recall that at timetthe option price is given by the following equation:
(11)
p
t
=
e
-
r
f
T
-
t
E
Q
G
S
1
T
,
…
,
S
4
T
∣
F
t
.We first estimate the price, defined asp
t
M, inserting the prices S
i
(
t
), observed in the market at time t, for i
=
1
,
…
,
4, into the option payoff equation. Then we use a Monte Carlo approach to estimate the price, defined as p
^
t
mc
(
Θ
Q
), by (1) and using the observation simulated from the vine structure. Finally, we calculate the set Θ
Q minimizing the sum of the quadratic error; namely,
(12)
min
⁡
Θ
Q
∑
i
=
1
n
‍
p
^
t
i
mc
-
p
t
i
M
2
,
where t
i
i are the past observations of the option prices. The question about the choice of the parameter set Θ
Q corresponds to the calibration of the pricing model.


## 1.1. Pair-Copula Decomposition

In this section, we show how to model the interdependence of the underlying assetsS
1
,
…
,
S
4, under the risk-neutral measure Q, by a copula-based approach. We would like to underline that in the d dimensional case for d
>
2 and differently from the bivariate one, the main difficulty is that of finding the right d-dimensional copula which properly reproduces the dependence structure existing between the single pairs of components. We solve such a problem by a multivariate copula construction method, namely, the pair-copula decomposition approach; see, for example, [8].Let us start introducing some basic concepts for the pair-copula representation in view of the analysis of the three-dimensional case; see [9].Definition 1 (tree).
A treeT
=
N
,
E is an acyclic graph, where N is its set of nodes and E is its set of edges (unordered pairs of nodes).Definition 2 (regular vine).
A regular vine tree ond variables is a structure of connected trees T
1
,
…
,
T
d
-
1, with nodes N
i and edges E
i, i
=
1
,
…
,
d
-
1, such that the following conditions are satisfied: (1)
T
1 is a tree with a set of nodes N
1
=
1
,
…
,
d and a set of edges denoted by E
1;
(2)
fori
=
2
,
…
,
d
-
1 the tree T
i has node set N
i
=
E
i
-
1 and a set of edges denoted by E
i;
(3)
two edges in treeT
i are joined in tree T
i
+
1 if they share a common node in tree T
i.Here we focus on two special cases of regular vines, namely, theD-vine and thecanonical vine; see, for example, [8] for details.Before going into details, we would like to underline that while from the pair-copula decomposition point of view the three-dimensional case can be easily treated (for example, the D-vine and the canonical vine specification coincide), difficulties arise when the number of correlated dimensions increases.In three dimensions the construction method proceeds as follows: letf denote the joint density probability function of three random variables X
1
,
X
2
,
X
3, and then f can be decomposed as follows:
(6)
f
x
1
,
x
2
,
x
3
=
f
1
x
1
f
2
x
2
f
3
x
3
c
12
F
1
x
1
,
F
2
x
2
×
c
23
F
2
x
2
,
F
3
x
3
c
13
∣
2
×
F
1
∣
2
x
1
∣
x
2
,
F
3
∣
2
x
3
∣
x
2
,
where, for every i
,
j
=
1
,
…
,
3, F
i stands for the marginal probability distribution function of X
i, f
i is the corresponding marginal density function, and c
i
j is the associated bivariate copula density, while c
13
∣
2 represents the bivariate copula density of X
1
,
X
3 conditional to the second component.In order to reproduce a density function through a pair-copula decomposition in dimensiond
>
3, we have first to select an appropriate decomposition structure. In particular, from a financial point of view, previous statement means that we first have to analyze the economic relationship between the variables of interest in order to find the variable that governs the interactions in the data set. If such a variable exists, then a canonical vine specification should be the best choice; otherwise, we will use the D-vine decomposition structure; see, for example, [8].Once the vine structure has been chosen, we are left with the selection of the order of the variables, the families of pair-copulas, and the specification of their parameters.The order of the variables is set according to the following rule: first we compute the dependence measure Kendall’s tau for each pair of variables and we order the variables in the first tree on the basis of their dependence structure; this means that the two variables with the strongest dependence are put in the first two nodes of the first tree and so on. The next trees are determined as a consequence; see [10]. Then we select the best-fitting pair-copula family; namely, we have to select the pair-copula specification that guarantees the best fitting with observed data. Such a problem can lead to different solutions depending on different setting. In our specific case, we are justified to consider only the Gaussian, the t-Student, the Clayton and the Gumbel copulas since they are the most used types in financial applications; see, for example, [1, 5–7]. The particular copula is selected exploiting standard information criteria such as Akaike information criterion (AIC) and the Bayesian information criterion (BIC) respectively; see, for example, [11–13] and references therein, for each edge in each tree. In particular, the choice is made in order to minimize the AIC and the BIC coefficient, respectively. Finally, the selected copula’s parameters are estimated using the maximum likelihood criterion for each pair.Once we have obtained the results for the first tree, we need to construct a sample for the conditional bivariate distribution, in order to find the pair-copula associated with the copula densityc
13
∣
2 in (6). Let us consider the following definition of bivariate conditional distribution function in terms of copula:
(7)
h
u
,
v
,
Θ
:
=
F
u
∣
v
=
∂
C
u
,
v
u
,
v
,
Θ
∂
v
,
where u and v are uniform random variables and Θ is the set of parameters characterizing the copula of the joint distribution of u and v.We would like to underline that the functionh plays a key role in the pair-copula decomposition approach in a d-dimensional setting when d
>
3, since it allows to reproduce the conditional behaviour of a random vector in terms of a bivariate function. Indeed, by an iterative algorithm, we can rewrite the conditional distribution using a proper choice of h and of the copula written on its marginals.Once the pair-copula method has been theoretically implemented, thus obtaining the joint distribution of interest, we are left with the need to calibrate the model in order to obtain satisfactory numerical results. The calibration procedure will be the main goal of the next section.


## 1.2. Calibration of the Model

The pair-copula decomposition and theh-function described in Section 1.1 can be calibrated starting from the definition of the Gaussian innovations η
i
,
t of the model, namely, a GARCH
(
1,1
), we have chosen for the marginal distributions of the underlying assets.In particular, if we are in ad-dimensional setting, d
≥
2, we define the standardized innovations as follows:
(8)
Z
1
,
s
,
Z
2
,
s
,
…
,
Z
d
,
s
s
:
=
η
1
,
s
σ
1
,
s
,
η
2
,
s
σ
2
,
s
,
…
,
η
d
,
s
σ
d
,
s
,
∀
(
i
,
s
)
∈
1
,
…
,
d
×
0
,
T
,
where σ
i
,
s is the standard deviation of the underlying i at time s. By the well-known results for the GARCH model (see, for example, [14, Part V, Sec.16]) the Z
i
,
s are i
.
i
.
d, standard Gaussian random variables, even if the stochastic processes Z
1
,
…
,
Z
d are not independent. By the Sklar’s Theorem we have that the joint distribution F (under the objective measure P) of the underlying innovations can be written in terms of its marginals; see, for example, [15, Sec.1]. In particular, since Z
1
,
…
,
Z
d are continuous stochastic processes, there exists an unique copula C such that
(9)
F
z
1
,
…
,
z
d
=
C
F
1
z
1
,
…
,
F
d
z
d
,
for all z
i
∈
R, i
=
1
,
…
,
d. Moreover, we assume that the copula C in (9) is a parametric copula of parameter θ
P. Latter assumption is justified since the dependence structure existing between the variables z
i is usually given by a copula function that depends on a particular parameter. Given the standardized innovations, we want to infer their dependence structure under P using the pair-copula decomposition method. Then, we will use the dependence structure, namely, the copula function that better reproduces the joint behaviour of the underlyings, to price an option written on four indexes which will be considered as underlyings.In order to use the pair-copula decomposition, and in particular theh-function defined in (7), we compute the corresponding estimated standardized innovations in the interval (
0,1
), applying the cumulative distribution function Φ of the standard normal distribution N
(
0,1
). Thus, the corresponding uniform variables U
i, for i
=
1
,
…
,
4 are defined as follows:
(10)
U
1
,
s
,
…
,
U
4
,
s
:
=
Φ
η
^
1
,
s
σ
^
1
,
s
,
…
,
Φ
η
^
4
,
s
σ
^
4
,
s
.
Hereafter, we will work using the random vector (
U
1
,
s
,
…
,
U
4
,
s
), with s
∈
0
,
T, where T
>
0 is the maturity date of ourinvestment.The next step is that of finding the copula on the joint distribution of the random vectorU
1
,
…
,
U
4. In particular, we have to choose the bivariate copulas that best fit our data. Such a copula structure will be then exploited to reproduce the 4-variate copula density and thus the joint distribution density, for the random vector U
1
,
…
,
U
4 under the objective measure P.


## 1.3. Option Pricing

In this section, we will consider (2) to solve the associated option pricing problem exploiting results obtained through previous sections. To reproduce the joint density function f
Q we use, as for the density under the objective probability P, a multivariate copula constructed through the pair-copula approach. In particular, we assume that the copula under Q belongs to the same family as the one determined under P, possibly being characterized by different parameters.The option pricing problem will be solved in two steps. First we simulate the underlying assetsS
1
,
…
,
S
4 under the risk-neutral measure Q; then, we simulate the option price and we calibrate the model to determine the set Θ
Q of the copulas’ parameters under Q. The set Θ
Q is composed of all the parameters of the pair-copulas used to create the proper vine specification. For example, in the 3-dimensional case, we have Θ
Q
=
(
θ
12
Q
,
θ
13
Q
,
θ
23
∣
1
Q
), for such a particular vine structure.In order to simulate the assetsS
1
,
…
,
S
4 under Q we consider the risk-neutral transformation of the returns defined in (5); see [4] for details. This turns out to be a recursive method to estimate both the returns and the volatility terms under the risk-neutral measure. Moreover, in order to maintain the dependence structure obtained from the market data, we simulate, using the vine specification previously defined, the variables U
i, i
=
1
,
…
,
4, according to the algorithm described in [10].Concerning the option price estimation, let us recall that at timetthe option price is given by the following equation:
(11)
p
t
=
e
-
r
f
T
-
t
E
Q
G
S
1
T
,
…
,
S
4
T
∣
F
t
.We first estimate the price, defined asp
t
M, inserting the prices S
i
(
t
), observed in the market at time t, for i
=
1
,
…
,
4, into the option payoff equation. Then we use a Monte Carlo approach to estimate the price, defined as p
^
t
mc
(
Θ
Q
), by (1) and using the observation simulated from the vine structure. Finally, we calculate the set Θ
Q minimizing the sum of the quadratic error; namely,
(12)
min
⁡
Θ
Q
∑
i
=
1
n
‍
p
^
t
i
mc
-
p
t
i
M
2
,
where t
i
i are the past observations of the option prices. The question about the choice of the parameter set Θ
Q corresponds to the calibration of the pricing model.


## 2. Numerical Implementation

In our implementation, we work with a dataset that comes from the U.S market. In particular, we consider the following indexes: the Standard and Poor’s (S&P) 500 index, the Nasdaq 100 index, the Nasdaq Composite index, and the New York Stock Exchange (Nyse) Composite index; see, for example, the related Wikipedia occurrences for detailed definitions of these market indexes. We have considered such indexes between January 1, 2012 and December 28, 2012, for a total of 249 days resulting in the same amount ofclosing levels for each index. On these indexes we write two option contracts, X
1 and X
2, respectively, which are defined by the following payoffs at maturity T
>
0:
(13)
X
1
T
:
=
max
⁡
max
⁡
S
i
T
-
K
,
0
,
(14)
X
2
T
=
max
⁡
S
i
T
-
min
⁡
S
i
T
,
respectively, where K
>
0 is the strike price of the option X
1 and i
=
1
,
…
,
4. Note that K
=
1000 in our numerical implementation. In Figure 1 we plot the four data sets.Figure 1
Plot of the data sets.Notice that the plots of the closing prices of the four indexes are not comparable, since they do not fluctuate around a common mean value. Hereafter, we will consider the log-return of each index instead of time series of the prices, as in (3). By an analysis of the time series of the returns, in Figure 2, we can infer an heteroskedasticity characteristic for the considered setting; hence, a GARCH model turns out to be considered as a natural choice.Figure 2
Plot of the indexes returns.The GARCH parameters estimation is the first step of our implementation. In particular, once the time series of the underlying assets prices are given, we compute the returns estimation. Before analyzing the results obtaining using the GARCH model, let us underline some statistics for each return series together with the response given by the Jarque-Bera normality test; see, for example, [16], for each series. From Table 1 we see that the Jarque-Bera statistics are quite high, especially for the S&P 500 series; thus we can infer that the returns are not normally distributed.Table 1
Descriptive statistics.
NYSE
Nasdaq100
Nasdaq Com
S&P500
Mean
0.350
0.466
0.448
0.377
Standard deviation
0.0085
0.0096
0.0093
0.0079
Kurtosis
3.549
3.685
3.629
3.915
JB statistics
3.18
5.07
4.09
8.68Moreover, both the kurtosis coefficients (greater than 3 which is the characteristic kurtosis value for a Gaussian distribution) and the quantile-quantile plot of the returns, as in Figure3, strengthen previous deduction.Figure 3
Quantile-quantile plot of the indexes returns.We would also like to underline the presence of fat tail for the S&P 500 time series of returns, in agreement with the (non-Gaussianity) results given by the Jarque-Bera statistics. We point out that using theGARCH
(
1,1
) model we obtain a highly satisfactory description of our market behaviour; indeed all the processes have a log-likelihood value rather high, between 800 and 850. Moreover, it results that the most significant variable in the model is the one corresponding to the variance past value, with the latter being in agreement with the use of a heteroskedastic model. In Table 2 we present the estimated coefficients and the corresponding estimation errors of the coefficients of our GARCH
(
1,1
) process for each index.Table 2
Estimated coefficients and corresponding standard errors.
NYSE
Nasdaq100
Nasdaq Com
S&P500
μ
^
5.03
e
-
04
9.38
e
-
04
7.47
e
-
04
5.29
e
-
04
(
5.43
e
-
04
)
(
6.03
e
-
04
)
(
5.98
e
-
04
)
(
5.13
e
-
04
)
γ
^
3.64
e
-
06
2.93
e
-
06
3.19
e
-
06
3.02
e
-
06
(
4.53
e
-
06
)
(
3.57
e
-
06
)
(
4.24
e
-
06
)
(
3.60
e
-
06
)
α
^
0.0471
0.0488
0.0419
0.0415
(
0.0372
)
(
0.0315
)
(
0.0324
)
(
0.0318
)
β
^
0.9035
0.9220
0.9229
0.9127
(
0.0899
)
(
0.0592
)
(
0.0710
)
(
0.0796
)In order to construct a proper pair-copula structure for the multivariate density function, we first have to find the vine structure that best fits our data, and then we will properly order the variables in the chosen vine. We first investigate the dependence structure between pairs out of {Nyse, Nasdaq 100, Nasdaq Com, S&P 500}, and the result, obtained using the sample Kendall’s tau, is given in Table3.Table 3
Sample Kendall’s tau.
NYSE
Nasdaq100
Nasdaq Com
S&P500
NYSE
1
Nasdaq100
0.5224
1
Nasdaq Com
0.6018
0.8985
1
S&P500
0.7286
0.7320
0.7548
1Using the data collected in Table3, we can infer that there is not a variable that governs the interactions. Thus we will consider a D-vine structure for the pair-copula decomposition. Moreover, we can see that the strongest dependence appears to be, as expected, between the Nasdaq 100 Index and the Nasdaq Composite Index, followed by the one between the S&P 500 Index and the Nasdaq Composite Index. The D-vine structure for our data set is given in Figure 4.Figure 4
D-vine structure with four variables.Each edge of the D-vines tree corresponds to a pair-copula density; therefore, in our implementation, the pair-copula decomposition of the joint density reads as follows:(15)
f
(
x
N
,
x
N
c
,
x
S
,
x
N
y
)
=
f
N
(
x
N
)
·
f
N
c
(
x
N
c
)
·
f
S
(
x
S
)
·
f
N
y
(
x
N
y
)
·
c
N
,
N
c
F
N
(
x
N
)
,
F
N
c
(
x
N
c
)
·
c
N
c
,
S
F
N
c
(
x
N
c
)
,
F
S
(
x
S
)
·
c
S
,
N
y
F
S
(
x
S
)
,
F
N
y
(
x
N
y
)
·
c
N
,
S
∣
N
c
F
N
∣
N
c
x
N
∣
x
N
c
,
F
S
∣
N
c
x
S
∣
x
N
c
·
c
N
c
,
N
y
∣
S
F
N
c
∣
S
x
N
c
∣
x
S
,
F
N
y
∣
S
x
N
y
∣
x
S
·
c
N
,
N
y
∣
N
c
,
S
F
N
∣
N
c
,
S
x
N
∣
x
N
c
,
x
S
,
·
c
N
,
N
y
N
c
,
S
i
i
i
·
c
N
,
N
F
N
y
∣
N
c
,
S
x
N
y
∣
x
N
c
,
x
S
·
c
N
,
N
y
N
c
,
S
i
i
i
·
c
N
,
N
F
N
y
∣
N
c
,
S
x
N
y
∣
x
N
c
,
x
S
.Given the pair-copula decomposition, we further investigate the dependence structure constructing the bivariate scatter plots of the copula dataU
i, for i
=
1
,
…
,
4, where U
i is the uniform vector in (10). Graphs in Figures 5 and 6 show the existence of tail dependence between almost all the pairs. In particular, we can easily notice that there is significant dependence between different time series concerning their tails. The latter fact is rather obvious since all the used financial indexes belong to the USA market, and hence, even if their behaviour is also influenced by non-USA assets, a quite strong link between them is exactly what we have to take into account. We would also like to underline that we may expect similar results, namely, the existence of tail dependence between almost all the pairs, even in the case of pairs constituted by indexes belonging to different economies, for example, DAX and NYSE, at least when rather global financial events occur, as in the case of the worldwide economic crisis of 2007/2008. It is worth to mention that, in oder to study such type of dependencies, one can make use of the so-called upper/lower dependence indexes, which can be considered as indicators of the dependence values in the tails of the bivariate distributions.Figure 5
Scatter plots of the copula data sets.Figure 6
Scatter plots of the copula data sets.Previous consideration is the first step in the copula choosing process. The next step, actually the most important one, is that of comparing the AIC values for the four possible copula alternatives, that is, the Gaussian, theT-Student, the Clayton, and the Gumbel copula. Let us first analyze the bivariate copulas of the first tree, and then we will estimate the corresponding copula functions’ parameters. From Figure 4 we infer that we have to compute the following copulas: C
N
,
N
c, C
N
c
,
S, C
S
,
N
y, C
N
,
S
∣
N
c, C
N
c
,
N
y
∣
S, and C
N
,
N
y
∣
N
c
,
S.The first three copulas are estimated using the copula data of the vectorU and comparing the AIC values for each possible choice of the copula function. The results are as follows:(i)
C
N
,
N
c is a Clayton copula;
(ii)
C
N
c
,
S is a Clayton copula too;
(iii)
C
S
,
N
y is a Gumbel copula.Once the copula functions have been estimated, we can compute the corresponding parameters, which are reported in Table4.Table 4
Copula parameters estimation.
θ
^
C
N
,
N
c
9.2084
C
N
c
,
S
3.5953
C
S
,
N
y
7.2419
C
N
,
S
∣
N
c
1.3708
C
N
c
,
N
y
∣
S
1.00
C
N
,
N
y
∣
N
c
,
S
1.00From the second tree of Figure4 it turns out that we have to compute C
N
,
S
∣
N
c and C
N
c
,
N
y
∣
S. Since they are conditional pair copula we have to infer first the proper h-function. Then, we can estimate the conditional copula using the sample created by the h-functions. In a four-variable D-vine structure we have to compute three h-functions, since we need to estimate the following conditional distribution functions: F
(
x
N
c
∣
x
S
), F
(
x
N
y
∣
x
S
) and F
(
x
N
∣
x
N
c
). Hence we simulate F
(
x
N
c
∣
x
S
) using the Clayton h-function h
C
(
x
N
c
,
x
S
), while for F
(
x
N
y
∣
x
S
) we use the Gumbel h-function h
G
u
(
x
N
y
,
x
S
), and finally the corresponding h-function for F
(
x
N
∣
x
N
c
) is again the Clayton one. Once the conditional distributions have been simulated, we can estimate the best-fitting copula functions for the pairs F
x
N
∣
x
N
c
,
F
(
x
S
∣
x
N
c
) and F
x
N
c
∣
x
S
,
F
x
N
y
∣
x
S, and we have that the selected copula is for both of the pairs the Gumbel copula, with parameters as in Table 4. Finally, the copula C
N
,
N
y
∣
N
c
,
S has to be computed. Again we have to use the h-function to simulate a sample for the conditional distribution functions F
x
N
∣
x
N
c
,
x
S and F
x
N
y
∣
x
N
c
,
x
S. The last pair copula is a Gumbel copula function too. Note that for both C
N
c
,
N
y
∣
S and C
N
,
N
y
∣
N
c
,
S we have θ
^
=
1, which implies that the variables in the distribution function are conditionally independent so that, for the copula C
N
c
,
N
y
∣
S and given the S&P Index, the Nasdaq Composite and the Nasdaq 100 are independent. Similarly, given the Nasdaq Composite and the S&P, the Nasdaq 100 and the NYSE are conditionally independent. Moreover, we have that the density copula function c
N
,
N
y
∣
N
c
,
S is equal to 1; see [10]; therefore, such copula components do not play a role in our pair-copula decomposition structure.Exploiting previous considerations about the estimated parameters, we obtain that the joint density function can be decomposed as follows:(16)
f
(
x
N
,
x
N
c
,
x
S
,
x
N
y
)
=
f
N
(
x
N
)
·
f
N
c
(
x
N
c
)
·
f
S
(
x
S
)
·
f
N
y
(
x
N
y
)
·
c
N
,
N
c
Clayton
F
N
(
x
N
)
,
F
N
c
(
x
N
c
)
·
c
N
c
,
S
Clayton
F
N
c
(
x
N
c
)
,
F
S
(
x
S
)
·
c
S
,
N
y
Gumbel
F
S
(
x
S
)
,
F
N
y
(
x
N
y
)
·
c
N
,
S
∣
N
c
Gumbel
F
N
∣
N
c
x
N
∣
x
N
c
,
F
S
∣
N
c
x
S
∣
x
N
c
.


### 2.1. Option Pricing with Four Underlying Assets

In this section we analyze how the option pricing problem can be solved using our pair-copula decomposition.We first consider the option pricing problem for an option written on four underlyings whit payoff given as in (13). Then we consider the same problem, but for an option with payoff as in (14). Following what we have developed in the previous sections from the theoretical model, we first compute the option price using the market prices at the final time and then we compare it with the prices under the risk-neutral probability Q using the same pair-copula decomposition, but with proper parameters. To find the optimal parameters set we use a sort of confidence interval for the possible values θ
i
Q. In particular, starting from the corresponding θ
i
P computed under the objective probability measure P, we look for the values that belong to a predetermined interval such that they minimize the difference between the market price and the estimated one, as in (12). Hereafter, we will face the problem of finding the copula parameters which allow obtaining an option price, computed under the risk-neutral measure Q
~
P, as close as possible to the market option price computed exploiting the observed data. We consider first an option with payoff given by (13) and we compute its market price at the end of our investment period, say p
M
(
T
). Then we compute the option price under the risk-neutral probability measure Q, and related numerical results (prices) are shown in Table 6. Before computing the option price under the risk-neutral measure Q, we have to derive the assets prices under the probability Q. Let us note that, by assumption, the dependence structure under Q is the same as the one under the objective probability P. The first step in the method used to construct the proper prices processes is that of simulating a 4-dimensional, uniformly distributed vector in [
0,1
]
4, say U
Q
=
[
U
1
Q
,
U
2
Q
,
U
3
Q
,
U
4
Q
]. Such a simulation is iterated up to 248 times, in order to create a sample of data as long as the data set of observations. Then we construct the returns series according to the specification given in (5). To do this we use the coefficients computed with the GARCH
(
1,1
) model in Table 2, we set an initial volatility level equal to γ
i
/
1
-
α
i
-
β
i, and we exploit the standardized innovations Z
i
Q given by the following identity:
(17)
Z
i
,
t
Q
=
Φ
-
1
U
i
,
t
Q
.
The next step consists in recursively computing the returns series under Q. This procedure allows us to obtain the fact (see [4]) that the final asset price under the risk-neutral probability is defined as follows:
(18)
S
i
,
T
=
S
i
,
t
exp
⁡
(
T
-
t
)
r
f
-
1
2
∑
s
=
t
+
1
T
‍
σ
i
,
s
2
+
∑
s
=
t
+
1
T
‍
η
i
,
s
,
where T is the final observation time, r
f is the constant risk-free rate, and σ
i
,
t
2 is the volatility of the asset i at time t, computed using the second line of (5).The final step consists in an option price estimation made using both Monte Carlo method, with 10000 simulated trajectories for each underlying asset, and the sample data series previously constructed. The estimate procedure is then completed finding the copula parameters that minimize the differences between the market option price at timeT, that is, p
M
(
T
), and the Monte Carlo simulated price, say p
MC
(
T
). In Table 5, we present the estimated copula parameters that minimize the difference between p
M and p
MC both for an option with payoff as in (13), say θ
^
1, and for an option with payoff as in (14), say θ
^
2.Table 5
Copula parameters underQ.
θ
^
1
θ
^
2
C
N
,
N
c
10.2084
10.2084
C
N
c
,
S
3.5953
2.5953
C
S
,
N
y
7.2419
6.2419
C
N
,
S
∣
N
c
1.3708
0.8708
C
N
c
,
N
y
∣
S
1.00
1.00
C
N
,
N
y
∣
N
c
,
S
1.00
1.00Table 6
Market price and Monte Carlo price.
Option1
Option2
Market pricep
M (T)
7.32
e
+
003
6.91
e
+
003
Monte Carlo estimatep
MC (T)
7.27
e
+
003
6.67
e
+
003Table6 reports the prices computed for an option with payoff as in (13), say Option 1, respectively, for a second option with payoff as in (14), say Option 2.Using data collected in Table6, we can infer that the price computed using the Monte Carlo simulation is rather similar to the market price. This means that the copula-based option pricing is a good estimation method to find the option price. The difference between the two prices can be seen as the witnessing of the lack of arbitrage free hypotheses. An alternative explanation relies on the fact that even if Q
~
P we do not have that the copula decomposition structure remains the same under these two probability measures.We would also like to underline that, from a purely computational point of view, the simulation of a sample from the D-vine dependence structure in four dimensions is rather more complicated than in three dimensions because of the increased algorithm complexity. Previous results show that also for the four-dimensional option pricing problem the copula approach gives satisfactory results. Latter approach is quite general and it can be adapted to a larger set of assets distribution. In fact we can use different specification of the copula functions and we can exploit the flexibility given by the pair-copula decomposition.


## 2.1. Option Pricing with Four Underlying Assets

In this section we analyze how the option pricing problem can be solved using our pair-copula decomposition.We first consider the option pricing problem for an option written on four underlyings whit payoff given as in (13). Then we consider the same problem, but for an option with payoff as in (14). Following what we have developed in the previous sections from the theoretical model, we first compute the option price using the market prices at the final time and then we compare it with the prices under the risk-neutral probability Q using the same pair-copula decomposition, but with proper parameters. To find the optimal parameters set we use a sort of confidence interval for the possible values θ
i
Q. In particular, starting from the corresponding θ
i
P computed under the objective probability measure P, we look for the values that belong to a predetermined interval such that they minimize the difference between the market price and the estimated one, as in (12). Hereafter, we will face the problem of finding the copula parameters which allow obtaining an option price, computed under the risk-neutral measure Q
~
P, as close as possible to the market option price computed exploiting the observed data. We consider first an option with payoff given by (13) and we compute its market price at the end of our investment period, say p
M
(
T
). Then we compute the option price under the risk-neutral probability measure Q, and related numerical results (prices) are shown in Table 6. Before computing the option price under the risk-neutral measure Q, we have to derive the assets prices under the probability Q. Let us note that, by assumption, the dependence structure under Q is the same as the one under the objective probability P. The first step in the method used to construct the proper prices processes is that of simulating a 4-dimensional, uniformly distributed vector in [
0,1
]
4, say U
Q
=
[
U
1
Q
,
U
2
Q
,
U
3
Q
,
U
4
Q
]. Such a simulation is iterated up to 248 times, in order to create a sample of data as long as the data set of observations. Then we construct the returns series according to the specification given in (5). To do this we use the coefficients computed with the GARCH
(
1,1
) model in Table 2, we set an initial volatility level equal to γ
i
/
1
-
α
i
-
β
i, and we exploit the standardized innovations Z
i
Q given by the following identity:
(17)
Z
i
,
t
Q
=
Φ
-
1
U
i
,
t
Q
.
The next step consists in recursively computing the returns series under Q. This procedure allows us to obtain the fact (see [4]) that the final asset price under the risk-neutral probability is defined as follows:
(18)
S
i
,
T
=
S
i
,
t
exp
⁡
(
T
-
t
)
r
f
-
1
2
∑
s
=
t
+
1
T
‍
σ
i
,
s
2
+
∑
s
=
t
+
1
T
‍
η
i
,
s
,
where T is the final observation time, r
f is the constant risk-free rate, and σ
i
,
t
2 is the volatility of the asset i at time t, computed using the second line of (5).The final step consists in an option price estimation made using both Monte Carlo method, with 10000 simulated trajectories for each underlying asset, and the sample data series previously constructed. The estimate procedure is then completed finding the copula parameters that minimize the differences between the market option price at timeT, that is, p
M
(
T
), and the Monte Carlo simulated price, say p
MC
(
T
). In Table 5, we present the estimated copula parameters that minimize the difference between p
M and p
MC both for an option with payoff as in (13), say θ
^
1, and for an option with payoff as in (14), say θ
^
2.Table 5
Copula parameters underQ.
θ
^
1
θ
^
2
C
N
,
N
c
10.2084
10.2084
C
N
c
,
S
3.5953
2.5953
C
S
,
N
y
7.2419
6.2419
C
N
,
S
∣
N
c
1.3708
0.8708
C
N
c
,
N
y
∣
S
1.00
1.00
C
N
,
N
y
∣
N
c
,
S
1.00
1.00Table 6
Market price and Monte Carlo price.
Option1
Option2
Market pricep
M (T)
7.32
e
+
003
6.91
e
+
003
Monte Carlo estimatep
MC (T)
7.27
e
+
003
6.67
e
+
003Table6 reports the prices computed for an option with payoff as in (13), say Option 1, respectively, for a second option with payoff as in (14), say Option 2.Using data collected in Table6, we can infer that the price computed using the Monte Carlo simulation is rather similar to the market price. This means that the copula-based option pricing is a good estimation method to find the option price. The difference between the two prices can be seen as the witnessing of the lack of arbitrage free hypotheses. An alternative explanation relies on the fact that even if Q
~
P we do not have that the copula decomposition structure remains the same under these two probability measures.We would also like to underline that, from a purely computational point of view, the simulation of a sample from the D-vine dependence structure in four dimensions is rather more complicated than in three dimensions because of the increased algorithm complexity. Previous results show that also for the four-dimensional option pricing problem the copula approach gives satisfactory results. Latter approach is quite general and it can be adapted to a larger set of assets distribution. In fact we can use different specification of the copula functions and we can exploit the flexibility given by the pair-copula decomposition.


## 3. Conclusions

In this paper we have faced a multivariate option pricing problem considering options written on 4 underlying assets by copula decomposition techniques. In particular, inspired by the method proposed in [9, 10], we have strictly generalized similar results already present in literature, see, for example, [1, 5, 6], which are limited to the trivariate case.The main difficulty to overcome passing from the three-dimensional to the four-dimensional setting is that in the latter case we have to select the appropriate vine structure which turns to be of D-vine type, contrary to what happens in the trivariate setting where the canonical vine is sufficient.The latter means that a deeper analysis of the relationship between the 4 underlying assets has to be done in order to give correct vine specification.We have also shown, by accurate simulations, that the particular pair-copula decomposition method which we proposed to solve option pricing problem produces excellent numerical results, without any particular restrictive assumption on the assets behaviour or on the dependence structure existing between them.Moreover, a copula-based approach has the necessary flexibility needed to face the most different scenarios characterizing present financial markets, and indeed such a technique allows us to model almost every dependence structure existing between the market assets in a precise and effective way.



---

*Source: 839204-2014-12-25.xml*