The Canadian Journal of Statistics Vol. 35, No. 2, 2007, Pages ???-??? La revue canadienne de statistique
Nonparametric estimation of copula functions
for dependence modelling
Song Xi CHEN and Tzee-Ming HUANG
Key words and phrases: boundary bias; copula models; dependence; kernel estimator. MSC 2000 : Primary 62G05; secondary 62E20.
Abstract: Copulas are full measures of dependence among components of random vectors. Unlike the
marginal and the joint distributions, which are directly observable, a copula is a hidden dependence structure that couples a joint distribution with its marginals. This makes the task of proposing a parametric copula model non-trivial and is where a nonparametric estimator can play a significant role. In this paper, we propose a kernel estimator which is mean square consistent everywhere in the support of the copula function. The bias and variance of the copula estimator are derived which reveal the effects of kernel smoothing on the copula estimation. A smoothing bandwidth selection rule based on the derived bias and variance is proposed. The theoretical findings are confirmed by a simulation study. The kernel estimator is then used to formulate a goodness-of-fit test for parametric copula models.
Title in French: we can supply this
R´esum´e : Ins´erer votre r´esum´e ici. Copulas are full measures of dependence among components of random vectors. Unlike the marginal and the joint distributions which are directly observable, a copula is a hidden dependence structure that couples a joint distribution with its marginals. This makes the task of proposing a parametric copula model non-trivial and is where a nonparametric estimator can play a significant role. In this paper, we propose a kernel estimator which is mean square consistent everywhere in the support of the copula function. The bias and variance of the copula estimator are derived which reveal the effects of kernel smoothing on the copula estimation. A smoothing bandwidth selection rule based on the derived bias and variance is proposed. The theoretical findings are confirmed by a simulation study. The kernel estimator is then used to formulate a goodness-of-fit test for parametric copula models.
1. INTRODUCTION
Quantifying the dependence among two or more random variables has been an enduring task for statisticians. A rich set of dependence measures has been proposed, including the well-known Pearson’s correlation coefficients, Kendall’s tau and Spearman’s rho for bivariate random variables. While these measures are simple and can be easily computed, they are designed to capture only certain aspects of dependence. Indeed, it is rather unreasonable to expect a single scalar measure to have the capability to quantify all the dependence existing among the random variables.
Copula is a device that fully quantifies the dependence among random variables. Let X = (X1, . . . , Xd)τ be a random vector, and F be the distribution function of X with marginal
existence of a multivariate distribution function C on [0, 1]d such that
F (x1, . . . , xd) = C{F1(x1), . . . Fd(xd)}.
The function C is called the copula associated with X and couples the joint distribution F with its marginals. See Nelsen (1998) for a comprehensive overview of copulas and their mathematical properties.
The implication of the Sklar’s Theorem is that, after standardizing the effects of marginals, the dependence among components of X is fully described by the copula. Indeed, most conven-tional measures of dependence can be explicitly expressed in terms of the copula. For example,
the Kendall’s tau between X1 and X2 is 4 [0,1]2C(u, v)dC(u, v)− 1, and Spearman’s rho is
12 [0,1]2C(u, v)dudv− 3. Copulas can also be used to describe tail dependence (Joe 1997), an important notion in risk management. The interest there is in the dependence between two extreme (risky) events.
An early statistical application of copulas can be found in Clayton (1978), where the dependence
between two survival times in a multiple events study is modelled by C(u, v) ={u−1/θ+ v−1/θ−
1}1/(θ−1), the so-called Clayton copula. In later research into copulas, a driving force has been in financial risk management for modelling dependence among different assets in a portfolio; see Embrechts, Lindskog & McNeil (2003) for a comprehensive review.
Estimation of copulas can be achieved fully parametrically by assuming parametric models for both the copula and the marginals and then performing maximum likelihood estimation (Oakes 1982, in the context of Clayton copula). Semiparametric estimation that specifies a parametric copula while leaving the marginals nonparametric is proposed in Genest, Ghoudi & Rivest (1995) and Chen & Fan (2006). Estimation for Archimedean copulas in Genest & Rivest (1993) can be considered semiparametric as well. Recently, parametric copula models have been used for more general situations of inference. For instance, Fine & Jiang (2000) considered estimating the parameter of a Clayton copula with covariates in the marginal distributions in the context of life time data.
A nonparametric estimation of copula treats both the copula and the marginals parameter-free and thus offers the greatest generality. Unlike the marginal and the joint distributions which are directly observable, a copula is a hidden dependence structure. This makes the task of proposing a suitable parametric copula model non-trivial and is where a nonparametric estimator can play a significant role. Indeed, a nonparametric copula estimator can provide initial information needed in revealing and subsequent formulation of an underlying parametric copula model.
Nonparametric estimation of copulas dates back to Deheuvels (1979), who proposed an es-timator based on a multivariate empirical distribution on the marginal empirical distributions. Smoother estimators based on the kernel method have been proposed in the literature. Gijbels & Mielnicnuk (1990) proposed a kernel estimator for a bivariate copula density, which is employed by Fermanian (2005) in a goodness-of-fit test. Another approach of kernel estimation is to estimate a copula function directly as explored in Fermanian & Scaillet (2003). This approach is the focus of our investigation in this paper. The advantage of targeting directly on a copula instead of its
density is that the estimation has a faster rate of convergence at n−1/2, which is not affected by
the dimension of the copula. This faster rate can lead to a less variable goodness-of-fit test statistic that is formulated based on the kernel copula estimator. Of course, whether to estimate a copula or its density depends on the particular statistical inference problem we encounter. One approach complements rather than rule out the prospects of the other.
A critical issue confronted with both kernel copula and copula density estimation is that a
copula and its density are defined on a compact cube [0, 1]d. This means that the boundary bias
associated with kernel curve estimation will be present. Indeed, the multivariate boundaries require more cares than the standard univariate boundary region. An analysis for kernel estimation with multivariate boundary regions is given in M¨uller & Stadtm¨uller (1999). For the bivariate case, it is
necessary to ensure consistent estimation of a copula function over entire [0, 1]2 especially near the corners (0, 0) and (1, 1). This is important in studying dependence between two extreme events, which is a target application of copula in risk management.
We propose in this paper a new bivariate kernel copula estimator based on local linear kernels and a simple mathematical correction that removes the boundary bias. We then derive the bias and variance of this estimator, which reveal that the kernel smoothing produces a second order reduction in both the variance and mean square error as compared with the unsmoothed empirical estimator of Deheuvels (1979). Our analysis also identifies which part of the kernel smoothing is the source of this variance reduction and provides a practical guideline for smoothing bandwidth selection.
The paper is organized as follows. The kernel estimator is proposed in Section 2. Its bias and variance are reported in Section 3. Section 4 considers bandwidth selection. Simulation results are reported in Section 5. Section 6 reports an empirical study. All the technical details are given in the appendix.
2. A KERNEL ESTIMATOR
The basic thrust for our kernel copula estimator is the fact that, when the two marginal distribu-tions are continuous, the copula C is the unique joint distribution of F1(X1) and F2(X2) as implied
by Sklar’s theorem. As copulas are not directly observable, a nonparametric copula estimator has
to be formed in two stages: estimate the two marginals (F1(X1), F2(X2)) first and then estimate
the copula based on the estimated marginals.
Let K be a symmetric probability density supported on [−1, 1] and G(x) =−∞x K(t)dt be the
distribution of K. In the first stage the marginal distribution Flis estimated by
ˆ Fl(x) = n−1 n i=1 G{(x − Xil)/bl}
with a bandwidth blfor l = 1 and 2; see Bowman, Hall & Prvan (1998) for details on this kernel
distribution estimator.
To prevent the boundary bias, we use in the second stage Ku,h(x) = K(x){a2(u, h)− a1(u, h)x}
a0(u, h)a2(u, h)− a21(u, h)
,
a local linear version of K, to smooth at a u ∈ [0, 1] with a bandwidth h > 0. Here al(u, h) =
u/h
(u−1)/htK(t)dt for l = 0, 1 and 2, which was proposed by Lejeune and Sarda (1992) and Jones
(1993) and was designed to remove the boundary bias in univariate density estimation. It is easy to check that Ku,h= K for u∈ [h, 1 − h].
Let Gu,h(t) =
t
−∞Ku,h(x)dx and Tu,h= Gu,h{(u − 1)/h}. A seemingly natural estimator of
C(u, v) would be n−1 n i=1 Gu,h u− ˆF1(Xi1) h Gv,h v− ˆF2(Xi2) h .
The estimator considered by Fermanian & Scaillet (2003) is of the above form based on the original kernel K rather than the local linear kernels Ku,h.
It can be seen from the bias expression given in (2) in the next section that the use of the local linear kernel removes the boundary bias near u = 0 and v = 0. However, it still incurs a
bias uTv,h+ vTu,h+ Tu,hTv,h near u = 1 or v = 1 due to the fact that each marginal distribution
kernel and h, the bias can be easily removed by subtraction. This leads to the proposed kernel copula estimator ˆ C(u, v) = n−1 n i=1 Gu,h u− ˆF1(Xi1) h Gv,h v− ˆF2(Xi2) h
−(uTv,h+ vTu,h+ Tu,hTv,h). (1)
It is noted that a single bandwidth h is used to smooth ˆFl(Xil) for l = 1 and 2 in the second stage, as the quantile transformation has already achieved a uniform data standardization.
The above boundary bias correction is different from the one used in Gijbels & Mielnicnuk (1990) for copula density estimation. The correction carried out there was rather involved via generating eight extra copies of the sample by reflecting the original data with respect to the four
edges and four corners of the unit square [0, 1]2 respectively, followed by applying the standard
kernel density estimation technique to certain linear transformations of the estimated marginals on the extended sample.
3. MAIN RESULTS
The study of the copula estimator faces two challenges. One is that the estimator is based on the estimated marginal ˆFl(Xi) instead of Fl(Xi). The other is that using the local linear kernels further increases the labor of derivations. However, we are able to obtain tractable expressions for the bias and variance, which add to the existing convergence results for copula density estimation (Gijbels & Mielnicnuk 1990; Fermanian, Radulovic & Wegkamp 2004) and for copula estimation (Fermanian & Scaillet 2003), and provide a finer scale description of the sampling properties of the copula estimation.
The following conditions are assumed in our analysis:
A1: K is a symmetric and continuous probability density supported on [−1, 1], and the
bandwidths satisfy h = O(n−1/3) and bl= O(h) for l = 1 and 2.
A2: For l = 1 or 2, Xl has a probability density function flsuch that f (1)
l and f
(2)
l , the first
and second derivatives of fl, are bounded and vanish at±∞. Furthermore, each of the two
limits limx→−∞fl(2)(x)/fl(1)(x) and limx→∞fl(2)(x)/fl(1)(x) is either a finite number or±∞.
A3: The copula C has a probability density function f and there exists a C∞ function g
such that f = g on [0, 1]2.
We note that Condition A3 does not imply that the copula density is infinitely differentiable. Rather there is an infinitely differentiable function g which is identical to f on [0, 1]2. For instance, for the independent copula, g is the constant 1. However, this condition does put a strong restriction on the smoothness of copula near the boundaries, as pointed out by one referee. In particular, if the copula density is unbounded, then the condition does not hold. For instance, Condition A3 does not hold for the Gumbel copula used in our simulation study.
Let Cu(u, v) and Cv(u, v) be the first and Cuu(u, v) and Cvv(u, v) be the second partial
deriva-tives of C(u, v) with respective to u and v respectively. Let ν(u, h) =(u−1)/hu/h s2dGu,h(s), which equals σ2K =: s2K(s)ds for u∈ [h, 1 − h]. Also, for l = 1 and 2, define
μl(v, h, λ) = v h v−1 h v h v−1 h max{r1+ fl(Fl−1(v))λw1, r2+ fl(Fl−1(v))λw2} × ×dGv,h(r2)dGv,h(r1)dG(w1)dG(w2).
and μ∗l(v, h, λ) = vh v−1 h v h v−1 h max{t, r + fl(F −1 l (v))λw}dGv,h(r)dGv,h(t)dG(w).
We have the following proposition, which quantifies the bias and variance of ˆC(u, v).
Proposition 1. Under Conditions A1-A3, for any (u, v) ∈ [0, 1]2,
E{ ˆC(u, v)}
= C(u, v) + 12h2{Cuu(u, v)ν(u, h) + Cvv(u, v)ν(v, h)} + o(h2) (2)
− 1
2σ2K
{Cu(u, v) + Tv,h}f (1)
1 {F1−1(u)}b21+{Cv(u, v) + Tu,h}f (1)
2 {F2−1(v)}b22
; and V ar{ ˆC(u, v)}
= n−1V ar{I(U ≤ u, V ≤ v) − Cu(u, v)I(U ≤ u) − Cv(u, v)I(V ≤ v)}
− hn−1 {C
u(u, v) + Tv,h}2μ1(u, h, b1/h) +{Cv(u, v) + Tu,h}2μ2(v, h, b2/h)
+ 2hn−1
{Cv(u, v) + Tu,h}2μ∗2(v, h, b2/h) +{Cu(u, v) + Tv,h}2μ∗1(u, h, b1/h)
− hn−1{C u(u, v)(1 + 2Tv,h) + Tv,h2 } u h u−1 h sdG2u,h(s) − hn−1{C
v(u, v)(1 + 2Tu,h) + Tu,h2 }
v
h v−1
h
tdG2v,h(t) + o(hn−1). (3)
A sketch of proof of Proposition 1 is given in the appendix.
While the bias (2) conveys a simple story that both the first and the second stage smoothing contribute to the bias, the variance given in (3) requires a further analysis. Let us concentrate on (u, v)∈ [h, 1 − h]2, the interior region as we use a compact kernel supported in [−1, 1].
Let xu= f1{F1−1(u)}, yv= f2{F2−1(v)} and bK =
tdG2(t). Then (3) can be simplified to
V ar{ ˆC(u, v)} = n−1V ar{I(U ≤ u, V ≤ v) − C
u(u, v)I(U ≤ u) − Cv(u, v)I(V ≤ v)}
+ hn−1{Cu2(u, v)ρ∗(xub1/h) + Cv2(u, v)ρ∗(yvb2/h)}
− hn−1{C
u(u, v)− Cu2(u, v) + Cv(u, v)− Cv2(u, v)}bK+ o(hn−1), (4)
where ρ∗(λ) = 2μ∗(λ)− μ(λ) − bK, μ∗(λ) = max{t, r + λw}dG(r)dG(t)dG(w) and μ(λ) = max{r1+ λw1, r2+ λw2}dG(r2)dG(r1)dG(w1)dG(w2).
A key fact needed in understanding (4) is that
ρ∗(λ)≥ 0 for any λ ≥ 0 and it is minimized at λ = 0. (5) In order to achieve the largest variance reduction, we need to minimize the second term on
the right hand side of (4) which involves the ρ∗-function. Due to (5), our strategy is to choose
bl= o(h) for l = 1 and 2 so that both ρ∗(xub1/h) and ρ∗(yvb2/h) are o(1), and hence
V ar{ ˆC(u, v)} = n−1V ar{I(U ≤ u, V ≤ v) − C
u(u, v)I(U ≤ u) − Cv(u, v)I(V ≤ v)}
− hn−1b
K{Cu(u, v)− Cu2(u, v) + Cv(u, v)− Cv2(u, v)} + o(hn−1). (6)
This indicates a second order variance reduction by the second stage smoothing by noting that the
by Deheuvels (1979): ˜ C(u, v) = n−1 n i=1 I( ˆUi≤ u, ˆVi≤ v) (7) where ˆUi = n−1 n j=1I(Xj1 ≤ Xi1) and ˆVi = n−1 n
j=1I(Xj2 ≤ Xi2). Despite the variance
reduction happens in the interior region only, it leads to a net reduction in the overall MISE over [0, 1]2 as shown in the next section.
A drawback of ˜C(u, v) is its lack of continuity. This lack of continuity, as indicated in a
simulation study reported in Section 5, can produce as twice the mean integrated square error (MISE) as that of the proposed kernel estimator. This shows that the variance and MISE reductions by the kernel estimator are significant in finite samples.
Finally, we note that the results in (2) and (3) apply to a kernel copula estimator that uses the marginal empirical distributions instead of the kernel distribution estimators in (1) by simply setting b1= b2= 0. In this case, (2) and (3) are simplified to
E{ ˆC(u, v)} = C(u, v) +1
2h2{Cuu(u, v)ν(u, h) + Cvv(u, v)ν(v, h)} + o(h2)
and (6) respectively.
4. BANDWIDTH SELECTION
The findings of the previous section suggest that we should undersmooth in the first stage to
reduce the variance, namely bl should be o(h) for l = 1 and 2. This reduces the bias from the first
stage smoothing too. This strategy largely simplifies the expressions of (2) and (4) and leads to a
tractable expression for the mean square error (MSE) for (u, v)∈ [h, 1 − h]2
M SE{ ˆC(u, v)} = n−1V ar{I(U ≤ u, V ≤ v) − Cu(u, v)I(U ≤ u) − Cv(u, v)I(V ≤ v)} − hn−1b
K
Cu(u, v){1 − Cu(u, v)} + Cv(u, v){1 − Cv(u, v)}
+ 14h4σ4k{Cuu(u, v) + Cvv(u, v)}2+ o(h4+ hn−1).
As the area of the boundary regions are of O(h) and the leading variance term is valid through-out the entire [0, 1]2, the MISE of ˆC is
M ISE( ˆC) = n−1
1
0
1
0 V ar{I(U ≤ u, V ≤ v) − Cu(u, v)I(U ≤ u)
−Cv(u, v)I(V ≤ v)}dudv − hn−1α +14h4σk4β + o(hn−1+ h4)
where β =0101{Cuu(u, v) + Cvv(u, v)}2dudv and
α = bK 1 0 1 0
Cu(u, v){1 − Cu(u, v)} + Cv(u, v){1 − Cv(u, v)}
dudv. The optimal h that minimizes the above MISE is then
h∗= σK4/3(α/β)1/3n−1/3. (8)
Various plug-in bandwidth selection rules that have been used in kernel smoothing can be employed here to attain an estimate for the optimal bandwidth. A simple approach is to assume a parametric family for the copula function which then leads to parametric expressions for α and β. This is similar to the reference rule suggested by Silverman (1986) for kernel density estimation. The parameters of the parametric copula can be estimated by either the pseudo-maximum likelihood estimation or the method of moments. Then we plug-in the α and β estimates from the reference rule to (8) to obtain an estimate of the optimal bandwidth.
We propose the following T -copula as the reference copula: C(u, v) = t−1ν (u) −∞ t−1ν (v) −∞ 1 2π(1− ρ2)1/2{1 + s2− 2ρst + t2 ν(1− ρ2) } −(ν+2)/2dsdv. (9)
It has two parameters: the degree of freedom ν and the correlation coefficient ρ. Here t−1ν is the
marginal quantile function of the univariate T -distribution. It contains the normal copula as its limit and accommodates a wide range of tail-thickness and tail-dependence. The T-copula has been shown to be a popular parametric model in empirical finance applications as reported in Embrechts, Lindskog & McNeil (2003).
5. SIMULATION STUDIES
We report results from simulation studies which are designed to confirm the theoretical findings in Section 3 and the proposed bandwidth selection method in Section 4. To demonstrate the advantage of kernel smoothing, the kernel estimator is compared with the unsmoothed estimator
˜
C given in (7).
Three copulas are considered in the simulation study, which are respectively
C(u, v) = uv
1− θ(1 − u)(1 − v), (10)
the Ali-Mikhail-Haq (AMH) family with θ = 1;
C(u, v) = exp(−[(− log u)θ+ (− log v)θ]1/θ) (11)
the Gumbel copula with θ = 2; and C(u, v) = Φ−1(u) −∞ Φ−1(v) −∞ 1 2π(1− ρ2)1/2exp{− s2− 2ρst + t2 2(1− ρ2) }dsdt
the normal copula where Φ is the standard normal distribution function and ρ is the correlation coefficient and was set at 0.5 in the simulation.
For each copula model, we first generate independent and identically distributed uniform ran-dom variables{Ui}n
i=1. Then, generate Vifrom the conditional copula distribution given Ui, which
is known under each model. The sample sizes considered are n = 50 and 100. We choose b1= b2= b
as the marginals have already been standardized.
The first simulation study is designed to check on the effect of smoothing at each of the two stages. Twenty equally spaced bandwidths are chosen for b and h respectively. For each given pair (b, h), the MISE and mean integrated variance (MIV) of the kernel and the unsmoothed
copula estimators are evaluated over 40×40 equally spaced grid points within [0, 1]2based on 1000
simulations. We present in Figure 1 only the MISEs and MIVs for the Gumbel copula as those for the other two follow the same pattern.
The results conveyed by Figures 1 can be summarized as follows. First of all, the smoothing at the first stage has little effect on the variance of the kernel estimator for all the three copulas and sample sizes considered. Indeed, this is shown for sample size as small as 50, which is the minimum sample size tried in the simulation. In particular, the shapes of the MISE and MIV contours coincide with our early predictions that (i) the role of first stage smoothing is in the bias
and has little affect on the variance as ρ∗ is slow varying and (ii) variance reduction is largely
due to the second stage smoothing. The simulation also shows that kernel smoothing leads to a substantial improvement in estimation accuracy as compared with the unsmoothed estimator. Indeed, for n = 50 and each of the copula models considered, the MISE of the kernel estimator is
b h 0.1 0.2 0.3 0.4 0.1 0.2 0.3 0.4 MISE,N=50 b h 0.1 0.2 0.3 0.4 0.1 0.2 0.3 0.4 VARIANCE,N=50 b h 0.05 0.15 0.25 0.05 0.15 0.25 MISE,N=100 b h 0.05 0.15 0.25 0.05 0.15 0.25 VARIANCE,N=100
Figure 1: MISE and Integrated Variance (IVAR) of kernel estimator for Gumbel copula. The MISE and IVAR for the unsmoothed estimator (7) are 0.000649 and 0.000309 for n = 50; and 0.000427 and 0.000292 for n = 100, respectively.
nearly half of the unsmoothed estimator. Although the gap between the two estimators is reduced, there is still around 30% advantage for the kernel estimator when the sample size is 100.
To evaluate the practical performance of the proposed reference to the T -copula rule for selecting h-bandwidth, we conducted simulations for the same three copula models to obtain the MISEs of the kernel copula estimator using (i) the prescribed reference-rule and (ii) a set of fixed bandwidths,
respectively, while setting b1 = b2 = 10−4 to realize the strategy of undersmoothing in the first
stage. The results of the simulation, displayed in Figure 2, show that the reference rule is able to achieve a level of MISE which is consistently close to the minimum MISE of the fixed bandwidth estimator. This is particularly encouraging as all the three copulas are not the T -copula and indicates that the proposed rule is robust against mis-specifying the copula model in bandwidth selection. We also used the marginal empirical distributions instead of smoothing in the first stage,
effectively setting b1 = b2 = 0. The performance of the copula estimator was almost identical to
the first stage undersmoothed estimator (with b1= b2= 10−4) in the interior of [0, 1]2. However,
it had some noticeable increase of bias near the upper and rights edges of the unit square, as compared to the undersmoothed estimator. This indicates some finite sample benefits by carrying out some smoothing in the first stage.
6. EMPIRICAL STUDY
We carry out an empirical study on a set of Uranium exploration data collected from water samples in the Montrose quadrangle in Colorado, which was originally studied in Cook & Johnson (1981).
The same dataset has been analyzed by Genest & Rivest (1993) and Genest, Quessy & R´emillard
(2006) to demonstrate a semiparametric inference for Archimedean copulas and a goodness-of-fit test. The dataset contains 655 log-concentrations of seven chemical elements including Uranium, Caesium and Lithium. A primary interest is to understand the dependence in concentrations
0.2 0.3 0.4 0.5 0.00038 0.00042 h MISE Fixed h Reference rule AMH, N=50 0.15 0.25 0.35 0.00028 0.00034 h MISE Fixed h Reference rule AMH, N=100 0.1 0.2 0.3 0.4 0.00036 0.00042 h MISE Fixed h Reference rule Gumbel, N=50 0.05 0.15 0.25 0.00026 0.00034 h MISE Fixed h Reference rule Gumbel, N=100 0.2 0.3 0.4 0.5 0.00047 0.00052 h MISE Fixed h Reference rule Normal, N=50 0.15 0.25 0.35 0.00030 0.00036 h MISE Fixed h Reference rule Normal, N=100
Figure 2: MISE of the kernel estimators with the plug-in bandwidth (dashed lines) and fixed bandwidths (solid lines).
between an actinide metal Uranium and two alkali metals, Caesium and Lithium.
Figure 3 displays the original data in panel (a) for Uranium versus Cesium and in panel (c) for Uranium versus Lithium. The kernel copula estimators are displayed in panels (b) and (d) with the h-bandwidth chosen by the proposed reference rule which assigns h = 0.176 for Uranium
versus Caesium, and h = 0.143 for Uranium versus Lithium, whereas b1= b2= 10−4.
The objective of the empirical study is to find a copula model for the two pairs of chemical elements which best describes the underlying dependence structure. We considered four parametric copulas which are respectively the AMH copula (10), the Gumbel copula (11), the Clayton copula Cθ(u, v) = max{
u−θ+ v−θ− 1 −1/θ, 0} which was used in Cook and Johnson (1981)’s original study and the T -copula (9).
The parameter of each copula needs to be estimated before we can check on the adequacy of each copula model. The first three copulas are members of the Archimedean family (Nelson,
1998), which can be expressed as C(u, v) = φ−1{φ(u) + φ(v)} for a convex decreasing function
φ (the generator) such that φ(1) = 0. The generator φ(t) is log{1−θ(1−t)t } for the AMH copula, {− log(t)}θ for the Gumbel copula and (t−θ− 1)/θ for the Clayton copula. We use the method of
moment estimator proposed in Genest & Rivest (1993) which is based on the following equation regarding Kendall’s tau
τ (X1, X2) = 4 C(u, v)dC(u, v)− 1 = 4 1 0 φθ(u) φθ(u)du + 1
after replacing τ (X1, X2) by its sample version. The parameters of the T -copula are estimated by
the method of moments too. Specifically, ρ is estimated by the sample correlation coefficient and ˆ
ν = max{ˆν1, ˆν2} where ˆνi = (4m4i+ 6m2i2)(m4i− 3m22i)−1 and mki = n−1
n
j=1(Xji− ¯Xi)k for
i = 1 and 2 and k = 2 and 4. Here ˆνi is the method of moment estimator of ν based on the i-th
0.5 1.0 1.5 2.0 1.5 2.0 2.5 (a) Uranium Caesium (b) Uranium Caesium 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.4 0.8 0.5 1.0 1.5 2.0 1.0 1.5 2.0 (c) Uranium Lithium (d) Uranium Lithium 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.4 0.8
Figure 3: Log-concentrations of Uranium versus Caesium in (a) and of Uranium versus Lithium in (c); the kernel copula estimators with the plug-in bandwidth (solid lines) of Uranium versus Caesium in (c) and of Uranium versus Lithium in (d).
Let ˆθ be the method of moments estimator and Cθˆbe the estimated parametric copula function.
Figures 4 and 5 displays the four parametric copulas at ˆθ and the kernel copula estimate for
Uranium versus Lithium. Copulas are monotone non-decreasing with respect to each variable and
in particular the contour curves {(u, v)|C(u, v) = t} are all confined in a triangle with vertices
(t, t), (1, t) and (t, 1). These features make the copula estimates look rather similar to each other. To check on the goodness-of-fit of a parametric copula model, a formal test procedure is needed as visual diagnostics can hardly detect the differences.
Let Cθˆbe the estimated parametric copula model with a parameter estimate ˆθ, and ˆC be the
kernel estimator based on a smoothing bandwidth h. We propose the following Cram´er-Von Mises
type test statistic
Tn = n 1 0 1 0 { ˆC(u, v)− Cθˆ(u, v)} 2 dudv = n−1 1 0 1 0 n i,j {Gu,h u− ˆF1(Xi1) h Gv,h v− ˆF2(Xi2) h − Cθˆ(u, v)− B(u, v)} × ×{Gu,h u− ˆF1(Xj1) h Gv,h v− ˆF2(Xj2) h
− Cθˆ(u, v)− B(u, v)}dudv
where B(u, v) = (uTv,h+ vTu,h+ Tu,hTv,h). Clearly, Tn is a L2-distance between the kernel
Uranium Caesium 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.4 0.8 AMH Uranium Caesium 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.4 0.8 Gumbel Uranium Caesium 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.4 0.8 Clayton Uranium Caesium 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.4 0.8 T−Copula
Figure 4: Copulas implied by the parametric models (dashed lines) and the kernel estimator (solid lines) for Uranium versus Caesium.
Let cαbe the upper-α quantile of the test statistic Tnat a level of significance α. The following
bootstrap procedure is employed to obtain an estimate of cα.
Step 1: Generate {Xi1∗}n
i=1 from Fn1, the empirical distribution of {Xn1}ni=1 by sampling
with replacement, and let Ui∗= Fn1(Xi1∗) for i = 1, . . . , n.
Step 2: Generate Vi∗from Cθ,Uˆ ∗ i(v) =
∂Cθˆ(u,v)
∂u |u=Ui∗, the conditional distribution of V given
U = Ui∗, and let Xi2∗ = Fn2−1(Vi∗) where Fn2 is the empirical distribution of{Xi2}ni=1. Then
{(X∗
i1, Xi2∗}ni=1 constitutes a bootstrap resample which respects both the parametric copula
and the two marginals.
Step 3: Construct ˆC∗(u, v), the kernel estimator based on the bootstrap resample using the
same h as in Tn and let Tn∗= n0101{ ˆC(u, v)− Cθˆ∗(u, v)}2dudv where ˆθ∗ is the parameter estimate based on the resample.
Step 4: Repeat the above steps B times for a large integer B and obtain, without loss of generality, Tn1∗ ≤ · · · ≤ TnB∗ . Estimate cα by Tn[B(1−α)]+1∗ , the upper α-th order statistic.
The proposed goodness-of-fit test rejects the parametric copula model as conveyed in H0 if
Tn≥ ˆcα. As the study of the kernel copula estimation has occupied large amount of space in the
present paper, the theoretical properties of the proposed goodness-of-fit test will be reported in a future paper. The purpose of presenting the goodness-of-fit test is to demonstrate the usefulness of the kernel copula estimator.
We apply the above procedure to test for the four copulas for Uranium versus Caesium and Uranium versus Lithium respectively, and the results are summarized in Table 1. It is found that both AMH and Clayton copulas are overwhelmingly rejected for both pairs despite that Clayton copula was the one used in Cook & Johnson (1981)’s original study. Gumbel copula is rejected
Uranium Lithium 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.4 0.8 AMH Uranium Lithium 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.4 0.8 Gumbel Uranium Lithium 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.4 0.8 Clayton Uranium Lithium 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.4 0.8 T−Copula
Figure 5: Copulas implied by the parametric models (dashed lines) and the kernel estimator (solid lines) for Uranium versus Lithium.
for Uranium versus Caesium but not for Uranium versus Lithium. The T-copula seems to provide the best dependence description for both pairs of data especially for the pair of Uranium versus Lithium. The goodness-of-fit offered by the T-copula echoes promising results in empirical finance (Embrechts, Lindskog & McNeil, 2003), which has been shown to be robust in fitting financial data, which typically have heavy tails and tail dependence.
APPENDIX
Proof of Proposition 1. We first outline some key steps in the proof of Proposition 1. A proof of (2) and a summary of the main results in the calculation of (3) are followed. Complete derivations can be found in a technical report (Chen and Huang, 2005).
For i = 1, ..., n, let Δ1,i= F1(Xi1)− ˆF1(Xi1), Δ2,i= F2(Xi2)− ˆF2(Xi2), Ij,k(s, t) = G(j)u,h{(u −
s)/h}G(k)v,h{(v − t)/h} for j, k ≥ 0 and S = {(j, k) : j ≥ 0, k ≥ 0, j + k ≤ 5}. A Taylor expansion for ˆC(u, v) is ˆ C(u, v) = (k,j)∈S Ak,j+ Rn− (vTu,h+ uTv,h+ Tu,hTv,h), (12) where Aj,k= n−1 n i=1 j!k!1 Ij,k(F1(Xi1), F2(Xi2)) (Δ1,i/h) j
(Δ2,i/h)k and, for some θ∈ [0, 1],
Rn= 6 j=0 1 (n)j!(6− j)! n i=1
Ij,6−j{F1(Xi1)− θΔ1,i, F2(Xi2)− θΔ2,i} (Δ1,i/h)j(Δ2,i/h)6−j.
Derivations given in the technical report show that E(Rn2) is o(h2/n), which means that Rn is
Table 1: Testing results for the four copula models (a) Uranium versus Caesium
Model Test Statistic 5% critical value p-value parameter estimate
AMH 0.360 0.0254 < 0.001 θ = 1
Gumbel 0.0484 0.0194 < 0.001 θ = 1.88
Clayton 0.173 0.0215 < 0.001 θ = 1.76
T 0.065 0.107 0.283 ρ = .60, m = 59
(b) Uranium versus Lithium
Model Test Statistic 5% critical value p-value parameter estimate
AMH 0.1334 0.0689 < 0.001 θ = 0.7675
Gumbel 0.0137 0.0210 0.221 θ = 1.1512
Clayton 0.0338 0.0179 < 0.001 θ = 0.3024
T 0.0212 0.0549 0.605 ρ = 0.17, m = 59
can be shown that E(Aj,k) = o(h2) for j + k≥ 2 and
E(A0,0) = C(u, v) +12h2{Cuu(u, v)ν(u, h) + Cvv(u, v)ν(v, h)}
+vTu,h+ uTv,h+ Tu,hTv,h+ o(h2),
E(A1,0) = −12σ2Kb21{Cu(u, v) + Tv,h}f1(1){F1−1(u)} + o(h2),
E(A0,1) = −12σ2Kb22{Cv(u, v) + Tu,h}f2(1){F2−1(v)} + o(h2).
These lead to (2). To establish (3), we note that
Cov(Aj,k, A˜j,˜k) = o(h4) if j + k + ˜j + ˜k > 2.
Therefore,
V ar{ ˆC(u, v)} = V ar(A0,0) + 2Cov(A0,0, A0,1) + 2Cov(A0,0, A1,0)
+ V ar(A0,1) + V ar(A1,0) + 2Cov(A0,1, A1,0) + o(h4). (13) The variance expression is attained after deriving each covariance term above.
The proof of Proposition 1 relies on the Taylor expansion of ˆC(u, v) in (12). To make the
expansion legitimate, we assume in this section that sup
x |K
(j)(x)| < ∞ for j = 0, 1, ..., 6. (14)
However, (14) can be removed based on the following arguments. Write ˆCW(u, v) for the copula
estimator based on a kernel W , and Lbias,W(u, v) and Lvar,W(u, v) be the leading order terms of
the bias and variance of ˆCW(u, v) respectively so that E{ ˆCW(u, v)} = Lbias,W(u, v) + o(h2) and
V ar{ ˆCW(u, v)} = Lvar,W(u, v) + o(h/n). Then, for any kernel K satisfying Condition A1, there is
a symmetric bounded kernel K supported on [−1, 1] that satisfies (14) and approximates K well
enough in L1 norm so that
sup
u,v,n,b1,b2,h|Lbias,K
(u, v)− Lbias,K∗(u, v)| + |Lvar,K(u, v)− Lvar,K∗(u, v)|
+|E ˆCK(u, v)− E ˆCK∗(u, v)| + |V ar{ ˆCK(u, v)} − V ar{ ˆCK∗(u, v)}| = o(n−2).
Thus (2) and (3) hold for K if they hold true for K. Therefore, it is sufficient to prove (2) and
Below we will give a detailed proof for (2) and a summary of key results in the proof for (3). We first state a lemma that gives the orders of the Δ1,i’s and Δ2,i’s.
Lemma 1. For nonnegative integers a1, b1, a2, b2,
sup x1,y1,x2,y2 EΔa1 1,1Δb2,11 Δa1,22Δb1,22 (X11, X12, X21, X22) = (x1, y1, x2, y2) = O n−(a1+b1+a2+b2)/2.
Proof of Lemma 1. Note that
E Δa1 1,1Δb2,11 Δa1,22Δb1,22 (X11, X12, X21, X22) = (x1, y1, x2, y2) = E 2 j=1 Aj− n−1 n i=3 Uj,i aj Bj− n−1 n i=3 Vj,i bj , where Aj = F1(xj)− EG xj− X11 b1 + 2n−1EG xj− X11 b1 − n−1G(0)− n−1Gxj− x3−j b1 , Uj,i = G xj− Xi1 b1 − EG xj− X11 b1 , Bj = F2(yj)− EG yj− X12 b2 + 2n−1EG yj− X12 b2 − n−1G(0)− n−1Gyj− y3−j b2 , and Vj,i = G yj− Xi2 b2 − EG yj− X12 b2 . Since EG x−X11 b1
− F1(x) =F1(x− sb1)K(s)ds− F1(x), the first derivative of f1 is bounded,
and b1= O(n−1/3), we have
sup
x1,y1,x2,y2
A4aj j = O(b8a1 j) = O n−2aj and for some positive constants C1and C2,
sup x1,y1,x2,y2 E Aj− n−1 n i=3 Uj,i 4aj ≤ C1 sup x1,y1,x2,y2 A4aj j + C2 ⎛ ⎝ sup x1,y1,x2,y2 En−1 n i=3 Uj,i 4aj⎞ ⎠ = O n−2aj . Similarly, supx1,y1,x2,y2E Bj− n−1 n i=3Vj,i 4bj
= O n−2bj . From Cauchy Schwartz Inequality
and the above calculation, sup x1,y1,x2,y2 EΔa1 1,1Δb2,11 Δa1,22Δb1,22 (X11, X12, X21, X22) = (x1, y1, x2, y2) ≤ ⎛ ⎝2 j=1 O n−2aj O n−2bj ⎞ ⎠ 1/4 = O n−(a1+b1+a2+b2)/2.
From (12), E{ ˆC(u, v)− Rn} = (j,k)∈SE(Aj,k)− vTu,h− uTv,h− Tu,hTv,h. From Lemma 1 and the fact that
sup 0<h<0.5h −1 1 0 G(m) u,h u− s h
ds < ∞ for u ∈ [0,1] and for 1 ≤ m ≤ 5. (15)
E(Aj,k) = o(h2) for j + k > 2 or (j, k) = (1, 1). Thus we will only compute E(Aj,k) for other cases to complete the proof of (2), which requires the following lemma.
Lemma 2. Suppose that G1 and H1 are absolutely continuous functions on (−∞, ∞) and g is a
continuous function on [0, 1]2. Let C∗(u, v) =0u0vg(s, t)f (s, t)dtds, then for u, v in [0, 1], (i) h−1 G1 u− F1(x) h H1 v− F2(y) h g(F1(x), F2(y))dF (x, y) = u h u−1 h v h v−1 h Cv∗(u− sh, v − th)dH1(t)dG1(s) + u−1 h −∞ v h v−1 h Cv∗(1, v− th)dH1(t)dG1(s) and
(ii) if0ug(s, 0)f (s, 0)ds = 0 =0ug(s, 1)f (s, 1)ds and H1 and H1 have compact supports, then
G1 u−F1(x) h H1 v−F2(y) h g(F1(x), F2(y))dF (x, y) = o(h).
Proof of (15) and Lemma 2. First of all, (15) can be derived by integration by parts and a fact that sup h≤0.5,0≤c≤1 a2(c, h) +|a1(c, h)| |a0(c, h)a2(c, h)− a21(c, h)| <∞. (16)
The argument for (16) is as follows. By the symmetry of K, a0(c, h) and a2(c, h) are minimized at
c = 0, 1 and maximized at c = 0.5, and|a1(c, h)| is maximized at c = 0, 1. Therefore,
a2(c, h) +|a1(c, h)|
|a0(c, h)a2(c, h)− a21(c, h)| ≤
σK2 + σK
a0(1, 0.5)a2(1, 0.5)− a21(1, 0.5)
and we have (16).
Lemma 2 is obtained by taking the partial derivative with respect to v of both sides of the following equation: G1 u− F1(x) h H1 v− F2(y) h g(F1(x), F2(y))dF (x, y) = u h u−1 h v h v−1 h C∗(u− sh, v − th)dH1(t)dG1(s) + (u−1)/h −∞ v h v−1 h C∗(1, v− th)dH1(t)dG1(s) + u h u−1 h v−1 h −∞ C∗(u− sh, 1)dH1(t)dG1(s) + (u−1)/h −∞ v−1 h −∞ C∗(1, 1)dH1(t)dG1(s).
We are now ready to finish the proof of (2) by computing E(Aj,k) for j+k≤ 2 and (j, k) = (1, 1). For (j, k) = (0, 2),
2E(A0,2) = h−2
where E(Δ22,1|X11= x, X12= y) = E F2(y)− EG y− X12 b2 + n−1EG y− X12 b2 − n−1G(0)− n−1n i=2 Vi 2 and Vi = G y− Xi2 b2 − EG y− X12 b2 . Since F2(y)− EG y− X12 b2 = F2(y)− F2(y− sb2)K(s)ds = −σ 2 Kb22f (1) 2 (y) 2 + o(b 3 2), 2E(A0,2) = o(h4) +σ 4 Kb42 4h2 I0,2(F1(x), F2(y)) f2(1)(y) 2 dF (x, y) +σ 2 Kb 2 2 nh2
I0,2(F1(x), F2(y))f2(1)(y)(G(0)− F2(y))dF (x, y)
+ (nh2)−1 I0,2(F1(x), F2(y))V ar G y− X12 b2 dF (x, y).(17) From (15) and the fact that
V ar G y− X12 b2 = F2(y− sb2)dG2(s)− F2(y− sb2)K(s)ds 2 ,
(17) gives 2E(A0,2) = o(h2) + (nh2)−1I0,2(F1(x), F2(y))F2(y)(1− F2(y))dF (x, y), which is o(h2)
by Lemma 2 part (ii).
Similar derivations yield E(A2,0) = o(h2), and
E(A0,0) = C(u, v) +12h2{Cuu(u, v)ν(u, h) + Cvv(u, v)ν(v, h)}
+vTu,h+ uTv,h+ Tu,hTv,h+ o(h2),
E(A1,0) = −12σ2Kb21{Cu(u, v) + Tv,h}f1(1){F1−1(u)} + o(h2),
E(A0,1) = −12σ2Kb22{Cv(u, v) + Tu,h}f2(1){F2−1(v)} + o(h2).
These lead to (2) since
E{ ˆC(u, v)} = E(A0,0) + E(A0,1) + E(A1,0)− (vTu,h+ uTv,h+ Tu,hTv,h) + o(h2). Below we summarize the main results in the proof of (3). It can be shown that
V ar(A0,1) = n−1v(1− v){Cv(u, v) + Tu,h}2− hn−1{Cv(u, v) + Tu,h}2μ2(v, h, b2/h) + o(h4),
V ar(A1,0) = n−1u(1− u){Cu(u, v) + Tv,h}2− hn−1{Cu(u, v) + Tv,h}2μ1(u, h, b1/h) + o(h4);
Cov(A0,0, A0,1) = −n−1{Cv(u, v) + Tu,h} (1− v)C(u, v) + v(1 − v)Tu,h
+ {C(u, v) − uv}Tv,h+ hn−1{Cv(u, v) + Tu,h}2μ∗2(v, h, b2/h) + o(h4),
Cov(A0,0, A1,0) = −n−1{Cu(u, v) + Tv,h} (1− u)C(u, v) + u(1 − u)Tv,h + {C(u, v) − uv}Tu,h
+ hn−1{Cu(u, v) + Tv,h}2μ∗1(u, h, b1/h) + o(h4); Cov(A1,0, A0,1) = n−1{C(u, v) − uv}{Cv(u, v) + Tu,h}{Cu(u, v) + Tv,h} + o(h4);
and
V ar(A0,0) = n−1C(u, v)(1 + 2Tu,h)(1 + 2Tv,h)
+ n−1{Tu,h2 v(1 + 2Tv,h) + Tv,h2 u(1 + 2Tu,h) + Tu,h2 T 2 v,h} − n−1h{C u(u, v)(1 + 2Tv,h) + Tv,h2 } u h u−1 h sdG2u,h(s) − n−1h{C
v(u, v)(1 + 2Tu,h) + Tu,h2 }
v
h v−1
h
tdG2v,h(t)− n−1{E(A0,0)}2+ o(h4).
Substituting these into (13), (3) can be established. ACKNOWLEDGEMENT
We thank two referees and the editor for constructive comments and suggestions that improve the presentation of the paper, and Cheng Yong Tang for valuable computation assistance and Professor Christian Genest for making the Uranium exploration data available. We also acknowledge the support of an Iowa State University Start-up grant.
REFERENCES
A. Bowman, P. Hall & T. Prvan (1998). Bandwidth selection for the smoothing of distribution functions.
Biometrika, 85, 799–808.
R. D. Cook & M. E. Johnson (1981). A family of distributions for modelling nonelliptically symmetric multivariate data. Journal of the Royal Statistical Society Series B, 43, 210–218.
Chen, S.X. and Huang, T. (2005), Bias and variance calculation for a copula kernel estimator, Technical Report 2005-08, Department of Statistics, Iowa State University.
S. X. Chen & Y. Fan (2006). Estimation of copula-based semiparametric time series models. Journal of
Econometrics, 130, 307–335.
D. G. Clayton (1978). A model for association in bivariate life tables and its applications in epidemio-logical studies of familial tendency in chronic disease incident. Biometrika, 65, 141–151.
P. Deheuvels (1979). La fonction de d´ependance empirique et ses propri´et´es. Un test non param´erique
d’ind´ependance. Acad. Roy. Belg. Bull. Cl. Sci., 65, 274–292.
P. Embrechts, F. Lindskog, & A. McNeil (2003). Modelling dependence with copulas and applications to risk management. In Handbook of Heavy Tailed Distributions in Finance, ed. S. Rachev, Elsevier, 329–384.
J. D. Fermanian (2005). Goodness-of-fit tests for copulas. Journal of Multivariate Analysis, 95, 119–152. J. Fermanian, D. Radulovic & M. Wegkamp (2004). Weak convergence of empirical copula processes.
Bernoulli, 10, 847–860.
J. D. Fermanian & O. Scaillet (2003). Nonparametric estimation of copulas for time series. Journal of
Risk, 5, 25–54.
J. P. Fine & H. Jiang (2000). On association in a copula with time transformations. Biometrika, 87, 559–571.
C. Genest, K. Ghoudi & L.-P. Rivest (1995). A semiparametric estimation procedure of dependence parameters in multivariate families of distributions. Biometrika, 82, 543–552.
C. Genest, J.-F. Quessy & B. R´emillard (2006). Goodness-of-fit procedures for copula models based on
the probability integral transformation. Scandinavian Journal of Statistics, 33, 337–366.
C. Genest & L.-P. Rivest (1993). Statistical inference procedures for bivariate Archimedean copulas.
Journal of American Statistical Association, 88, 1034–1043.
I. Gijbels & J. Mielniczuk (1990). Estimating the density of a copula function. Comm. Statist. Theory
Methods, 19, 445–464.
H. Joe (1997). Multivariate Models and Dependence Concepts, Chapman and Hall, London.
M. Lejeune & P. Sarda (1992). Smooth estimators of distribution and density functions. Computational
Statistics & Data Analysis, 14, 457-471.
M. C. Jones (1993). Simple boundary correction for kernel density estimation. Statistic Computing, 3, 135–146.
H. G. M¨uller & U. Stadtm¨uller (1999). Multivariate boundary kernels and a continuous least squares
principle. Journal of the Royal Statistical Society Series B, 61, 439–458. R. B. Nelsen (1997). An Introduction to Copulas, Springer, New York.
D. Oakes (1982). A model for association in bivariate survival data. Journal of the Royal Statistical
Society Series B, 44, 414–422.
B. Schweizer & A. Sklar (1983). Probabilistic metric spaces, North-Holland Series in Probability and Applied Mathematics, North-Holland Publishing, New York.
A. Sklar (1959). Fonctions de r´epartition ´a n dimensions et leurs marges. Publ. Inst. Statist. Univ.
Paris, 8, 229–231.
B. W. Silverman (1986). Density Estimation for Statistics and Data Analysis, Chapman and Hall, London.
Received ??? Song Xi CHEN: songchen@iastate.edu
Accepted ??? Department of Statistics, Iowa State University
Ames, IA USA, 50011-1210 Tzee-Ming HUANG: tmhuang@nccu.edu.tw Department of Statistics, National Chengchi University Taipei Taiwan-ROC, 11605
