Multivariate Option Pricing with Pair-Copulas

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 decompositionmethod 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.


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: where  > 0 is the maturity time of the option,   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, ( 1 (), . . .,  4 ()) is the payoff of the option written on the four assets  1 , . . .,  4 , whose prices, at any time  ∈ [0, ], are   () for  = 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  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: where  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 parameters  = 1 and  = 1, namely, a GARCH(, ) 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,  1 , . . .,  4 and a discrete set of times whose elements represent a day between day 0, beginning of transactions, and day  which can be taken as the end of our daily observations of the asset's prices.The quantity   (), for  ∈ {0, . . ., } and  ∈ {1, . . ., 4}, stands for the closing price of the underlying  at the trading day .The one-day log-return of the th asset is defined as (, ) ∈ {1, . . ., 4} × {0, . . .,  − 1} . ( Note that the objective 4-variate distribution of the logreturns is specified conditional to all return information available at time , that is conditional to the sigma algebra F  := (( 1, , . . .,  4, ) :  ≤ ).
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 a GARCH(1, 1) process with Gaussian innovations; see, for example, [4]; namely, where, for every  = 1, . . ., 4,   is the expected daily logreturn of the asset   , while  ,+1 is the Gaussian innovation under the objective, or real world, probability P linked to the return  ,+1 .In particular,  ,+1 , conditioned to all returns information available at time , that is conditioned to F  , has mean zero and variance  2 , which evolves as a GARCH(1, 1) processes with parameters   ,   ,   > 0 such that   +   < 1.
Since the marginal distributions, at a given time , are specified conditional to a common set of information, that is, with respect to the -algebra F  , 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 riskneutral 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 probability Q satisfies a local risk-neutral valuation relationship (LRNVR) (see [4, Def.2.1, Th.2.2]), then the law of the returns under Q is given by where  * ,+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 assets  1 , . . .,  4 , under the risk-neutral measure Q, by a copulabased approach.We would like to underline that in the  dimensional case for  > 2 and differently from the bivariate one, the main difficulty is that of finding the right -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 paircopula representation in view of the analysis of the threedimensional case; see [9].Definition 1 (tree).A tree  = {, } is an acyclic graph, where  is its set of nodes and  is its set of edges (unordered pairs of nodes).
Definition 2 (regular vine).A regular vine tree on  variables is a structure of connected trees  1 , . . .,  −1 , with nodes   and edges   ,  = 1, . . .,  − 1, such that the following conditions are satisfied: (1)  1 is a tree with a set of nodes  1 = {1, . . ., } and a set of edges denoted by  1 ; (2) for  = 2, . . .,  − 1 the tree   has node set   =  −1 and a set of edges denoted by   ; (3) two edges in tree   are joined in tree  +1 if they share a common node in tree   .
Here we focus on two special cases of regular vines, namely, the D-vine and the canonical 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: let  denote the joint density probability function of three random variables  1 ,  2 ,  3 , and then  can be decomposed as follows: where, for every ,  = 1, . . ., 3,   stands for the marginal probability distribution function of   ,   is the corresponding marginal density function, and   is the associated bivariate copula density, while  13|2 represents the bivariate copula density of  1 ,  3 conditional to the second component.
In order to reproduce a density function through a paircopula decomposition in dimension  > 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 paircopulas, 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 -Student, the Clayton and the Gumbel copulas since they are the most used types in financial applications; see, for example, [1,[5][6][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][12][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 density  13|2 in (6).Let us consider the following definition of bivariate conditional distribution function in terms of copula: where  and V are uniform random variables and Θ is the set of parameters characterizing the copula of the joint distribution of  and V.
We would like to underline that the function ℎ plays a key role in the pair-copula decomposition approach in a dimensional setting when  > 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 ℎ 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.

Calibration of the Model.
The pair-copula decomposition and the ℎ-function described in Section 1.1 can be calibrated starting from the definition of the Gaussian innovations  , 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 a -dimensional setting,  ≥ 2, we define the standardized innovations as follows: for all   ∈ R,  = 1, . . ., .Moreover, we assume that the copula  in ( 9) is a parametric copula of parameter  P .
Latter assumption is justified since the dependence structure existing between the variables   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.
Hereafter, we will work using the random vector ( 1, , . . .,  4, ), with  ∈ [0, ], where  > 0 is the maturity date of our investment.The next step is that of finding the copula on the joint distribution of the random vector  1 , . . .,  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  1 , . . .,  4 under the objective measure P.

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  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 assets  1 , . . .,  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 , for such a particular vine structure.
In order to simulate the assets  1 , . . .,  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   ,  = 1, . . ., 4, according to the algorithm described in [10].
Concerning the option price estimation, let us recall that at time  the option price is given by the following equation: We first estimate the price, defined as    , inserting the prices   (), observed in the market at time , for  = 1, . . ., 4, into the option payoff equation.Then we use a Monte Carlo approach to estimate the price, defined as pmc  (Θ 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, min where {  }  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.

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 of closing levels for each index.On these indexes we write two option contracts,  1 and  2 , respectively, which are defined by the following payoffs at maturity  > 0: respectively, where  > 0 is the strike price of the option  1 and  = 1, . . ., 4. Note that  = 1000 in our numerical implementation.In Figure 1 we plot the four 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.
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.
Moreover, 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 Figure 3, strengthen previous deduction.
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 the GARCH(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.
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 Table 3.Using the data collected in Table 3, 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.Each edge of the D-vines tree corresponds to a paircopula density; therefore, in our implementation, the paircopula decomposition of the joint density reads as follows: Given the pair-copula decomposition, we further investigate the dependence structure constructing the bivariate scatter plots of the copula data   , for  = 1, . . ., 4, where   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.
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, the -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:  , ,  , ,  , ,  ,| ,  ,| , and  ,|, .
The first three copulas are estimated using the copula data of the vector  and comparing the AIC values for each possible choice of the copula function.The results are as follows: (i)  , is a Clayton copula; (ii)  , is a Clayton copula too; (iii)  , is a Gumbel copula.
Once the copula functions have been estimated, we can compute the corresponding parameters, which are reported in Table 4.
From the second tree of Figure 4 it turns out that we have to compute  ,| and  ,| .Since they are conditional pair copula we have to infer first the proper ℎ-function.Then, we can estimate the conditional copula using the sample created by the ℎ-functions.In a four-variable Dvine structure we have to compute three ℎ-functions, since  The last pair copula is a Gumbel copula function too.Note that for both  ,| and  ,|, we have θ = 1, which implies that the variables in the distribution function are conditionally independent so that, for the copula  ,| 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  ,|, 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:

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  Q  .In particular, starting from the corresponding  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   ().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 . 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   /(1−  −  ), and we exploit the standardized innovations  Q  given by the following identity: 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:

Figure 1 :
Figure 1: Plot of the data sets.

Figure 2 :
Figure 2: Plot of the indexes returns.

𝑆 2 𝑇𝜂
, =  , exp [( − )   − 1 , ] , 2, , . . .,  , ) [15,e  , is the standard deviation of the underlying  at time .By the well-known results for the GARCH model (see, for example, [14, Part V, Sec.16]) the  , are .., standard Gaussian random variables, even if the stochastic processes  1 , ...,   are not independent.By the Sklar's Theorem we have that the joint distribution  (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  1 , . . .,   are continuous stochastic processes, there exists an unique copula  such that

Table 2 :
Estimated coefficients and corresponding standard errors.

Table 4 :
Copula parameters estimation.(  |   ), (  |   ) and (  |   ).Hence we simulate (  |   ) using the Clayton ℎ-function ℎ(  ,   ), while for (  |   ) we use the Gumbel ℎ-function ℎ(  ,   ), and finally the corresponding ℎfunction for (  |   ) is again the Clayton one.Once the conditional distributions have been simulated, we can estimate the best-fitting copula functions for the pairs ((  | ), (  |   )) and ((  |   ), (  |   )), and we have that the selected copula is for both of the pairs the Gumbel copula, with parameters as in Table4.Finally, the copula  ,|, has to be computed.Again we have to use the ℎ-function to simulate a sample for the conditional distribution functions (  |   ,   ) and (  |   ,   ).