7.2.1 Goodness-of-Fit Tests

Traditional goodness-of-fit approach is about testing the null hypothesis H₀: X ∼ F₀(x) versus H₁: otherwise. In general setting, we define a nonnegative “distance” (which may or may not be a proper mathematical distance or metric, satisfying conditions of identity of indiscernibles, symmetry, and triangle inequality) d(F*_n(x), F₀(x)) between the theoretical c.d.f. F₀(x) suggested by the model and the empirical c.d.f.

(7.1)

of the data sample X = (X₁, …, X_n), where I_A(x) is the standard indicator function of the set A. Notice the slightly nonstandard use of the denominator n + 1 in (7.1) instead of more common n, in which we follow Genest and Rivest [19] and others, who find it convenient that the empirical c.d.f. never takes the value 1. We will briefly discuss two common measures of goodness-of-fit, which are most popular in application to copula structures. These measures are Kolmogorov–Smirnov distance

(7.2)

and Cramer–von Mises distance

(7.3)

defining two corresponding criteria of closeness.

The tests based on these two distances with rejection regions of the form d_KS > c and ω² > c are convenient because the distributions of the test statistics and thus the p-values of the tested data are available under the assumptions of H₀ being true, at least for one-dimensional case. In order to determine d_KS one has to compare the vertical distance between two functions (empirical and theoretical), only at the sample points, which represent the jumps of empirical c.d.f., so that the search for the maximum is simplified. Similarly, instead of integral ω² the following test statistic can be calculated:

where x_(i) is the ith order statistic.

One way both of these measures can be applied beyond testing the null hypothesis is to compare several theoretical candidate models, selecting those with the lower values of d_KS or ω².

Let us notice that while goodness-of-fit tests are essentially non-Bayesian, both Kolmogorov–Smirnov and Cramer–von Mises distances may be used in a Bayesian context as we will see later on. We should also be aware that in case of copula models F₀(x) is a joint distribution of vector X with two or more components, which makes practical application of goodness-of-fit tests much more technically difficult. In particular, distributions of test statistics are no longer distribution free and finding the maximum for Kolmogorov–Smirnov statistic is not as trivial as in one dimension. That is why the determination of even approximate p-values of the tests requires a substantial amount of work such as Monte Carlo generation of multiple samples from the theoretical distribution [18]. However, this is not an issue for model comparison purposes.

The real issue is: using goodness-of-fit criteria for model selection, we do not put any barriers against overfitting. Which means that the theoretical model with a c.d.f. close to the empirical c.d.f. is always preferable. It is always easy to imagine a theoretical model with way too many parameters (overparameterized), which will provide a really good fit for the empirical c.d.f. However this model will take all observed features of the sample for granted, even those ”features” which have completely random origin. At the same time, it will not allow for the features not observed. It is like looking at a picture of a human face which has a wart on the nose, an ink stain in the middle of the forehead, and one ear cut-off. Should we take all other human faces to have warts and ink stains in appropriate places, but just one ear? Such a model will provide a good description of a particular portrait, but a poor prediction of what most humans look like.

7.2.2 Posterior Predictive p-values

Before moving on to the measures of model fit, which will address overfitting, let us introduce a fully Bayesian concept which may serve the purpose of testing goodness-of-fit hypotheses without applying Bayesian factors as was done in Chapter 2. To be more precise, we will follow the idea of Bayesian predictive checking introduced by Rubin [37].

Let us review a parametric version of the general goodness-of-fit test and the frequentist concept of p-value, which is consistently criticized by Bayesian statisticians and lately (along with confidence intervals) became a sore spot of applied statistics as one of the most intuitively misunderstood and misused traditional statistical techniques. This critique was epitomized in an editorial by Trafimov and Marks [44] in a popular applied journal Basic and Applied Social Psychology and further developed by Wasserstein and Lazar [47].

For simplicity, let us consider only the continuous case treating all data as one vector and parameter as another. Let us denote by X the data (it could be a finite i.i.d. sample or a different object taking values x on the sample space S, X = x ∈ S) with the distribution density (likelihood) f(x; θ) depending on parameter θ. Consider a null hypothesis H₀: θ = θ₀ to be tested against the composite alternative H₁: θ ≠ θ₀. For test statistic T(X) and critical region of the form R_T(x) = {x ∈ S: T(x) ≥ c}, define the test:

“Reject H₀” if T(x) ≥ c or “Do not reject H₀” if T(x) < c.

Significance level approach suggests determining first the size or significance level of the test as the tolerable level of the Type I Error α = P(T(X) ≥ c|θ₀), so that the critical value c(α) depends on α. However a more popular approach is to determine the p-value of the data, p = P(T(Y) ≥ T(x)|x, θ₀), which may be interpreted as the tail probability of statistic T(Y) for the given x and θ₀. We also may rewrite it as

(7.4)

Notice that under the null hypothesis p-value as a function of x is uniformly distributed on [0, 1], so that exactly 5% of all data from S(θ₀) will have p ≤ 0.05. General recommendation is to reject the null hypothesis for data with low p-values.

If the choice of statistic T(X) is smart, this test will have high power defined as rejection probability under the alternative: k = P(T(X) ≥ c|H₁). Famous Lindley’s paradox [30] guarantees that for large enough sample sizes, virtually all values θ ≠ θ₀ will be eventually rejected. However high rejection probability is not the goal of Bayesian hypothesis testing. If the model is false (and all models are false, according to George Box), but fits our data well, we do not want to reject it. This is the motivation to introduce a Bayesian concept, which will measure the feasibility not of a single value of parameter θ₀, but rather of the entire posterior distribution, see also Meng [31].

In order to formulate a Bayesian alternative to the p-value approach, we need to introduce Bayesian predictive distributions, which we ignored so far being focused on estimation. In estimation framework we have considered prior π(θ) and posterior π(θ|x) densities. Now let us consider a new sample, y ∈ S. What is the distribution of y given our information on θ, provided both by the prior and the data x? We can define Bayesian predictive density as the integral over posterior

(7.5)

It is also possible to define prior predictive density as an extreme case with no data available (with prior replacing posterior):

It is easy to see that the prior predictive density g(y) coincides with the marginal density m(y) in the denominator of the formula (2.26) playing the central role in Bayesian estimation.

Let us use the posterior predictive density to define the posterior predictive p-value similar to (7.4) as

(7.6)

It is important that the latter does not depend on a particular value of θ₀, but integrates the information from the posterior. The role of pp(x) is similar to traditional p-values: reject null when pp(x) is small; let us say, if it does not exceed usual benchmark of 0.1, 0.05, or 0.01.

In order to illustrate the difference between (7.4) and (7.6), let us consider a slightly paraphrased example introduced by Gelman [15].

Example 7.2.1 Let us test a hypothesis regarding the normal mean using a single data point: H₀: θ = 0 versus a two-sided alternative H₁: θ ≠ 0. Assume f(x; θ) ∼ N(θ, 1) and π(θ) ∼ N(0, τ²). It is easy to see (compare to Exercise 2.8) that

Defining T(X) = |X| and taking into account that one tail of the test statistic translates into two symmetric tails of the standard normal distribution Φ(u), we obtain from (7.4) and (7.6) the following expressions:

What can we say about particular data values which are extreme from the point of view of the likelihood? Consider, for instance, x = 3, three standard deviations away from the mean. It gives ample evidence against H₀ from frequentist point of view, bringing about p-value of 0.0027. However, if we consider a vague prior with τ = 10, pp(3) = 0.9832 and H₀ should hardly be rejected from a Bayesian standpoint. Even if x = 30, which is three standard deviations of prior predictive distribution away from the mean, p(30, 0) ≈ 0, but pp(30) = 0.8333 and hypothesis H₀ stays highly feasible and should not be rejected.

But why do we not want our hypothesis to be rejected when the observation x = 30 seriously contradicts the prior? Yes, in this case the prior is far away from the likelihood, but somehow the posterior is O’K: π(θ|30) ∼ N(29.7, 0.99), which means that the Bayesian model works too well to be discarded. The difference between the classical decision to reject H₀ and the Bayesian decision not to reject corresponds to the difference between the way we treat the null hypothesis: in classical setting it characterizes the point θ₀ of the parametric space, while in Bayesian setting it corresponds to the entire posterior distribution.

7.2.3 Information Criteria

In parametric setting, a special role when the theoretical c.d.f. is defined as is played by the likelihood function L(x; θ) defined in (2.13). Let us define deviance of a parametric model M with the k-dimensional parameter value θ as

(7.7)

The smaller values of the deviance correspond to the higher likelihood. The maximum likelihood estimator will provide the smallest possible deviance. We will define Akaike Information Criterion or AIC [1] as

(7.8)

Information criterion represents the relative loss of information which occurs when we use model M instead of the true model [2]. The lower values of AIC correspond to the better quality of the model. Most of it is determined by likelihood, but additional term 2k represents a penalty for overfitting: larger models with more parameters (higher dimension of vector θ) will be penalized, and smaller parsimonious models will be preferred.

Another information criterion was suggested by Schwarz [39]. It penalizes overparameterized models even more substantially. It is called Bayesian Information Criterion or BIC and is defined as

(7.9)

The names might be misleading, and there is nothing particularly Bayesian in BIC. Both AIC and BIC are based on maximizing the likelihood function, and only the sizes of the penalty for the parametric dimension are different, for BIC increasing with the sample size. The main advantage of both AIC and BIC and their later versions [2] becomes evident when we compare two or more parametric models of different parametric dimensions or even different parametrizations with regard to overfitting. Let us say, there exist two models M₁ with parameter θ₁ of dimension k₁ and M₂ with parameter θ₂ of dimension k₂. Then we say that a model M₁ is better than M₂ if AIC(M₁) < AIC(M₂) or BIC(M₁) < BIC(M₂), depending on your preference of AIC or BIC. Models M₁ and M₂ may be non-nested (e.g., Clayton copula with Weibull margins vs. Gumbel–Hougaard copula with Gompertz margins).

None of these criteria are especially convenient for comparison of models within Bayesian framework. First, they are based on the value of MLE, which might not be directly available in Bayesian estimation. Second, they do not utilize the entire form of the posterior distribution. Deviance Information Criterion introduced by Speielgelhalter et al. [43] resolves these problems. For a model M with parameter θ it may be defined as

(7.10)

where is the deviance at the posterior mean and is the posterior mean of the deviance D(θ). In this presentation the first term characterizes the fit of the model, and the second term implicitly penalizes overparameterization. Another version of DIC suggested by Gelman et al. in Reference [16]

(7.11)

is based on the posterior variance. The most attractive feature of DIC is its straightforward application to MCMC procedures of Bayesian estimation, where all three characteristics , , and can be estimated directly from the generated sample. As with other information criteria, the lower values of DIC indicate the better model.

7.2.4 Concordance Measures

What is the most important feature of pair copula models? It is the way they make it possible to model joint distribution of two variables P(X ≤ x, Y ≤ y) = C_α(F(x), G(y)), and also to separate the estimation of dependence from the estimation of marginal distributions F(x) and G(y). A better copula model should accurately estimate the strength of dependence whatever the marginals are. The strength of dependence suggested by model C_α within a given copula family is determined by the association parameter(s). The strength of dependence between components of the random vector (X, Y) can be also expressed in terms of Spearman’s coefficient of rank correlation ρ* or Kendall’s concordance τ defined in (5.4). It is important that both the rank correlation and concordance are invariant with respect to marginals F and G and are fully determined by the copula function C_α(u, v). Relationships between the copula function and its concordance measures are well established:

(7.12)

and

(7.13)

In particular, for Archimedean copula with generator ϕ_α(t) it holds that

(7.14)

For instance, see [20], for Clayton’s family

(7.15)

for Gumbel–Hougaard’s family (also, its survival copula version)

(7.16)

For Frank’s family the formula is not as simple, though we still can express concordance

(7.17)

through Debye’s integral

(7.18)

and approximate it numerically. Also, for FGM family

and for all elliptical copulas (including Gaussian and t-copulas, see [9])

(7.19)

For two-parametric BB1 copulas (6.19) concordance may be expressed through both association parameters:

Using the fact that Kendall’s tau may be estimated with the help of its empirical value according to (5.5), we can suggest two options for one-parametric copulas. First, may be used to directly obtain an association parameter estimate within a given family of copulas. Second, if several models with different values of association are considered, the difference between the values of τ induced by models via the formulae above and the empirical value may be treated as a concordance-based criterion of model comparison.

7.2.5 Tail Dependence

In applications of copula models to risk management, the attention is often concentrated on the tails of the joint distributions indicating extreme co-movements of two or more variables of interest. These tails correspond to catastrophic events. That is why it makes sense to construct the models with especially good fit at the tails. Coefficients of lower tail dependence

(7.20)

and upper tail dependence

(7.21)

where is the survival copula defined in (6.4), characterize respectively the behavior of the lower left and upper right tails of the copula function. If tail dependence coefficient is zero, we can talk about ”no tail dependence.” As well as the concordance measures from the previous subsection, tail dependence is invariant with respect to marginal distributions.

For popular elliptical and Archimedean copulas tail dependence can be expressed in terms of association parameter, namely, for Clayton’s family

(7.22)

for Gumbel–Hougaard copula

(7.23)

and for its survival version

(7.24)

Frank and Gaussian copula have no tail dependence, while t-copula has both lower and upper ones:

(7.25)

where T_ν(t) is the c.d.f. of t-distribution with ν degrees of freedom [9]. Two-parametric BB1 copulas allow for two different tail dependence coefficients:

(7.26)

A better copula model should give a good correspondence of the model-induced tail dependence with the empirical tail dependence. The only problem is caused by a certain ambiguity in estimating empirical tail dependence.

7.3.1 Parametric, Semiparametric, or Nonparametric?

If we choose the two-step IFM approach, it allows for one more choice. In general, one may choose nonparametric approach for both steps: estimating marginal parameters and building a nonparametric copula, see, for example, Genest and Rivest [19] and Nicoloutsopoulos [34]. In the book dedicated to parametric estimation we will choose not to consider a fully nonparametric approach and concentrate on parametric copulas. However it still allows for both parametric and nonparametric estimation of the marginals thus bringing about either fully parametric estimation (for both marginal and association parameters) or semiparametric estimation (nonparametric marginal and parametric copula model).

Example 7.3.2 In order to apply IFM approach to the data in Table 7.1, we have to estimate the margins first. We will start with parametric estimation. Using standard procedure for exponential distribution, we can use either MLE or the method of moments with the same results: and . Corresponding values of survival functions for sample data are listed in the first and the second rows of Table 7.2. Similarly, we can obtain empirical estimates for the same survival functions using empirical c.d.f. as defined in (7.1), see the third and the fourth rows of Table 7.2.

In general, fully parametric models will be less data dependent and will be expected to have higher predictive quality. On the other hand, semiparametric estimation is fully theoretically justifiable as discussed by Genest et al. [17] and may be more appropriate for the comparison of copula models and copula model selection, eliminating the effect of possible misspecification of the marginal distributions discussed by Kim et al. [24]. From this point on we will concentrate our attention on the copula parameter of association, estimating marginal distributions parametrically or nonparametrically, separately or simultaneously with the copula parameter as dictated by the context and convenience.

7.3.2 Method of Moments

Generalized method of moments or plug-in approach suggests (a) finding functional relationships between the parameters of interest and some distribution characteristics, such as moments; (b) using sample moments to estimate distribution moments; and (c) estimating the parameters by plugging in sample moments instead of distribution moments into the formulas established in (a). In case of copula parameters the role of moments may be played by concordance measures: Kendall’s and Spearman’s rank correlation. Formulas introduced in Subsection 7.2.4 provide functional relationships for (a) and can be used to estimate association for one-parameter Archimedean copulas and correlation for elliptical copulas. It is less obvious what moment-like measures should be used in case of two-parameter families.

Example 7.3.3 It is easy to estimate Kendall’s concordance τ for the initial sample from Table 7.2 or even from ranks in Table 7.1. Kendall’s concordance is based on ranks, so parametric (rows 1 and 2) and nonparametric (rows 3 and 4) estimates of the margins in Table 7.2 will bring about the same result . Solving Equation (7.16) for α, we obtain the method of moment estimate

(7.27)

which is pretty close to the true value. This degree of precision cannot be guaranteed for such a small sample, but confirms that the method of moments provides a “fast and dirty” procedure for the estimation of copula parameters.

7.3.3 Minimum Distance

This method suggests minimizing certain distance between the empirical distribution and distributions in the chosen family of copulas. The distribution from the family minimizing the distance is supposed to be the best choice. Difficulties with this approach may be related both to the proper choice of a distance and the technical task of finding its minimum. For certain classes of Archimedean copulas this approach may utilize the concept of Kendall’s distribution defined in Section 6.6 and was proven to be very productive by Nicoloutsopoulos [34].

To illustrate this point, let us define pseudo-observations or empirical copula for i = 1, …, n as

(7.28)

the number of sample points both strictly to the left and below of (X_i, Y_i). Slightly unusual normalizing factor n − 1 may be explained by convenience of the following representation of sample concordance in terms of pseudo-observations: . Then define for any t ∈ [0, 1] the empirical Kendall’s function

(7.29)

which can serve as an empirical analog (and estimator) of the theoretical Kendall’s function K_α(t) = P(C_α(U, V) ≤ t) corresponding to the true model. The values of Z_i and K_n(Z_i) calculated by the ranks are summarized in Table 7.3.

Then we can compare the distances

(7.30)

for different values of association. The value minimizing d_K may serve as an estimate of α.

Example 7.3.4 In the case of Gumbel–Hougaard family, according to (6.31), Kendall’s function may be represented as

(7.31)

and the numerical minimization of d_K brings about .

7.3.4 MLE and MPLE

The most popular method used in applications is probably still the method of maximum likelihood estimation (MLE). In the IFM setting [22], when on the second step we form pseudo-likelihood function by plugging in estimates of the marginals instead of their values and maximize this function to obtain the estimate for the association, it is more appropriate to talk about maximum pseudo-likelihood estimation (MPLE). In case of empirical estimates of the marginals, this method was discussed and justified by Genest et al. [18] and is also known as canonical maximum likelihood (CMLE).

Example 7.3.5 We will use both MPLE in IFM setting and full MLE to obtain copula parameter estimate for the data in Table 7.1. Pseudo-likelihood function for IFM can be obtained as

(7.32)

where c_α(u, v) is the copula density defined in (6.27); and for parametric estimation and and for semiparametric estimation. Maximizing in (7.32) brings about the estimates for the association parameter: for parametric and for semiparametric estimation (CMLE).

Full likelihood function for MLE can be expressed as

(7.33)

and maximized with respect to all three parameters at a time. The numerical results yield: , , and .

It will be presumptuous to make conclusions based on a small sample numerical example. However, it seems to be true in general that the method of moments using Kendall’s tau is the simplest and fastest way to obtain numerical value of association. Maximum likelihood approach is preferable for large samples because of its good large sample properties (asymptotic optimality). Using either nonparametric or parametric margins in IFM two-step context is possible, but numerical results may differ. The role of marginal estimation and possible marginal misspecification may be significant as our examples demonstrate. Full MLE and MPLE in parametric version typically bring about close numerical results, but it has to be decided whether the change of the marginal estimates from one-dimensional to two-dimensional case is an advantage or a drawback, and whether it is appropriate at all in the context of complex research projects with many pairs of variables being analyzed.

7.3.5 Bayesian Estimation

In Bayesian framework, we can certainly utilize both one-step and two-step approaches, and also use nonparametric estimates for the marginal. The difference between Bayesian approach and such classical methods as MLE, MPLE, or the method of moments as it was stated in Chapter 2 is more than just the choice of the best fit for our data within the given family of copulas. The benefit, as usual, is the use of entire posterior distribution instead of one representative of a family, providing more stable procedures and fully utilizing the benefits of integration versus optimization. The problem, as usual, is the choice of priors. In the presence of prior information we can elicit prior distributions using subjective considerations or use relevant data to implement empirical Bayes approach. In the absence of prior information we use noninformative priors. Examples and case studies in Sections 7.5 and 7.6 and Chapter 8 demonstrate a variety of ways in which Bayesian approach can be implemented to copula estimation. Here we will provide one particular implementation of Bayesian estimation via Metropolis algorithm using the familiar 10-point example and assuming a lack of reliable prior information.

Example 7.3.6 Let us consider a two-step IFM scheme, where both marginals are estimated first. This approach will allow us to avoid the problem of defining a joint prior for all three parameters and to assign priors independently. Assume very weak priors for the inverses of marginal parameters λ_i = 1/θ_i, i = 1, 2, π(λ₁) ∼ Gamma(ϵ, ϵ) and π(λ₂) ∼ Gamma(ϵ, ϵ), where ϵ is chosen as small as possible in order not to sacrifice computational stability. This way, the prior means are equal to 1, and prior variances are both equal to 1/ϵ, which makes priors as nonintrusive as possible. The resulting chains with 10,000 iterations produce the left and central trace plots in Figure 7.1 with Bayes estimates of the posterior mean and .

At the second stage, using the estimated marginals, assume a conservative prior , penalizing higher values of association. The resulting trace plot is shown in the bottom part of Figure 7.1 and the estimate of posterior mean is , which as can be expected for weak priors is similar to MLE.

Using the same independent priors for one-step analysis with simultaneous block update of all three parameters, we can expect a lower acceptance rate. Corresponding calculations bring about numerical estimates of , , and α = 1.747. The corresponding trace plots are demonstrated in Figure 7.2.

In this example we cannot see any distinctive advantage of using Bayesian approach. First of all, we use no prior information, and therefore the final numerical estimates are not too far from MLE or method of moments estimates. Second, there is no technical problem with appying MLE for such a simple setting, therefore Bayesian approach is not necessary. However, it is definitely applicable, and leads to reasonable numerical results.

7.4.1 Hybrid Approach

This approach will combine two methods of estimation of the copula parameter: say, method of moments and MPLE. Suppose we consider two different families of copulas and want to choose the better model. To reduce the role of marginal specification, we will use the same estimate for the marginal distributions in both cases. Let us determine MPLE for the parameter of association for each of the two families. Then we can use the idea of generalized moments or concordances. Calculate concordance values implied by MPLEs in Subsection 7.2.3 for each family. Then compare these values with empirical concordance, and the smaller deviation will determine the best fit. For instance, model-induced values of τ obtained from estimated parameter values via 7.15 and 7.16 can be compared with the empirical value , and a smaller discrepancy will indicate a better model fit.

We can also use the hybrid between the minimum distance method and MLE or MPLE. In this case we choose two best representatives of two different copula families, and then compare them from the point of view of a distance from the empirical measure. The smaller distance wins. This approach is similar to the minimum distance method from the previous section, but now the problem is much simplified by the need to compare two numbers instead of running a minimization procedure over entire family. For instance, we can calculate the values of distances d_K defined in (7.30) for two best representatives of two families and choose the one granting the smaller value. The role of such distances in hybrid approach can be also played by goodness-of-fit measures discussed in Subsection 7.2.1. However there exists some danger of choosing an overfitted model if we have to choose between two families of different parametric complexity.

7.4.2 Information Criteria

Information criteria defined in Subsection 7.2.4 penalize model complexity and overfitting, and their application makes more sense for comparison of models with a different number of parameters. These criteria are based on the values of likelihood function, so they cannot be used “across the board” and we have to compare separately two-step to two-step, IFM to IFM, or CMLE to CMLE estimates obtained for two families. These criteria are not necessarily good for nonlikelihood-based methods of estimation.

Example 7.4.1 Let us extend the example with data from Table 7.1 to a different family of copulas. Clayton’s family is often applied in the same context as survival Gumbel–Hougaard, when the modeling of the lower tails is of primary importance. Let us use Clayton copula to obtain MPLE estimates for α in parametric and nonparametric setting and also full MLE. Table 7.4 will summarize all information on numerical estimates obtained by full one-step MLE, parametric IFM (MPLE), and semiparametric method CMLE. We also calculate the model-induced values of tau and Kendall’s distance d_K to implement the hybrid approach to model selection. In order to calculate the distance in (7.30) for the Clayton’s family we can use the following formula obtained from (6.31):

(7.34)

Copula	Gumbel–Hougaard			Clayton
method	MLE	IFM	CMLE	MLE	IFM	CMLE
	2.64	2.73	*	2.72	2.73	*
	3.70	3.96	*	3.85	3.96	*
	1.69	1.69	2.03	0.86	0.86	1.48
τ	0.41	0.41	0.51	0.30	0.30	0.43
dK	0.37	0.37	0.36	0.47	0.47	0.46
Log like	−42.06	1.72	2.45	−42.51	1.48	1.98
AIC	90.12	−1.44	−2.90	91.02	−0.96	−1.96
BIC	91.03	−1.14	−2.60	91.93	−0.66	−1.66
Tail	0.49	0.49	0.59	0.45	0.45	0.63

log-likelihood function, information criteria values, and tail dependence according to Subsection 7.2.5. We will leave the derivation of Bayes estimates for Clayton’s family and comparison of two Bayes estimates for the end-of-the-chapter exercises.

It is clear that no parameter estimates are obtained for marginals when semiparametric estimation is applied. Also, we can see slight discrepancies between the values of marginal parameters obtained by full one-step MLE and MPLE with IFM. Comparing values of α between Gumbel–Hougaard and Clayton model is in general not practical, because the meaning of these parameters is quite different. In this particular case τ = 0.5 corresponds to α = 2 for both families, but this is a very rare occasion.

However, according to the hybrid approach to estimation, comparison of corresponding values of τ is always informative. The model providing concordance estimates closer to the empirical value is more adequate. We should notice though that is not a very reliable estimate of the theoretical τ, especially for small samples. For instance, in our example data were generated using Gumbel–Hougaard copula with association parameter α = 2, corresponding by (7.16) to theoretical concordance value of τ = 0.5. In this case survival Gumbel–Hougaard model performs somewhat better. This is to be expected, because the data were generated using this model, and Clayton copula should not necessarily provide a good fit.

From row 7, corresponding to d_K, it also looks like Gumbel–Hougaard family provides a better fit than Clayton: values of d_K are smaller for all pairs of estimates: MLE versus MLE, IFM versus IFM, and semiparametric CMLE versus CMLE.

Comparison of two families using AIC and BIC based on log-likelihood function should also be done pairwise. Full MLE will bring about the largest values of information criteria, which reflects the more complex structure of the model. Semiparametric model by definition will provide a better fit than the parametric IFM. We can compare though the pairs provided by similar settings for different copula choices. In all the three cases (MLE, IFM, CMLE) the Gumbel–Hougaard model ends up slightly on top. It would be interesting to compare the tail behavior of two copula families, but it is much harder to do so for the lack of reliable empirical measures of tail dependence, especially for small samples.

Bayesian estimation for Clayton family and comparison of Bayesian estimates for two families of copulas using information criterion DIC specific for Bayesian estimation is left for the end-of-chapter exercises in Chapter 8.

7.4.3 Bayesian Model Selection

Following Bretthorst [7] and Huard et al. [20], we suggest to compare the data fit provided by several pair copula models not at a single value of association parameter(s) obtained by MPLE, but rather over the entire range of possible association values. This can be accomplished by specifying a prior distribution for association parameter(s) and integrating the likelihood with respect to the prior distribution. The problem is the difference of meaning and ranges of association parameters for different copula classes. If we want to compare several classes of copulas in Bayesian framework, we need to establish the common basis of comparison. For that purpose we need to suggest a universal parameter, which can be evaluated for all classes of copulas under consideration.

One such universal parameter is Kendall’s concordance τ, which as we have seen can be expressed in terms of association for many copula families. Sample concordance is a reasonable nonparametric estimator of τ [19]. Using formulas expressing concordance through association parameters, for example, (7.15), (7.16), and (7.19), we can calculate values of τ induced by MPLE for parameters of elliptic or Archimedean copulas and compare them to the sample values of , which is done in Table 7.4. Proximity of model-induced values of τ to the sample value may serve as a measure of the model fit and help to compare the model performance as in hybrid approach to model selection discussed above in Example 7.4.1. However, this comparison is still using single values to represent entire families.

We will assume that the classes of copulas we choose correspond to exhaustive and mutually exclusive hypotheses H₁, H₂, …, H_m. Posterior probabilities of hypotheses H_k, k = 1, …, m, for data D may be rewritten as

(7.35)

where we will consider all m hypotheses a priori equally likely. If the dependence between variables is positive for all hypotheses, we can assume τ ≥ 0. In this case the natural choice of prior for τ is beta distribution, and the choice of parameters for the prior can be subjective or noninformative objective. However, in the presence of relevant additional data, similar in nature, the prior might be suggested by sample concordance for this data consistently with empirical Bayes approach: . If we have multiple pairs of variables included in the study, estimates of parameters of the beta distribution for empirical Bayes can be obtained from all pairs of components.

We will not need to calculate the denominator of the posterior in (7.35). It suffices to calculate the weights

(7.36)

or, using Monte Carlo approach and drawing samples from the Beta prior, evaluate

(7.37)

Then we select the class with the highest weight and obtain the Bayes estimate of the association parameter using MCMC.

7.6.1 Estimation of Association Parameters

The pseudo-likelihood function for the copula estimation can be written as

where and are marginal empirical distribution functions. Using a numerical implementation of the maximum likelihood method [6, 25], obtain for each class of copluas H1–H5 the CMLE (also for BB1 and and for t-copula). We present these values for each of the 10 pairs of five components in Table 7.6.

	H1	H2	H3	H4	H5	H5	H6	H6
AB	3.02	3.12	2.95	3.14	0.62	2.48	0.92	2.00
AC	0.58	1.36	0.57	1.37	0.20	1.26	0.42	2.44
AD	1.05	1.61	0.93	1.64	0.52	1.33	0.54	2.00
AE	0.62	1.44	0.71	1.42	0.15	1.36	0.48	2.34
BC	0.74	1.41	0.58	1.44	0.44	1.19	0.46	2.37
BD	0.93	1.64	0.95	1.63	0.26	1.48	0.67	2.00
BE	0.81	1.57	0.87	1.54	0.27	1.41	0.57	3.46
CD	1.12	1.74	1.10	1.75	0.42	1.48	0.64	2.98
CE	0.62	1.46	0.73	1.42	0.15	1.37	0.49	2.51
DE	0.43	1.21	0.30	1.24	0.32	1.07	0.28	3.43

Notice that the number of degrees of freedom η was set to be at least 2, and this limit was attained for AB, AD, and BD.

The choice of a specific copula class from the set H1–H6 is important. Let us consider the example of early failure of both components A (spark plug assembly) and B (ignition coil assembly) in the first year of the vehicle use. There was a total of K failures for the first 365 days out of N vehicles on the record, corresponding to the relative frequency of the first-year failure of . If we use the parameter values from Table 7.6, Clayton’s copula H1 would suggest the probability of failure p₁ = 0.365, copula H2 would yield p₂ = 0.371, BB1 copula H5 would yield p₅ = 0.373, and t-copula gives p₆ = 0.387.

7.6.2 Comparison of Copula Classes

We will first use information criteria to determine which of the six copula classes H1–H6 provides the best fit. Applying Akaike information criterion (AIC) and Bayes information criterion (BIC) brings about the results in Tables 7.7 and 7.8. The larger negative values correspond to the better models. The best values of the criteria are boldfaced. Two-parametric class H6 provides the best fit in most of the cases.

The Kolmogorov–Smirnov distance defined in (7.2) measures the maximum distance between a cumulative distribution function (c.d.f.) and its empirical cumulative distribution function (e.c.d.f.), and is applicable to joint distributions as well. Table 7.9 summarizes the Kolmogorov–Smirnov distances between the joint e.c.d.f.s and model c.d.f.s corresponding to the MPLE of parameters for each copula class.

From the computed Kolmogorov–Smirnov distance values, H5 and H6 are two best choices by far when compared to the other hypotheses. The distance values for H5 and H6 range from 0.0225 up to 0.06, whereas the other four combined range from 0.0653 up to 0.498. Each of H1–H4 misrepresents the e.c.d.f.s at some point, but H5 and H6 accurately represent the entirety of the data. H5 and H6 can be concluded as the obvious selections when comparing Kolmogorov–Smirnov distance values.

Tail dependence describes extreme co-movements between a pair of random variables in the tails of the distributions. The tail dependence coefficients for copulas C(u, v) was defined in (7.20) and (7.21). Certain copulas, such as the Gaussian copula, do not admit tail dependence (their tail dependence coefficients are equal to 0). Some copulas (e.g., H1–H4) exhibit only either upper or lower tail dependence, and the t-copula and BB1 family exhibit both. Thus, when selecting an accurate copula model for data it is important to consider whether the data displays upper and/or lower tail dependence. Formulas (7.22)–(7.26), from Section 7.2 can be used to calculate the lower and/or upper tail dependences of the selected five copula models: H1: λ_l = 2^{− 1/α}, H2: λ_u = 2 − 2^1/α, H3: λ_u = 2^{− 1/α}, H4: λ_l = 2 − 2^1/α, H5 reflects both lower and upper tail dependences, and they are allowed to be asymmetric: λ_l = 2^{− 1/αβ}, and λ_u = 2 − 2^1/β, H6: .

All 10 pairs in Table 7.10 exhibit both upper and lower tail dependence. Since H5 and H6 are the only two classes that allow for both, they demonstrate the best fit, though the symmetry of the tails is forced for H6. Lower values for H5, especially at the lower tails, may indicate some problems.

All of these approaches (using AIC, BIC, Kolmogorov–Smirnov statistic, or tail dependence) can help us to determine the best class of pair copula models. However they share the same weakness: the analysis is restricted to the comparison of the single representatives of each class obtained by CMLE. The following section discusses a possibility to make conclusions based on multiple representatives of six classes H1–H6.

7.6.3 Bayesian Model Selection

We suggest Kendall’s concordance τ as a universal parameter, which can be evaluated for all six classes H1–H6 defined earlier. Sample concordance is a reasonable nonparametric estimator of τ. Using parametric estimates from Table 7.5 and formulas expressing concordance through association parameters (see definitions of classes H1 and H2, (7.15) and (7.16)), we can calculate the values of τ induced by CMLE for parameters of copula classes H1–H6, and compare them to the sample values of , which is done in Table 7.11.

		H1	H2	H3	H4	H5	H6
AB	0.707	0.601	0.679	0.596	0.682	0.692	0.737
AC	0.269	0.224	0.264	0.222	0.269	0.278	0.278
AD	0.375	0.343	0.377	0.318	0.389	0.403	0.364
AE	0.314	0.235	0.307	0.262	0.291	0.316	0.320
BC	0.307	0.271	0.290	0.226	0.304	0.311	0.306
BD	0.400	0.318	0.389	0.323	0.387	0.402	0.471
BE	0.386	0.289	0.364	0.302	0.350	0.375	0.382
CD	0.452	0.360	0.425	0.354	0.429	0.442	0.443
CE	0.318	0.237	0.315	0.268	0.295	0.321	0.324
DE	0.188	0.177	0.174	0.130	0.193	0.194	0.182

Proximity of model-induced values of τ to the sample value may serve as a measure of the model fit and help to compare the model performance. However, this comparison is still using single values representing entire families.

Following Bretthorst [7] and Huard [20], in Section 7.4 we suggested to compare the data fit provided by models H1–H6 not at a single value of association parameter(s) obtained by MLE, but rather over the entire range of possible association values. This can be accomplished by specifying a prior distribution for association parameter(s) and integrating the likelihood with respect to the prior distribution. The problem is the difference of meaning and ranges of association parameters for different copula classes. However this problem is resolved by specifying a prior on τ.

There exist convenient formulas expressing association as a function of τ for all four one-parametric copulas H1–H4, inverting (7.15) and (7.16):

for H1 and H3, and

for H2 and H4.

In the case of the two-parametric class H5 we will use the conditional relationship between α and τ treating β as the nuisance parameter and using its MPLE estimate :

This approach can be justified by the higher risks of early related failures corresponding to the lower tails of the joint distribution represented by α, which becomes the critical parameter for the class H5. Notice that in the two-parametric class H6, there is an additional parameter η (degrees of freedom), which is not directly related to τ. We will apply (7.19) and use CMLE in the further analysis.

We will assume that six classes H1–H6 represent exhaustive and mutually exclusive hypotheses. Posterior probabilities of these hypotheses may be rewritten as

(7.39)

where we will consider all six hypotheses a priori equally likely, the dependence between variables being positive which suggests τ ≥ 0. In this case the natural choice of prior for τ is Beta distribution, and the choice of parameters for the prior is suggested by sample concordance for the entire dataset consistently with empirical Bayes approach: . Estimates of parameters of the Beta distribution are obtained from 10 pairs of components as and .

We will not need to calculate the denominator of the posterior in (7.39). It suffices to calculate the weights

(7.40)

or, using Monte Carlo approach and drawing samples from the Beta prior, evaluate

(7.41)

Then we choose the class with the highest weight and obtain the Bayes estimate of the association parameter using MCMC.

Table 7.12 gives the values of the natural logarithms of weights calculated by N = 20, 000 runs of Monte Carlo sampling from the Beta prior. The highest values for each pair of components are boldfaced.

7.6.4 Conclusions

As we can see from the warranty claim data, a substantial dependence is observed between the TTFs for automotive components related to engine subassembly of Hyundai Accent vehicles. This dependence cannot be ignored, because the assumption of independence would lead to grossly underestimated related failure risks and warranty costs. However it can be addressed using pair copula models. Open source software packages in R environment [6, 25] allow for a straightforward implementation of parametric and semiparametric estimation for different classes of copulas: direct and survival Archimedean copulas (Clayton and Gumbel–Hougaard classes), two-parameter Archimedean copulas (BB1), and elliptical copulas (Student t-class).

The problem of model selection a.k.a. an adequate choice of a copula class stays central, because misspecification of the class of copulas can be critical. For model comparison, one can consider using information criteria or Kolmogorov–Smirnov statistic as a good measure of overall fit, or studying tail dependence to assess the model quality. A convenient tool of model selection is provided by Bayesian framework using Kendall’s tau as the common parameter for different classes of copulas.

For TTFs of automotive components most of the tools of model selection indicate t-copulas being superior to four one-parameter Archimedean copulas and even to BB1. This fact can be explained by t-copulas exhibiting tail dependence in both lower and upper tails of the joint TTF distribution, while direct or survival Archimedean copulas concentrate at one of the tails or, as seems to be the case with BB1, underestimates tail dependence. However, it is not clear whether the assumption of symmetry of the tails of t-copulas is plausible and not too restrictive. One possible suggestion for future work is to use more complex hybrid or mixed Archimedean copula models as an alternative to t-copulas (see also Reference [26]). Another alternative would be to consider skewed or asymmetric t-copulas [9] or other more complex extensions of elliptical copula models.