A multivariate Lévy process model with linear correlation

In this paper, we develop a multivariate risk-neutral Lévy process model and discuss its applicability in the context of the volatility smile of multiple assets. Our formulation is based upon a linear combination of independent univariate Lévy processes and can easily be calibrated to a set of one-dimensional marginal distributions and a given linear correlation matrix. We derive conditions for our formulation and the associated calibration procedure to be well-defined and provide some examples associated with particular Lévy processes permitting a closed-form characteristic function. Numerical results of the option premiums on three currencies are presented to illustrate the effectiveness of our formulation with different linear correlation structures.


Introduction
Pure-jump Le´vy processes and their marginal infinite divisibility have attracted considerable attention amongst practitioners and academics for the prime reason that the flexibility of their distributions, such as heavy tails and asymmetry, suits very well various practical subjects in, for example, finance, telecommunications, turbulence, to mention just a few. In particular, in mathematical finance, there exists a vast literature on Le´vy process modeling for a single asset (see, for example, Carr et al. 2002, Kou 2002, Madan et al. 1998, Prause 1999, Rydberg 1997, and Schoutens and Teugels 1998. In the recent structured products market, it has become quite normal for the coupons to be determined by more than a single index, for example two FX rates USDJPY and AUDJPY. However, once the pricing model for multi-asset exotics is required to go beyond the Black-Scholes framework, there is still no market consensus on how to model both volatility smiles and an intended correlation structure. In the literature, there has recently been increasing interest in multivariate Le´vy process modeling. For example, the Le´vy copula models of Kallsen and Tankov (2006) and Tankov (2005) completely characterize the law of a multivariate Le´vy process, and are applied in pricing basket options, while Luciano and Schoutens (2006) and Moosbrucker (2006) propose the production of some dependence among components in the variance gamma process frameworks by setting a (fully or partially) common time-changing stochastic process for every component. A general construction of multi-factor Le´vy models from Le´vy models on rays is studied by Boyarchenko and Levendorskiı˘(2002). Among other examples of the multivariate model related to Le´vy processes and infinitely divisible distributions are the multidimensional structure model of discount factors of Bakshi et al. (2007), and the various credit derivatives models of Baxter (2007), Hull and White (2004), O'Kane and Schloegl (2003) and Kalemanova et al. (2007).
In this paper we propose a multivariate risk-neutral Le´vy process model based on a linear combination of independent univariate Le´vy processes that can easily be calibrated to a set of one-dimensional marginal distributions and a given linear correlation matrix. Our formulation has three contributing features. First, applying Le´vy processes permitting a parametric form of the characteristic function, we can easily derive the conditions for both our formulation and the Carr-Madan method (Carr and Madan 1999) to be well-defined (Proposition 3.1). Second, our model is built so as to precisely satisfy a given linear correlation matrix. Regardless of the increasing interest in the tail dependence in extreme value theory, from the viewpoint of financial practitioners it is usually more important to make sure of attaining an intended linear correlation. Although, as claimed in the literature, the linear correlation may be too simple, it is in fact sufficiently complicated in the sense that its change is unhedgable, even in the Black-Scholes framework. Third, and perhaps most importantly, compared with existing models, our formulation provides a considerably easier way to simulate the resultant multivariate Le´vy process as a linear combination of independent univariate Le´vy processes.
The rest of the paper is organized as follows. Section 2 gives some general notation and motivates our model development. Section 3 formulates the model and derives conditions for the entire calibration procedure to be well-defined. We also provide some explicit examples associated with particular Le´vy processes permitting closed-form characteristic functions. Section 4 presents numerical results for the option premiums of three currencies, USDJPY, EURJPY and AUDJPY, with different linear correlation structures and also considers the extreme value dependence. Finally, section 5 concludes.

Preliminaries
Let us begin with some notation and definitions that will be used in the following. R d is the d-dimensional Euclidean space with norm kÁk, R d 0 :¼ R d n f0g, and BðR d 0 Þ is the Borel -field of R d 0 . L(X) denotes the law of the random vector X, while ¼ L is the equality in law.
We denote by G 2 R dÂd the linear correlation matrix for a d-dimensional random vector of the form which is symmetric positive definite and where for each i, j, j i,j j 1. A stochastic process {X t : t ! 0} in R d is called a Le´vy process if it has stationary and independent increments, if it is stochastically continuous, and if X 0 ¼ 0 almost surely. It is well-known that the marginal distributions of a Le´vy process are infinitely divisible, that is, by the Le´vy-Khintchine representation, their characteristic function is uniquely characterized as where 1 2 R d , A 2 R dÂd is symmetric non-negative definite, and is a Le´vy measure, that is, a -finite measure on R d 0 , satisfying R R d 0 ðkzk 2^1 ÞðdzÞ 5 þ1. In this paper, we will only consider Le´vy processes with a finite second moment. The above characteristic function can then be rewritten as where the Le´vy measure is now in L 2 ðR d 0 Þ, that is, R R d 0 kzk 2 ðdzÞ 5 þ1. It also holds that every infinitely divisible distribution with a finite second moment is uniquely characterized by the three parameters ( 2 , A, ) in the representation (1). We denote by ' the distinguished logarithm of the characteristic function (1), that is, '( y) :¼ ln E[e ihy,X 1 i ]. For convenience, we will write ' X for the distinguished logarithm for the random vector X.
Consider a non-trivial d-dimensional random vector whose law is known componentwise and such that E[kZk 2 ]5þ1. Define also the following notation: (i) the mean vector : The main scope of our model is to ease the calibration procedure in finding a set of (C, X), where C 2 R dÂd and X :¼ (X 1 , . . . , X d ) 0 is an infinitely divisible random vector in R d with independent components, in the following sense: (i) the law of CX approximates that of Z componentwise; and (ii) the correlation matrix of CX is approximately or precisely G, which we give as an exogenous factor.
To motivate our discussion, let us begin with the simplest, yet very illustrative, example. Suppose the random vector Z is purely Gaussian. Let X be a d-dimensional standard normal random vector (with independent components) and let C be the lowertriangular matrix obtained via the Cholesky decomposition of the variance-covariance matrix diag()G diag(). It is then evident that the law of CX þ is identical to (and not only an approximation of) that of Z, since the normal random vector is characterized by the mean vector and the variance-covariance matrix. In most cases, the matrix C is taken as the lower-triangular Cholesky matrix, since then CX does not require full matrix multiplication due to the lower-triangular structure of C. Note, however, that C can also be any matrix satisfying CC 0 ¼ A. For example, it can be a matrix obtained via singular value decomposition.
Let us now apply the same approach to the case when Z is infinitely divisible without a Gaussian component. For simplicity, we assume that the random vector Z is in R 2 and C 2 R 2Â2 is the lower-triangular Cholesky matrix of diag()G diag(), whose (i, j)-element is denoted by c i,j . For each component k, we denote by k the Le´vy measure of Z k , that is We wish to find a random vector X that gives The law of the first component can easily be recovered by setting the characteristic function Let T r be a transformation defined by (T r )(B) :¼ (r À1 B) for B 2 BðR d 0 Þ. For the second component, using the independence between X 1 and X 2 , we set where the division in (3) is well-defined since the characteristic function of infinitely divisible distributions is necessarily non-zero. A crucial drawback here is that the above measure might not be a Le´vy measure. More precisely, this is clearly in L 2 ðR d 0 Þ, while it might no longer be a nonnegative measure. If that is the case, then 2 is not necessarily a characteristic function of any distribution.
Example 2.1: The Meixner distribution of Schoutens and Teugels (1998) is infinitely divisible and its Le´vy measure is characterized by three parameters (a, b, d) in the form where a40, b 2 (À, ), and d40. Now, let Z 1 and Z 2 be identically distributed Meixner random variables whose common Le´vy measure is set with the parameters (a, b) ¼ (1, À2) and d ¼ 2(cos(b/2)/a) 2 . This formulation guarantees E[Z 1 ] ¼ 0 and Var(Z 1 ) ¼ 1. Following the aforementioned procedure, we find independent random variables X 1 and X 2 , for example, with the correlation matrix G ¼ 1 À0:6 À0:6 1 : A numerical experiment shows that the measure (4) falls below zero approximately for z40.27. Such a measure certainly cannot be a Le´vy measure.

Model formulation
Let us now proceed to our model formulation. To maintain the full generality of the discussion, we drop the assumption of the lower-triangular structure of the matrix C 2 R dÂd . Let {X t : t ! 0} be a Le´vy process in R d without a Gaussian component and with independent components. Moreover, we mean by {X t,k : t ! 0} the kth component of {X t : t ! 0}, and denote by k its Le´vy measure on R 0 . Due to the component independence, the Le´vy process can be characterized componentwise in full. We denote by ' k the distinguished logarithm of the kth component, that is, e ' k ( y) :¼ E[e iyX 1,k ]. Let us then put the following assumption in force: for some 40.
Next, let {Z t : t ! 0} be a stochastic process in R d , with its kth component written as {Z t,k : t ! 0} and defined by provided that, for all k and l, j' l ðÀic k,l Þj 5 þ1: ð7Þ By assumption (5), this formulation guarantees that the given linear correlation matrix G is attained, that is, for each t ! 0 and for k 1 , k 2 , Remark 1: The independence of the components implies that the Le´vy measure on R d 0 of {X t : t ! 0} is supported on the axes of R d and can be written as where 0 is the Dirac delta measure at zero. After the linear transformation by the matrix C, the Le´vy measure {a(c 1,k , . . . , c d,k ) 0 : a 2 R} k¼1,. . .,d . In particular, if the matrix C is that obtained via singular value decomposition, its support is the union of all the eigenspaces of the correlation matrix G.
Now let r be a non-negative (constant) risk-free rate and let {S t : t ! 0} be a stochastic process in R d , again with its kth component written as {S t,k : t ! 0} and defined by where S 0,k is a positive constant. It is immediate that, for each k, the discounted process of {S t,k : t ! 0} is a martingale with respect to the natural filtration generated by the Le´vy process {X t : t ! 0}, since, for each t40, where the last equality holds by the component independence of {X t : t ! 0}. Next, fix T40 and let b k be the characteristic function of the logarithm of the kth component of S T , that is When its closed form is available and when the maturity T is not very small, the vanilla call option price with strike K40, that is can be efficiently computed using the well-known method of Carr and Madan (1999), where is a positive constant satisfying Noting that, for each k, in order for both (7) and (10) to hold we need to have, for every k and l, The following rewrites conditions (7) and (11) as integrability conditions of the underlying Le´vy measure. Remark 2: In view of (6) and (8), an addition of the drift to each component {X t,k : t ! 0} cancels out. We may thus define the distinguished logarithm at our convenience, such as for some c40, whenever those integrals are well-defined.
Example 3.2 (Meixner process): Suppose that the lth component {X t,l : t ! 0} is a Meixner process with parameter (a, b, d), as described in example 2.1. Its variance is then given in closed form: while the distinguished logarithm is as simple as 'ð yÞ ¼ 2d ln cos b 2 À ln cosh ay À ib 2 ! : We can show that the condition (11) reduces to p À b a 4 ð þ 1Þ max 0, max k ðc k,l Þ , and p þ b a 4 À ð þ 1Þ min 0, min k ðc k,l Þ : Note that the above domain of (a, b) is not quite so simple for the calibration procedure, in the sense that the parameters a and b restrict each other.
Example 3.3 (variance gamma process): The variance gamma process of Madan et al. (1998) is a Le´vy process obtained by evaluating Brownian motion with drift at a random time given by a gamma process: where {W t : t ! 0} is a standard Brownian motion and where {Y t : t ! 0} is a gamma process whose marginal is characterized by E[e iyY t ] ¼ (1 À iy) Àt/ (see also Boyarchenko and Levendorskiı˘(2000)). Its Le´vy measure is given in the form where 40, 2 R, and 40. This also implies that the variance gamma process can be viewed as the difference of two independent gamma processes. Its variance is given in closed form: while the distinguished logarithm is as simple as 'ð yÞ ¼ À 1 ln 1 À iy þ 2 2 y 2 : We can show that the condition (11) holds when We observe that, as in example 3.2, the parameters , , and restrict each other in a complex manner.
Example 3.4 (CGMY process): The CGMY process proposed by Carr et al. (2002Carr et al. ( , 2003 is also a Le´vy process obtained by tempering the Le´vy measure of the non-Gaussian stable distribution, and its Le´vy measure is given in the form ðdzÞ ¼ C n e ÀGjzj jzj 1þY n 1ðz 5 0Þ þ C p e ÀMz z 1þY p 1ðz 4 0Þ ! dz, z 2 R 0 , where C n , C p , G, M40 and Y n , Y p 2 (À1, 2). Its variance is given in closed form: while the distinguished logarithm is as simple as provided that Y n 6 ¼ 1 and Y p 6 ¼ 1. We can show that condition (11)

Numerical illustration: multi-currency
We consider the 1-year currency options on the three foreign exchanges USDJPY, EURJPY, and AUDJPY, and model their mean-correcting forward forex dynamics by the following stochastic process: where r d and r f denote the continuously compounded domestic and foreign risk-free rates, respectively. In our case, 'domestic' means 'JPY', while 'foreign' corresponds to either one of 'USD', 'EUR', or 'AUD'. We assume for simplicity that both risk-free rates are deterministic. F 0 is the spot FX rate, while e À(r f Àr d )t F 0 (¼ E[F t ]) is the so-called forward FX rate at time t. For notational convenience, we henceforth replace the time indices of {F t : t ! 0} and {Z t : t ! 0} by FX indices, that is, F 1 , F 2 and F 3 denote the USDJPY, EURJPY and AUDJPY rates at 1 year, respectively, while Z 1 , Z 2 and Z 3 are defined analogously.

Underlying Le´vy process
For our numerical experiments, we employ a Le´vy process whose lth component has a Le´vy measure given by l ðdzÞ ¼ 1 l jzj ½e À l jzj 1ðz 5 0Þ þ e À l z 1ðz 4 0Þ dz, z 2 R 0 , where l 40, l 40, and l 40. We set the characteristic function as From its gamma structure, it immediately follows that e ' X t,l ð yÞ ¼ 1 þ iy l Àt= l 1 À iy l Àt= l As can be seen, our Le´vy process is very similar to the variance gamma process, and in fact it can even recover the variance gamma process by setting and where the parameters (, ) are as defined in example 3.3. Here, by unfastening the above dependence between and , we not only make the Le´vy measure more flexible, but also simplify the conditions for (11), compared with those derived in example 3.3, to l 4 À ð þ 1Þ minð0, min k ðc k,l ÞÞ, that is, each condition consists of a single parameter. Yet our Le´vy process exhibits the features of the variance gamma process: (i) its characteristic function is given in closed form; and (ii) it can be viewed as the difference of two independent gamma processes. One thing our Le´vy process does not inherit from the variance gamma process is the property of 'time-changed Brownian motion' of (12), since parameters and are no longer dependent as in (15) and (16). From a practical point of view, however, this property is not very important, as long as our Le´vy process can still be simulated by two independent gamma processes. For condition (5) to be satisfied, with the help of the closed-form variance (14), we fix the parameter for each l, Hence, in the calibration procedure, we will only control the remaining two parameters ( l , l ) for the lth component. Note that, under condition (18), the realized linear correlation matrix of the three FXs is assumed to be precisely G, whatever the calibration result.

Calibration with given correlation
We test for the following correlation structures: We also prepare G 4 , which is an artificial modification of G 1 , to test negative correlation. Here, we obtain the matrix C such that CC 0 ¼ G via singular value decomposition, and these are given by, respectively, We perform the calibration via the Nelder-Mead direct search method to minimize the difference between market premiums (market-pr) and model premiums (model-pr) in the sense of the root mean square error (rmse), defined as (rmse) :¼ 1 #(pr) X ½(market-pr) À (model-pr) 2 ! 1=2 : In our example, #(pr) ¼ 15 (that is, five premiums for each of the three currencies) and the calibration is a minimization problem with six parameters ( 1 , 1 , . . . , 3 , 3 ), due to the common variance condition (5), where those parameters can take values in the domain (17). During the calibration procedure, we compute the model premiums using the Carr-Madan method (9). It is well-known that the calibration of infinitely divisible marginals is an ill-posed inverse problem (see, for example, Cont and Tankov (2006)), and therefore local minima are always a problem. We thus perform the minimization problem many times, with different common variance 2 in (5), and from different randomly chosen starting points. The calibration results in table 4 are the best in each case. As can be seen from figure 1, our formulation is capable of fitting market premiums very well with different linear correlation structures.
Scatter plots of 500 points for currency pairs are given in figure 2. We observe that, compared with the ordinary Gaussian framework, our Le´vy process model can express a dependence in more-negative values. These realizations can easily be generated based on a linear combination of three random variables, X 1 , X 2 and X 3 . In particular, thanks to the gamma structure of the underlying Le´vy measure (13), each component X k is simulated as the difference of independent gamma random variables, for which a sample generation routine is included with most standard mathematical and statistical software. This simplicity in implementation should be a great advantage compared with existing multivariate Le´vy process models. For example, the stochastic volatility formulation of Luciano and Schoutens (2006) is based on a time-changed Brownian motion, and usually requires a Euler-type discretization of sample paths for simulation, while the Le´vy copula model (Tankov 2005) resorts to the series representation of Le´vy processes, which is wellknown to be computationally very expensive for most simulations.

Extreme value dependence
Our formulation can be used to approximate a set of onedimensional marginal distributions with a given linear correlation structure, and is not modeled to have control over the extreme value dependence. It is worth checking how our multivariate model behaves with extreme values. To this end, for a fixed T, and for each k, let F À1 k ðuÞ be the inverse cumulative distribution function of Z T,k , that is,  F À1 k ðuÞ :¼ inffv 2 R : PðZ T,k vÞ 4 ug. Then, define the measures of the extreme value dependence, We perform Monte Carlo simulation to estimate the above quantities based on 1e þ 5 independent replications, and report them in table 5. We observe quite direct results of a greater extreme value dependence for a greater correlation. Our multivariate model also seems to be capable of expressing strongly correlated downside jumps, which has often been claimed to be important in the economics literature. We can also see that the negative correlation induces a negative extreme value dependence.
Remark 1: The limiting values lim u#0 LL (u) and lim u#0 UU (u) are often called, respectively, the lower and upper tail-dependence coefficients, and their estimation methods and applications to financial markets have been studied intensively. For details, see, for example, Caillault and Gue´gan (2005), Frahm et al. (2005), Malevergne and Sornette (2005), and Poon et al. (2003).
Remark 2: We have chosen the matrix C to be the one obtained via singular value decomposition. However, as mentioned earlier, it could be any matrix satisfying CC 0 ¼ G, and it is natural to be tempted to choose instead the lower-triangular Cholesky matrix. For example, for the correlation matrix G 1 , the matrix C 1 is then Unfortunately, with this lower triangular matrix C 1 , we could not obtain a satisfactory result via a calibration problem with six parameters. The lower triangular structure of C 1 may reduce the simultaneous calibration to a bootstrap procedure, that is, the first component is calibrated with only 1X 1 , then the second with 0.9682X 2 , where 0.25X 1 is fixed, and so on. As can be seen, however, as this bootstrap procedure treats later components, the quality of the calibration seems to become poorer. Perhaps this is partially because our underlying Le´vy measure (13) of three parameters is not sufficiently flexible for this purpose. Indeed, our experiment indicates that the CGMY process with six parameters (Example 3.4) fits all three components very well. It is difficult, however, to justify the fact that the flexibility of the Le´vy measure always provides significant contributions to the calibration quality in the bootstrap procedure. Another concern is that the procedure can certainly be carried out in arbitrary order in which different components are calibrated. It is not clear to what extent this component ordering affects the calibration quality.

Concluding remarks
In this paper we have formulated a multivariate riskneutral Le´vy process model that can easily be calibrated to a set of one-dimensional marginal distributions and a given linear correlation matrix. Our model is built on a linear combination of independent univariate Le´vy processes so as to precisely realize the given linear correlation matrix. Applying Le´vy processes permitting closed-form characteristic functions, our formulation provides a very simple calibration procedure in the Carr-Madan framework. The resultant multivariate Le´vy processes can easily be generated for simulation use as a linear combination of independent univariate Le´vy processes. Numerical results indicate that our model is capable of simultaneously fitting three currency options, with different linear correlation structures, and strongly correlated downside jumps can also be well expressed.
As future research, it would be interesting to conduct further empirical analyses on different sets of currencies and different types of financial assets, for example a set of equity indices, and with different Le´vy processes, such as the CGMY process. Moreover, an extension to the multimaturity context should certainly be pursued. This may pave the way to the evaluation of popular, but very intricate, financial exotics in the Le´vy process framework, a subject that will be studied by Kawai (2009).