Calibration for Weak Variance-Alpha-Gamma Processes
Boris Buchmann Kevin W. Lu Dilip B. Madan
stat.ME q-fin.MF 60G51, 62F10, 60E10
	









	
											
												
												
													












𝒦
𝒱
𝒮
𝒥








 








































































































































equationsection





theoremTheorem
lemmaLemma


remark0Remark
remark

propProposition
corollaryCorollary
*reverseReverse direction
definitionDefinition
propositionProposition
theoremsection
propositionsection
lemmasection
corollarysection
remark0section
definitionsection







Calibration for Weak Variance-Alpha-Gamma Processes[This
research was partially supported by ARC grant DP160104037.]
    Boris BuchmannResearch School of Finance, Actuarial Studies & Statistics,
Australian National University,
ACT 0200,
Australia.
Email: boris.buchmann@anu.edu.au

Kevin W. LuMathematical Sciences Institute,
	Australian National University,
	ACT 0200,
	Australia.
	Email: kevin.lu@anu.edu.au
Dilip B. MadanRobert H. Smith School of Business, University of Maryland, College Park, MD. 20742, USA, Email: dbm@rhsmith.umd.edu

    December 30, 2023
========================================================================================================================================================================================================================================================================================================================================================================================================================================




The weak variance-alpha-gamma process is a multivariate Lévy process constructed by weakly subordinating Brownian motion, possibly with correlated components with an alpha-gamma subordinator. It generalises the variance-alpha-gamma process of Semeraro constructed by traditional subordination. We compare three calibration methods for the weak variance-alpha-gamma process, method of moments, maximum likelihood estimation (MLE) and digital moment estimation (DME). We derive a condition for Fourier invertibility needed to apply MLE and show in our simulations that MLE produces a better fit when this condition holds, while DME produces a better fit when it is violated. We also find that the weak variance-alpha-gamma process exhibits a wider range of dependence and produces a significantly better fit than the variance-alpha-gamma process on an S&P500-FTSE100 data set, and that DME produces the best fit in this situation.



2000 MSC Subject Classifications:  Primary: 60G51

 Secondary: 62F10, 60E10

Keywords: Brownian Motion, Gamma Process, Lévy Process, Subordination, Variance-Gamma, Variance-Alpha-Gamma, Self-Decomposability, Log-Return, Method of Moments, Maximum Likelihood Estimation, Digital Moment Estimation.





§ INTRODUCTION



The subordination of Brownian motion has important applications in mathematical finance, and acts as a time change that models the flow of information, measuring time in volume of trade, as opposed to real time. This idea was initiated by Madan and Seneta in <cit.> who introduced the variance-gamma (VG) process for modelling stock prices, where the subordinate is Brownian motion and the subordinator is a gamma process.

Subordination can be applied to model dependence in multivariate price processes. The multivariate VG process in <cit.> uses n-dimensional Brownian motion as its subordinate and a univariate gamma process as its subordinator, which gives it a restrictive dependence structure, where components cannot have idiosyncratic time changes and must have equal kurtosis when there is no skewness. Models based on linear combinations of independent Lévy processes <cit.> also do not account for both common and idiosyncratic time changes. These deficiencies are addressed by the use of an alpha-gamma subordinator, resulting in the variance-alpha-gamma (VAG) process which was introduced by Semeraro in <cit.> and also studied in <cit.>. However, in this case, the Brownian motion subordinate must have independent components, which also restricts the dependence structure.


To be precise, let B=(B_1,…,B_n), where T=(T_1,…,T_n) be independent n-dimensi­onal processes, where B is Brownian motion and T is a subordinator. Subordination is the operation that produces the process B∘ T defined by (B∘ T)(t):=(B_1(T_1(t)),…,B_n(T_n(t))),t≥ 0. Subordination in the case when T has indistinguishable components has been studied in <cit.>, and when B has independent components in <cit.>. In these cases, which we refer to as traditional subordination, B∘ T is a Lévy process, otherwise it may not be (see <cit.>, their Proposition 3.9). We refer the reader to <cit.> for a thorough discussion of traditional subordination and its applications.


In <cit.>, we introduced the weak subordination of B and T, an operation that extends traditional subordination and always produces a Lévy process B⊙ T. Then the weak variance-alpha-gamma (WVAG) process can be constructed using weak subordination instead of traditional subordination, while allowing for the Brownian motion to have possibly correlated components. The WVAG process exhibits a wider range of dependence while remaining parsimoniously parametrised, each component has both common and idiosyncratic time changes, it has VG marginals with independent levels of kurtosis, and the jump measure has full support.

Weak subordination also has applied in quantitative finance. In <cit.>, various marginal consistent dependence models have been constructed by weak subordination. In <cit.>, log return modelling based on the WVAG process was applied in instantaneous portfolio theory. In <cit.>, weak subordination using subordinators with arbitrary marginal components and dependence specified by a Lévy copula was studied in the context of financial information flows.

Maximum likelihood estimation (MLE) has been used to fit financial data to a univariate VG process in <cit.>, to a bivariate VG process in <cit.>, to a WVAG process in <cit.>, and to a factor-based subordinated Brownian motion in <cit.>, a generalisation of the WVAG process. Since the density function of the VAG and WVAG distribution is not explicitly known but its characteristic function is, the density function is computed using Fourier inversion.

In this paper, we derive a sufficient condition in terms of the parameters for Fourier invertibility, a problem that to our knowledge is not addressed in the existing literature. Then we compare MLE with method of moments (MOM) and digital moment estimation (DME) from <cit.>. Using simulations we find that MLE produces a better fit when the Fourier invertibility condition is satisfied but that DME is better when it is violated. In addition, we fit both the WVAG and VAG model to an S&P500-FTSE100 data set and show that the weak model has a significantly better fit, and that DME is the better method in this situation. Finally, using a condition for the self-decomposability of the WVAG process from <cit.>, we find that the log returns are self-decomposable.

This paper is structured as follows. In Section 2, we review the definition and properties of the WVAG process, and other preliminaries. In Section 3, we derive a condition for Fourier invertibility. In Section 4, we apply MOM, MLE, DME to simulated and real data, and discuss our findings. In Section 5, we conclude the paper.






§ WEAK VARIANCE-ALPHA-GAMMA PROCESS

Let ^n be n-dimensional Euclidean space whose elements are row vectors =(x_1,…,x_n) with canonical basis {_k:1≤k≤n}.

Let =' denote the Euclidean product, ^2=' denote the Euclidean norm, and let ^2_Σ:=Σ'.
For n-dimensional processes X and Y, X Y indicates that X and Y are identical in law, that is their systems of finite dimensional distributions are equal.


A overview of Lévy processes and weak subordination is given in the appendix. Throughout, B=(B_1,…,B_n)∼ BM^n(,Σ) refers to an n-dimensional Brownian motion with linear drift [B(t)]= t and covariance matrix Cov(B(t))=tΣ, t≥0.
An n-dimensional subordinator T=(T_1,…,T_n)∼ S^n() is an n-dimensional Lévy process with nondecreasing components, and its Lévy measure is denoted by .


Gamma subordinator. For a,b>0, a univariate subordinator G∼Γ_S(a,b) is a gamma subordinator if its marginal G(t), t≥0,
is gamma distributed with shape parameter at and rate parameter b. 

If a=b, we refer to G as a standard gamma subordinator, in short, G∼Γ_S(b):=Γ_S(b,b).


Alpha-gamma subordinator. Assume n≥2. Let =(α_1,…,α_n)∈(0,∞)^n and G_0,…,G_n be independent gamma subordinators such that G_0∼Γ_S(a,1), G_k∼Γ_S(β_k,1/α_k), where a>0, aα_k<1, β_k:=(1-aα_k)/α_k, 1≤ k≤ n. A process T∼ AG_S^n(a,) is an alpha-gamma (AG_S^n) subordinator <cit.> with parameters a, if T G_0 +(G_1,…,G_n). An alpha-gamma subordinator T has correlated components with marginals T_k∼Γ_S(1/α_k), 1≤ k≤ n.

Variance-gamma process. Let b>0, ∈^n and Σ∈^n× n be a covariance matrix. A process V∼ VG^n(b,,Σ) is a variance-gamma (VG^n) process <cit.> with parameters b,,Σ if V∼ BM^n(,Σ)∘ (Γ_S(b)), where :=(1,…,1)∈^n.

The characteristic exponent of V is (see <cit.>, their Formula (2.9))

    Ψ_V()=-bln{ 1-/b+^2_Σ/2b} ,  ∈^n ,

where ln:\(-∞,0]→ is the principal branch of the logarithm.


Strong variance-alpha-gamma process. Assume n≥2. Let ∈^n and Σ∈[0,∞)^n× n be a diagonal matrix. A process X∼ VAG^n(a,,,Σ) is a (strong) variance-alpha-gamma (VAG^n) process <cit.> with parameters a,,,Σ if X∼ BM^n(,Σ)∘ AG^n_S(a,).


In <cit.>, the weak VAG process was formulate using weak subordination, allowing B to have dependent components while remaining a Lévy process.


Weak variance-alpha-gamma process. Assume n≥2. Let ∈^n and Σ∈^n× n be an arbitrary covariance matrix. The process X∼ WVAG^n(a,,,Σ) is a weak variance-alpha-gamma process <cit.> with parameters a,,,Σ if X∼ BM^n(,Σ)⊙ AG^n_S(a,), where ⊙ denotes the weak subordination operation (see the appendix).


Next, we gather various known results about the WVAG process that will be useful later on. The notation  is defined in (<ref>) and self-decomposability is defined in the appendix.

	Let n≥2 and X∼ WVAG^n(a,,,Σ).
	
		
  (i) X is an n-dimensional Lévy process with Lévy exponent
		
    Ψ_X()=   -a ln{ 1-+1/2^2_Σ}
       -∑_k=1^nβ_kln{ 1-α_kμ_kθ_k+1/2α_kθ_k^2Σ_kk},  ∈^n  .

		
  (ii) Let V_0∼ VG^n( a,a,aΣ), V_k∼ VG^1(β_k,(1-aα_k)μ_k,(1-aα_k)Σ_kk), 1≤ k≤ n be independent. Then X V_0+∑_k=1^nV_k_k.
		
  (iii)For any c>0, (X(ct))_t≥0∼ WVAG^n(ca,/c,c,cΣ).
		
  (iv) For 1≤ k≤ n, X has marginal distribution X_k∼ VG^1(1/α_k,μ_k,Σ_kk).
		
  (v) If Σ is diagonal, then X∼ VAG^n(a,,,Σ).
		
  (vi) For 1≤ k≠l≤ n, (X_k(1),X_l(1))=a(α_k∧α_l)Σ_kl+aα_kα_lμ_kμ_l.




		
  (vii) If Σ is invertible, then X is self-decomposable if and only if = 0.
	


See <cit.> for (i), (ii), (iv)-(vi), <cit.> for (iii), and <cit.> for (vii).

The WVAG process exhibits a wider range of dependence than the VAG process. For example, it has an additional covariance term a(α_k∧α_l)Σ_kl from Proposition <ref> (vi).




§ FOURIER INVERTIBILITY


Let n≥2 and X∼ WVAG ^n(a,,,Σ) and Y=I + X, ∈^n, where I:[0,∞)→[0,∞) is the identity function. The density function of Y(t), t>0, which is needed for MLE, exists because Y(t) has an absolutely continuous distribution for. However, it is not explicitly known, so it is computed using Fourier inversion as

    f_Y(t)() = (2π)^-n∫_ℝ^nexp(--𝐦)Φ_X(t)() ,  ∈^n ,

where Φ_X(t)()=exp(tΨ_X()), ∈^n, and Ψ_X() from (<ref>), provided Φ_X(t)∈ L^1. If Φ_X(t)∈ L^1, we say that X(t) is Fourier invertible and we give a condition for this in terms of an inequality relating the parameters.


	Let B∼ BM^n(,Σ), B^* ∼ BM^n( 0,Σ), T ∼ S^n(,), X  B ⊙ T, X^*  B^*⊙ T, Y I+X, ∈^n. For all t≥0 and p>0, if Φ_X^*(t)∈ L^p, then Φ_Y(t)∈ L^p.

 
	
	For all t≥0, Φ_Y(t)()=e^ tΦ_X(t)(), ∈^n, so that |Φ_Y(t)()|=exp(tΨ_X()). Using (<ref>), we have
	
    Ψ_X()   =∫_[0,∞)_*^n (Φ_B()()-1) ()
       ≤∫_[0,∞)_*^n (|Φ_B()()|-1) ()
       =∫_[0,∞)_*^n( exp(-1/2^2_Σ)-1)  ()
       =Ψ_X^*().

	Therefore, |Φ_Y(t)()|≤|Φ_X^*(t)()|, from which the result follows.




	Let V∼ VG^n(b,,Σ), and assume that Σ is invertible. Let p>0. If pb>n/2, then Φ_V∈ L^p.


 Since variance-gamma processes are weakly subordinated processes (see (<ref>)), we can apply Lemma <ref>, which means that we can assume =. For V∼ VG^n(b,0,Σ), by (<ref>), V has characteristic function
	
    Φ_V()=(1+_Σ^2/2b)^-b,  ∈^n.

	
	
	
	
	
	
	Using the Cholesky decomposition, Σ=U'U, where U is a lower triangular matrix with positive elements on the diagonal. Let p>0. Making the transformation =(2b)^1/2 (U')^-1, noting that (U')^-1 exists, and hence the transformation is injective, we have
	
    ∫_^n |Φ_V()|^p =|(2b)^1/2U^-1|∫_^n( 1+^2) ^-pb .

	Using the polar decomposition (see Corollary B.7.7 in <cit.>) on the RHS of (<ref>), we have Φ_V∈ L^p if and only if
	
    ∫_0^∞ (1+r^2)^-pbr^n-1  r<∞,

	
	which is equivalent to pb>n/2.





	Let X∼ WVAG ^n(a,,,Σ) and Y  I + X, ∈^n. Assume that Σ is invertible. For t>0, if
	
    (a/n+min_1≤ k ≤ nβ_k)t>1/2,

	then Φ_X(t),Φ_Y(t)∈ L^1.


 By Proposition <ref> (iii), it suffices to prove the result for t=1, and by Lemma <ref>, we can assume = and =, so that Y∼ WVAG ^n(a,, 0,Σ).
	
	Let V_0∼ VG^n(a,0, a Σ), V_k∼ VG^1(β_k,0,(1-aα_k)Σ_kk), 1 ≤ k ≤ n, be independent, and let V^*:=(V_1,…,V_n). By Proposition <ref> (ii), Y has characteristic function Φ_Y()=Φ_V_0()Φ_V^*(), where Φ_V^*():=∏_k=1^nΦ_V_k(θ_k). For p^-1+q^-1=1, p,q>1, 
	Hölder's inequality gives
	
    ∫_^n|Φ_Y()|    ≤( ∫_^n|Φ_V_0()|^p )^1/p( ∫_^n|Φ_V^*()|^q )^1/q
       = ( ∫_^n|Φ_V_0()|^p )^1/p∏_k=1^n( ∫_|Φ_V_k(θ)|^q θ)^1/q.

	By Lemma <ref>, this integral is finite when pa>n/2, qβ_k>1/2 and p,q>1 for all 1≤ k ≤ n. Thus,
	
    1=1/p+1/q<2(a/n∧1/2) +2( min_1≤ k ≤ nβ_k∧1/2) ,

	which is equivalent to (<ref>).














































 Note that this condition for a VG^1(b,μ,Σ) distribution to be Fourier invertible is identical to the condition for its density function having no singularity in <cit.>, which is b>1/2.


 We see that for sufficiently small t>0, (<ref>) will not be satisfied. This means that using (<ref>) to compute the density function may not be valid when attempting parameter estimation for a WVAG process based on observations from such a sufficiently small sampling interval.




§ CALIBRATION


We now specialise to the case of n=2. Let X∼ WVAG ^2(a,,,Σ), Y:=(Y_1,Y_2)=I +X, 𝐦∈^2. Let (S_1, S_2) be a bivariate price process

    S_k(t)=S_k(0)exp(Y_k(t)),   t≥0,   k=1,2.

For N equally spaced discrete observations with sampling interval c>0, the log returns are

    𝐲_j:=(y_1j,y_2j):=(lnS_1(jc)/S_1((j-1)c),
    	lnS_2(jc)/S_2((j-1)c))D= Y(c),   j=1,…, N,

and are iid. We call this the WVAG model. If Σ_12=0, we called it the VAG model as X reduces to a VAG process by Proposition <ref> (v).



 §.§ Simulation method


The result in Proposition <ref> (ii) can be used to simulate X∼ WVAG ^2(a,,,Σ) in terms of VG^n and VG^1 processes.

For the sampling intervals c=1,0.1 and sample size N=1000, we make 100 simulations of Y, and estimate the parameters from the observations (𝐲_j)_j=1^N with true parameters a=1, =(0.8,0.6), =(0.1,-0.3), Σ=[1,0.6;0.6,1.2], =(-0.1,0.3).




 §.§ Calibration methods


We estimate the parameters (a,,,Σ,) from the observations (𝐲_j)_j=1^N using method of moments (MOM), which is quick and easy to implement, maximum likelihood estimation (MLE) from Micha­elsen & Szimayer <cit.>, which may be expected as being asymptotically optimal under the model, and a modification of digital moment estimation (DME) from Madan <cit.>, which is more robust to model misspecification.

Method of moments. The initial values of μ_k,α_k,Σ_kk, m_k, k=1,2, are obtained by least squares on the first four central moments (Y_k(c)), ((Y_k(c)-(Y_k(c)))^p), p=2,3,4, with the corresponding sample moments. The initial values of the joint parameters a, Σ_12 are obtained by least squares on ((Y_1(c)-(Y_1(c)))^p(Y_2(c)-(Y_2(c)))^p), p=1,2, with the corresponding sample moments, with p=1 excluded when fitting the VAG model. Using these initial values, least squares is solved over all parameters. Note that this last step has no effect when these moments can be matched exactly.

Moment formulas can be found in <cit.>.


Maximum likelihood estimation.
The density function of Y(c) is not explicitly known so it is numerically computed using Fourier inversion by (<ref>). The numerical optimisation needed to implement MLE requires initial values. The first initial values can be obtained by MOM. Using the first initial values, MLE is applied to each marginal observations to obtain the second initial values of μ_k,α_k,Σ_kk, m_k, k=1,2, and to the bivariate observations to obtain the second initial values of a, Σ_12. Finally, using the second initial values, MLE is applied on all parameters. For the VAG model, we apply the above method with the constraint Σ_12=0.

Digital moment estimation.
Let k=1,2, let 𝐪 be the vector of 10 equally spaced points from 0.05 to 0.95, and let ℙ_k be the empirical quantiles of the observations (y_kj)_j=1^N at the probabilities 𝐪. Let p_y(μ_k,α_k,Σ_kk,m_k) := (Y_k(c)≤ y), where y∈ℙ_k, Y_k∼ m_kI+VG^1(1/α_k,μ_k,Σ_kk) (see Proposition <ref> (iv)), and q_y is the corresponding empirical probability. Marginal parameters μ_k,α_k,Σ_kk, m_k are estimated by minimizing the error ∑_y∈ℙ_k (p_y(μ_k,α_k,Σ_kk,m_k)-q_y)^2.
		
With the estimated marginal parameters, let ρ:=Σ_12/(Σ_11Σ_22)^1/2, p_𝐲(a,ρ):= (Y_1(c)≤ y_1, Y_2(c)≤ y_2), where 𝐲:=(y_1,y_2)∈ℙ_1×ℙ_2, Y∼ I+WVAG ^2(a,,,Σ), and q_𝐲 is the corresponding empirical probability. Since p_𝐲(a,ρ) is computationally expensive to calculate directly, it is estimated by the empirical probability over 10000 simulations. The joint parameters a, ρ are estimated by minimizing the LOESS smooth <cit.> of the error ∑_𝐲∈ℙ_1×ℙ_2 (p_𝐲(a,ρ)-q_𝐲)^2. The predictor variables for the LOESS smooth are 100 equally spaced points on the feasible set of (a,ρ)∈(0,(1/α_1)∧(1/α_2))×(-1,1). For the VAG model, we apply the above method with the constraint ρ=0.




 §.§ Goodness of fit statistics


To assess the overall goodness of fit of each parameter estimation method, as opposed to assessing individual parameters, we consider 3 goodness of fit statistics, the negative log-likelihood (-log L), a chi-squared (χ^2) statistic, and a Kolmogorov-Smirnov (KS) statistic.

To compute χ^2, we apply the Rosenblatt transform <cit.> of the fitted distribution to the observations, which has a uniform distribution on [0,1]^2 if the fitted distribution coincides with the true distribution, and then we compute the χ^2 statistic for a test of uniformity over an equally spaced partition of [0,1]^2 into 100 cells. Since computing -log L and χ^2 requires Fourier inversion, it may not be possible to compute these statistics accurately when the Fourier invertibility condition does not hold, so they are not displayed in Table <ref>.

Therefore, we also consider the 2-dimensional, two-sample Kolmogorov-Smirnov statistic introduce by Peacock in <cit.>, and computed using the method of <cit.>. This is the statistic for testing equality of the fitted distribution and the true distribution based on a sample from the respective distributions, and therefore does not require the density function f_Y(c)() or Fourier inversion. When applied to real data in Subsection <ref>, we take the average of the KS statistics computed from the observations and 100 samples from the fitted distribution. When applied to simulated data in Subsection <ref>, the KS statistic is computed from the observations and a sample from the fitted distribution. All 3 goodness of fit statistics were averaged over the 100 simulations.



 §.§ Quantile choice for DME


Different choices of quantiles for DME are possible. Let 𝐪_1 be the vector of 10 equally spaced points from 0.05 to 0.95, 𝐪_2 be the vector of 10 equally spaced points from 0.01 to 0.99, 𝐪_3 be the vector of 10 equally spaced points from 0.1 to 0.9, 𝐪_4 be the vector of 20 equally spaced points from 0.05 to 0.95. For sampling interval c=1, Table <ref> shows the goodness of fit for 4 choices of quantiles. We find 𝐪=𝐪_1 as yielding the lowest RMSE for most variables and the lowest goodness of fit statistics. However, given that the results are so similar, these quantile choices make only a small difference to the overall goodness of fit.




 §.§ Simulated data results


For the sampling interval c=1, the Fourier invertibility condition is satisfied as the LHS of (<ref>) is 0.75>1/2. The calibration results for the WVAG model with c=1 is shown in Table <ref>. Here, we find that MLE gives the best fit with the lowest χ^2 statistic. The KS statistic for MLE and DME are approximately equal.



For the sampling interval c=0.1, the Fourier invertibility condition is violated as the LHS of (<ref>) is 0.08<1/2. The corresponding results are shown in Table <ref>. Here, we find that DME gives the best fit with the lowest KS statistic, however MLE still produces a good fit and does not break down. This suggests that the condition may not be necessary for the MLE to produce accurate parameter estimates. In both cases, c=1,0.1, the RMSE and goodness of fit statistics are highest for MOM.





 §.§ Real data results

Next, we fit the WVAG and VAG models to the S&P500 and FTSE100 indices as the bivariate price process (<ref>) for a 5 year period from 14 February 2011 to 12 February 2016 with daily closing price observations taking c=1. The estimated parameters, goodness of fit statistics and standard errors computed using 100 bootstrap samples are listed in Table <ref>. Contour plots of the fitted distributions and scatter plots of the bivariate log returns are shown in Figure <ref>.

[!htbp]
	
	

		
			    2cMOM     2cMLE     2cDME
 Parameter     WVAG     VAG     WVAG     VAG     WVAG     VAG 
 a                     -0.695 (0.199)     -0.696 (0.217)     -0.962 (0.139)     -1.087 (0.117)     -0.899 (0.138)     -1.114 (0.327) 

			α_1              -1.436 (0.398)     -1.436 (0.349)     -0.919 (0.105)     -0.908 (0.104)     2c-0.898 (0.143) 

			α_2              -0.933 (0.277)     -0.707 (0.120)     -0.837 (0.105)     -0.919 (0.085)     2c-0.878 (0.132) 

			1000μ_1             -1.153 (0.704)               -1.156 (0.595)               -0.403 (0.550)               -0.630 (0.694)               2c-0.562 (0.563) 

			1000μ_2             -1.246 (0.704)               -1.255 (0.789)               -0.891 (0.586)               -0.867 (0.578)               2c-1.166 (0.554) 

			10000Σ_11      -0.982 (0.070)     -0.980 (0.069)     -0.976 (0.061)     -0.935 (0.064)     2c-0.928 (0.071) 

			10000Σ_22      -1.006 (0.053)     -1.007 (0.062)     -1.028 (0.063)     -1.014 (0.060)     2c-1.051 (0.086) 

			10000Σ_12      -0.994 (0.072)     --                   -0.813 (0.089)     --                   -0.844 (0.095)     -- 

			1000m_1               -1.422 (0.604)     -1.426 (0.519)     -0.705 (0.430)     -0.879 (0.575)     2c-0.982 (0.444) 

			1000m_2               -1.198 (0.617)     -1.207 (0.722)     -0.847 (0.440)     -0.882 (0.470)     2c-1.066 (0.410) 
 -ln L                -8465.969                      -8192.168                      -8496.315                      -8239.301                      -8492.798                      -8237.767 

			χ^2                -144.034             -702.241             -118.574             -586.789             -99.359              -584.547 

			KS                      -0.073               -0.151               -0.050               -0.139               -0.048               -0.138 

			
		Parameter estimates and standard errors (in parentheses) when fitted to the S&P500-FTSE100 data set. Estimation for the marginal parameters of WVAG and VAG models are identical.
		
	

	




Note that the Fourier invertibility condition is satisfied for all fitted models. Based on the χ^2, KS statistic and contour plots, the WVAG model produces a better fit than the VAG model. In addition, for the WVAG model, DME gives a fit with a lower χ^2 and KS statistic than MLE and MOM.

Assuming that the log returns satisfies the WVAG model, a likelihood ratio test can be used to test the hypothesis H_0:Σ_12=0 versus H_1:Σ_12≠0.
The test statistic D=514.03 is asymptotically χ^2 distributed with 1 degree of freedom. The p-value is <10^-4, so the VAG model is rejected.
Indeed, the VAG model is not suited for modelling strong correlation since Cov((B∘ T)_1(1), (B∘ T)_2(1))=aα_1α_2μ_1μ_2 by Proposition <ref> (vi), which is approximately 0 when μ_1μ_2 is.

It has been suggested that log-returns should be self-decomposable <cit.>. Note that =(-0.0004,-0.0008) is very close to 0, which suggests that the log-returns process Y is indeed self-decomposable (see Proposition <ref> (vii)). A likelihood ratio test can be used to test this hypothesis, H_0:= 0 versus H_1:≠ 0. The test statistic D=4.11 is asymptotically χ^2 distributed with 2 degrees of freedom. The p-value is 0.128, so at a 5% significance level we cannot reject that Y is self-decomposable.



§ CONCLUSION


The WVAG process constructed by using weak subordination generalises the VAG process, and we obtain a condition for Fourier invertibility in Theorem <ref>. We have shown that MOM, MLE and DME can be used to estimate the parameters of a WVAG process, and find that in our simulations MLE produces a better fit when the Fourier invertibility condition holds, while DME produces a better fit when it is violated. However, MLE may still produce good parameter estimates even when the Fourier invertibility condition is violated. In all cases, MOM produces the worst fit. We find that the WVAG process exhibits a wider range of dependence and produces a significantly better fit than the VAG process when used to model the S&P500-FTSE100 data set, and that DME produces the best fit in this situation.



§ APPENDIX

Lévy process. The reader is referred to the monographs <cit.> for necessary material on Lé­vy pro­cesses, to <cit.> for financial applications, while
our notation follows <cit.>. For A⊆^n, let A_*:=A\{ 0} and let _A(ω) denote the indicator function. Let :={∈^n:≤ 1} be the Euclidean unit ball centred at the origin. The law of an n-dimensional Lévy process X=(X_1,…,X_n)=(X(t))_t≥ 0 is determined by its characteristic function Φ_X:=Φ_X(1), with

    Φ_X(t)() := exp(X(t)) = exp(tΨ_X()) ,      t≥0,

and Lévy exponent Ψ_X:=Ψ, where

    Ψ() := -1/2 ^2_Σ
    +∫_^n_*(e^-1-_()) () ,

∈^n, ∈^n, Σ∈^n× n is a covariance matrix,
and  is a nonnegative Borel measure on ^n_* such that
∫_^n_*(^2∧ 1) ()<∞. We write X∼ L^n(,Σ,) provided X is an n-dimensional Lévy process with canonical triplet (,Σ,).

A subordinator T∼S^n()=L^n(,0,) is drift-less if its drift -∫__* ()= 0. All subordinators considered in this paper are drift-less.

An n-dimensional random variable X is self-decomposable if for any 0<b<1, there exists a random variable Z_b, independent of X, such that X bX+Z_b. A Lévy process X is self-decomposable if X(1) is.


Strongly subordinated Brownian motion. Let B=(B_1,…,B_n)∼ BM^n(,Σ) be a Brownian motion and T=(T_1,…,T_n)∼S^n() be a drift-less subordinator. A process B∘ T is the traditional or strong subordination of X and T if (B∘ T)(t):=(B_1(T_1(t)),…,B_n(T_n(t))), t≥ 0.

Weakly subordinated Brownian motion. Let =(t_1,… t_n)∈[0,∞)^n, =(μ_1,…μ_n)∈^n and Σ=(Σ_kl)∈^n× n be a covariance matrix. Introduce the outer products ∈^n and Σ∈^n× n by

    :=(t_1μ_1,…,t_nμ_n)  and  (Σ)_kl:=Σ_kl(t_k∧ t_l),   1≤ k,l≤n.

Let B∼BM^n(,Σ) be an n-dimensional Brownian motion and T∼S^n() be an n-dimensional drift-less subordinator. A Lévy process B⊙ T is called the
weak subordination of B and T (see <cit.>, their Proposition 3.1) if it has Lévy exponent

    Ψ_B⊙ T()=∫_[0,∞)^n_*(exp( -1/2^2_Σ)-1) () ,

∈^n. Note that a more general definition of weak subordination and a proof of existence is given in <cit.>.

Assume that independent B and T. If T has indistinguishable components or B has independent components, then 

    B∘ T B⊙ T

Otherwise B∘ T may not a Lévy process, but B⊙ T always is (see <cit.>, their Proposition 3.3 and 3.9).



§ ACKNOWLEDGEMENT

B. Buchmann's research was supported by ARC grant DP160104737. K. Lu's research was supported by an Australian Government Research Training Program Scholarship.

9
Ap09 Applebaum, D. (2009). Lévy Processes & Stochastic Calculus.
Cambridge Studies in Advanced Mathematics, 116, 2nd ed, Cambridge University Press, Cambridge. MR2512800


BB14 Ballotta, L. & Bonfiglioli. E. (2016). Multivariate asset models using Lévy processes and applications. European J. Finance 22, 1320–1350.

BPS01 Barndorff-Nielsen, O.E., Pedersen, J. & Sato, K. (2001). Multivariate subordination, self-decomposability and stability. Adv. in Appl. Probab. 33, 160–187. MR1825321

BS10
Barndorff-Nielsen, O.E. & Shiryaev, A. (2010). Change of Time and Change of Measure. World Scientific Publishing Co. Pte. Ltd., Hackensack, NJ. MR2779876

bBertoin, J. (1996). Lévy Processes. Cambridge University Press, Cambridge. MR1406564

Bi06 Bingham, N.H. (2006). Lévy processes and self-decomposability in finance. Probability and Mathematical Statistics 26, 131–142.

BiKi02 Bingham, N.H. & Kiesel, R. (2002). Semi-parametric modelling in finance: Theoretical foundations. Quantitative Finance 2, 241–250.

BogaBogachev, V. I. (2007). Measure theory (Vol. 1). Springer Science & Business Media.

BKMS16 Buchmann, B., Kaehler, B., Maller, R., Szimayer, A. (2017). Multivariate subordination using generalised gamma convolutions with applications to variance gamma processes & option pricing. Stochastic Process. Appl. 127, 2208-–2242.

BLMa Buchmann, B., Lu, K., Madan, D. (2017). Weak subordination of multivariate Lévy processes and variance generalised gamma convolutions. To appear in Bernoulli. Available at https://arxiv.org/abs/1609.04481
BLMb Buchmann, B., Lu, K., Madan, D. (2017). Self-decomposability of variance generalised gamma convolutions Preprint. Australian National University and University of Maryland. Available at https://arxiv.org/abs/1712.03640

CS09 Cariboni, J. & Schoutens, W. (2009). Lévy Processes in Credit Risk. Wiley, New York.

CGMY07 Carr, P., Geman, H., Madan, D.B. & Yor, M. (2007). Self-decomposability and option pricing. Mathematical Finance 17, 31–57.

CGS Cleveland, W.S., Grosse, E. & Shyu, W.M. (1991). Local Regression Models. Chapter 8, Statistical Models in S, editors Chambers, J.M. & Hastie, T.J., Chapman & Hall/CRC, Boca Raton.

CT Cont, R. & Tankov, P. (2004). Financial Modelling with Jump Processes, Chapman & Hall. London, New York, Washington D.C. MR2042661.

FiS Finlay, E., & Seneta, E. (2008). Stationary-increment Variance-Gamma and t models: Simulation and parameter estimation. Int. Stat. Rev. 76, 167–186. MR2492088.

FS07 Fung, T., & Seneta, E. (2010). Modelling and estimation for bivariate financial returns. International statistical review 78, 117–133.

Gu13 Guillaume, F. (2013). The αVG model for multivariate asset pricing: calibration and extension. Rev. Deriv. Res. 16, 25–52.

Ka09 Kawai, R. (2009). A multivariate Lévy process model with linear correlation. Quantitative Finance 9, 597–606.

KuTaKüchler, U. & Tappe, S. (2008). On the shapes of bilateral Gamma densities. Statistics & Probability Letters 78, 2478–2484.

LMS16Luciano, E., Marena, M., & Semeraro, P. (2016). Dependence calibration and portfolio fit with factor-based subordinators. Quantitative Finance 16, 1037–1052.


LS10Luciano, E. & Semeraro, P. (2010). Multivariate time changes for Lévy asset models: Characterization and calibration. J. Comput. Appl. Math. 233, 1937–1953. MR2564029


Ma11 Madan, D.B. (2011). Joint risk-neutral laws and hedging. IIE Transactions 43, 840–850.

Ma15 Madan, D.B. (2015). Estimating parametric models of probability distributions. Methodol. Comput. Appl. Probab. 17, 823–831. MR3377863

Ma18 Madan, D.B. (2018). Instantaneous portfolio theory. Quantitative Finance. Available at https://doi.org/10.1080/14697688.2017.1420210

MCC98 Madan, D.B., Carr, P.P., & Chang, E.C. (1998). The variance gamma process and option pricing. Review of Finance 2, 79–105.

MaSe90 Madan, D.B. & Seneta, E. (1990).
The variance gamma (v.g.). model for share market returns. Journal of Business 63, 511–524.


Mi18 Michaelsen, M. (2018). Information flow dependence in financial markets. Preprint. Universität Hamburg. Available at https://ssrn.com/abstract=3051180

MiSz17 Michaelsen, M. & Szimayer, A. (2018). Marginal consistent dependence modelling using weak subordination for Brownian motions. Quantitative Finance. Available at https://doi.org/10.1080/14697688.2018.1439182

peac Peacock, J.A. (1983). Two-dimensional goodness-of-fit testing in astronomy. Mon. Not. R. astr. Soc. 202, 615–627.

R52Rosenblatt, M. (1952). Remarks on a multivariate transformation. Ann. Math. Statist. 23, 470–472. MR0049525

sas13 Sasvári. Z. (2013). Multivariate Characteristic and Correlation Functions. Walter de Gruyter, Berlin, 2013. MR3059796

s Sato, K. (1999). Lévy Processes and Infinitely Divisible Distributions.
Cambridge University Press, Cambridge. MR3185174

Se08 Semeraro, P. (2008). A multivariate variance gamma model for financial applications. Journal of Theoretical and Applied Finance 11, 1–18. MR2398464

Wa09 J. Wang. (2009). The Multivariate Variance Gamma Process and Its Applications in Multi-asset Option Pricing. PhD Thesis, University of Maryland.

xiao Xiao, Y. (2017). A fast algorithm for two-dimensional Kolmogorov-Smirnov two sample tests. Computational Statistics & Data Analysis 105, 53–58.







