Copula Simulation

1. Introduction

This lecture explains simulation by the method of copulas. We look first at generalities then some specific cases, in several categories:
1. Recap of Normal and T
2. Other copulas in 2D
3. Archimedian via transforms
For those of you who do not know about Laplace transforms, I shall give a cookbook approach for the 2D case so you can get going on the simulation. There are differing conventions here and there - I tend to follow Schmidt and Aas (see later).

2. Dependency Philosophy

In this framework the dependency is managed completely separately from the marginals. First we simulate the copula by (through some method we need to be careful about) constructing samples from the unit hypercube:

{(Subscript[U, 1]),(Subscript[U, 2]),…,(Subscript[U, n])}

where 0≤ (Subscript[U, i])≤ 1, with the given copula dependency structure. Then, given marginal quantile functions (Subscript[Q, i])(u) the sample of interest is constructed by

{(Subscript[Y, 1]),(Subscript[Y, 2]),…,(Subscript[Y, n])}={(Subscript[Q, 1])((Subscript[U, 1])),(Subscript[Q, 2])((Subscript[U, 2])),…,(Subscript[Q, n])((Subscript[U, n]))}

It is clear that the (Subscript[Q, i]) can be anything we choose, and that this relation makes clear the utility of quantile functions.

3. Copulas coming from a real multivariate density

So how do we make the list {(Subscript[U, 1]),(Subscript[U, 2]),…,(Subscript[U, n])}? If the copula is associated with a genuine multivariate probability density and furthermore,
(a) we can make samples {(Subscript[X, 1]),(Subscript[X, 2]),…,(Subscript[X, n])} from that multivariate density;
(b) the marginal CDFs for each univariate marginal are known as (Subscript[F, i])(x); 
then

{(Subscript[U, 1]),(Subscript[U, 2]),…,(Subscript[U, n])}={(Subscript[F, 1])((Subscript[X, 1])),(Subscript[F, 2])((Subscript[X, 2])),…,(Subscript[F, n])((Subscript[X, n]))}

4. The Gaussian copula

Let (Subscript[Z, i])be independent and standard N(0,1) random variables and consider

(Subscript[X, i])=(∑_(k=1)^n)(Subscript[A, i k]) (Subscript[Z, k])

for some matrix A. With the Einstein summation convention for repeated indices we can write this as:

(Subscript[X, i])=(Subscript[A, i k]) (Subscript[Z, k])

The covariance matrix of the (Subscript[X, i]) is then

(Subscript[C, i j])=Cov((Subscript[X, i]),(Subscript[X, j]))=Cov((Subscript[A, i k]) (Subscript[Z, k]),(Subscript[A, j m]) (Subscript[Z, m]))=(Subscript[A, i k]) (Subscript[A, j m]) Cov((Subscript[Z, k]),(Subscript[Z, m]))
=(Subscript[A, i k]) (Subscript[A, j m]) (Subscript[δ, k m])
= (Subscript[A, i k]) (Subscript[A, j k])
= (Subscript[(A.(A^T)), i j])

so that the desired covariance structure is recovered provided we can make a matrix A with the property that

C=A.(A^T)

There are two commonly used routes to this decomposition.

Cholesky Decomposition

The Cholesky decomposition of a symmetric positive definite matrix M is a factorization into a unique upper triangular G such that

M==(G^T).G

This factorization has a number of uses, one of which is that, because it is a triangular factorization, it can be used to solve systems of equations involving symmetric positive definite matrices.

The Cholesky factorization can be computed in Mathematica with the function CholeskyDecomposition.

If we take L = (G^T) then we obtain

M==L.(L^T)

where L  is lower triangular. Clearly this will do for the matrix A. There are good discussions and program implementations in all versions of Numerical Recipes, especially in Edition III for C++. There is also a decent discussion on wikipedia!
http://en.wikipedia.org/wiki/Cholesky_decomposition

2D Cholesky on the correlation

In two dimensions, if you first scale out the volatilities, it is easy to write down matters directly. Let

L = {{1,0},{ρ, Sqrt[1-ρ^2]}}

Then we have:

MatrixForm[L.Transpose[L]]

and the operation of L on standard normal variables preserve unit variance.

Diagonalization

This might seem a little over the top but I prefer to explicitly create the diagonalized system and check out the eigenvalues. This will tell you rather clearly if the covariance or correlation structure has been meddled with so that the matrices are no longer positive semi-definite. This is a good check. Let's pull out the key operations used in the last lecture to have another quick look at the key parts.

vols[n_]:=RandomReal[{0.1,0.5},n]

volatility = vols[5]

rawcorr[n_]:=Table[If[i==j,1,RandomReal[{-0.1,0.1}]],{i,1,n},{j,1,n}];
corr[n_]:=Module[{raw=rawcorr[n]},(1/2) (raw+Transpose[raw])]

Note - in general you absolutely cannot invent correlation matrices like this. We are taking small correlations so will stay positive definite!

Let's make a 5 by 5 version:

AssetCorr = corr[5];
TableForm[N[AssetCorr, 1]]

AssetCov =  Outer[Times, volatility, volatility]*AssetCorr

{LocalEvals, LocalEvecs} = Eigensystem[AssetCov];

A matrix that is zero apart from diagonal entries that are the pseudo-vols of the local independent processes is:

LocalDiag = Table[If[i==j, Sqrt[LocalEvals[[i]]], 0], {i, 5}, {j, 5}]

Acklam= Compile[{{u, _Real}},
Module[
{a = Reverse[{-39.69683028665376, 220.9460984245205, -275.9285104469687, 138.3577518672690, -30.66479806614716, 2.506628277459239}],
b = Reverse[{-54.47609879822406, 161.5858368580409,-155.6989798598866, 66.80131188771972,-13.28068155288572}],c=Reverse[{-0.007784894002430293,-0.3223964580411365,-2.400758277161838,-2.549732539343734,4.374664141464968,2.938163982698783}],
d=Reverse[{0.007784695709041462,0.3224671290700398,2.445134137142996,3.754408661907416}]},
Which[0.02435≤ u≤ 0.97575,
Module[{v = u-1/2, r = (u-1/2)^2},
v*(a[[1]]+r*(a[[2]]+r*(a[[3]]+r*(a[[4]]+r*(a[[5]]+r*a[[6]])))))/(1+r*(b[[1]]+r*(b[[2]]+r*(b[[3]]+r*(b[[4]]+r*b[[5]])))))],
u> 0.97575, 
Module[{q=Sqrt[-2*Log[1-u]]},-(c[[1]]+q*(c[[2]]+q*(c[[3]]+q*(c[[4]]+q*(c[[5]]+q*c[[6]])))))/(1+q*(d[[1]]+q*(d[[2]]+q*(d[[3]]+q*d[[4]]))))],
True,
Module[{q=Sqrt[-2*Log[u]]},(c[[1]]+q*(c[[2]]+q*(c[[3]]+q*(c[[4]]+q*(c[[5]]+q*c[[6]])))))/(1+q*(d[[1]]+q*(d[[2]]+q*(d[[3]]+q*d[[4]]))))]]]];

Sets of size n can then be made rather trivially (this is a nice thing about quantile maps!)

setgen[n_]:=Table[Acklam[RandomReal[]],{n}];

vals = Table[setgen[5], {1000}];

CorrSet = Map[Transpose[LocalEvecs].LocalDiag.#&, vals];

(* Needs["MultivariateStatistics`"];Needs["ComputationalGeometry`"]*)

simcov=Covariance[CorrSet]

So now let's look at the ratio of the simulated to the actual covariance:

simcov/AssetCov

The maximum absolute difference in the entries is:

Max[Flatten[Abs[simcov-AssetCov]]]

Checking for trouble

Let's be a bit more cavalier with some invented correlations:

rawcorr[n_]:=Table[If[i==j,1,RandomReal[{-1,1}]],{i,1,n},{j,1,n}];
corr[n_]:=Module[{raw=rawcorr[n]},(1/2) (raw+Transpose[raw])]

AssetCorr = corr[5];
TableForm[N[AssetCorr, 1]]

Looks plausible enough....

AssetCov =  Outer[Times, volatility, volatility]*AssetCorr

Ditto, but.....

CholeskyDecomposition[AssetCov]

You get more information this way:

Eigenvalues[AssetCov]

This cannot have come from real data.

From the distribution to the copula

Having made our samples we apply the CDF to get to the copula. Let's drop down to dimension two to see what is happening. In this case we can use the easy representation of the Cholesky formula to build the mixing matrices.

ρ = 0.5;
L = {{1,0},{ρ, Sqrt[1-ρ^2]}}

uncorrsamp = Table[setgen[2], {i, 1, 1000}];

corrsamp = Map[L.#&, uncorrsamp];

ListPlot[corrsamp]

Now we go get the CDF (recall we know how to do this in C++ with either a rational approximation or the Marsaglia formula) here we will be lazy and do it with the built in gadget:

Ncdf[x_]:=(1/2) (Erf[(x/(Sqrt[2]))]+1)

copsamp = Map[Ncdf[#]&, corrsamp];

ListPlot[copsamp,AspectRatio->1]

Other correlations

Keeping the same underlying random data, we can fiddle with the correlations and see the impact on the copula samples.

ρ=0.9`;
L={{1,0},{ρ,(Sqrt[1-(ρ^2)])}};
corrsamp=(L.#1&)/@uncorrsamp;
copsamp=(Ncdf[#1]&)/@corrsamp;
ListPlot[copsamp,AspectRatio->1]

ρ=-0.95`;
L={{1,0},{ρ,(Sqrt[1-(ρ^2)])}};
corrsamp=(L.#1&)/@uncorrsamp;
copsamp=(Ncdf[#1]&)/@corrsamp;
ListPlot[copsamp,AspectRatio->1]

ρ=0;
L={{1,0},{ρ,(Sqrt[1-(ρ^2)])}};
corrsamp=(L.#1&)/@uncorrsamp;
copsamp=(Ncdf[#1]&)/@corrsamp;
ListPlot[copsamp,AspectRatio->1]

5. The T variation(s)

The machinery developed above can all be re-jigged to make samples from a T distribution and its associated copula. In the good old-fashioned multivariate T, of often known as "the T copula", the only thing that changes is that having mixed up the normal samples, they are all divided by the square root of the same normalized chi-square variable. As shown in my paper with Kim Lee in the Journal of Multivariate Analysis, this may be generalized in all sorts of ways.

Needs["Statistics`ContinuousDistributions`"]

Needs["Statistics`NormalDistribution`"]

Random[ChiSquareDistribution[4]]

setgen[n_]:=((Table[Acklam[RandomReal[]],{n}])/(Sqrt[(1/2) Random[ChiSquareDistribution[2]]]));

uncorrsamp = Table[setgen[2], {i, 1, 10000}];

The T CDF

You can find examples of formulae for these things in my paper in JCF. It is paticularly easy when the degrees of freedom is an even integer, when there is a neat recurrence formula.

Clear[a]

a[0, n_] := Gamma[(n+1)/2]/Gamma[n/2]/Sqrt[n Pi];
a[k_, n_] := a[k,n] = (n-2 k)/n/(2k+1) a[k-1, n];

TCDF[n_, x_] := 
1/2 + x*Sum[a[p, n] x^(2p), {p, 0, n/2-1}]/(1+x^2/n)^(1/2(n-1))

TCDF[2,x]

ρ=0;
L={{1,0},{ρ,(Sqrt[1-(ρ^2)])}};
corrsamp=(L.#1&)/@uncorrsamp;
copsamp=(TCDF[2,#1]&)/@corrsamp;
ListPlot[copsamp,AspectRatio->1, PlotStyle -> PointSize[0.005]]

This plot illustrates the phenomenon of asymptotic tail dependence, which may be quantified.

6. Other copulas in two dimensions

It is useful to make corresponding analysis of some other copulas in dimension two, as then the simulation remains straightforward. Let's look at some examples. In this case we do not work from some other distribution, but go direct for samples from the unit square, i.e., directly from the copula. I will focus on simulation methods rather than dwell on the motivation. This, the choice, calibration etc are all headaches! Having thought about how what order to do things, I will fist supply you with some "cookbook recipes" for the 2D case, then delve a little deeper into where things come from and how to go to higher dimensions. Then at least you will be able to do some simulation right away. Note that the focus here will be on modern methods based more on the Laplace transform methodology. There are other options, particularly in 2D, that are based on the "conditional sampling" methodology. I will indeed quote the result for the Frank case for 2D. However, this gets messy and computationally intensive when you go to dimensions greater than two. For details of this you are referred to Section 6.3 of "Copula Methods in Finance", Cherubini et al.

7. Marshall-Olkin Copulas

This has been discussed by Schmidt in his nice summary paper (Algorithm 5.4), though there is at least one typo, and an apparent disagreement over conventions with Embrechts, Lindskog and McNeil. 
1. Generate independent random variables (Subscript[Z, i]) with expontial distributions with paramaters (Subscript[λ, i]), via the standard exponential (reversed-)quantile based simulation, with (Subscript[V, i])uniform and independent:

(Subscript[Z, i])=-((log((Subscript[V, i])))/(Subscript[λ, i]))

2. Return, for j=1,2,

(Subscript[U, j])=(e^(-((Subscript[λ, j])+(Subscript[λ, 3])) Min((Subscript[Z, j]),(Subscript[Z, 3]))))

MarshallOlkin[lamone_,lamtwo_,lamthree_,n_]:=Module[{raw=Table[{-((Log[RandomReal[]])/lamone),-((Log[RandomReal[]])/lamtwo),-((Log[RandomReal[]])/lamthree)},{n}],maxed},maxed=({Min[#1[[1]],#1[[3]]],Min[#1[[2]],#1[[3]]]}&)/@raw;({(E^(-(lamone+lamthree) #1[[1]])),(E^(-(lamtwo+lamthree) #1[[2]]))}&)/@maxed]

MOdata = MarshallOlkin[1.1,0.2, 0.6, 2000];

ListPlot[MOdata,AspectRatio->1,PlotRange->All]

Compare this with figure 4.1 of Embrechts, Lindskog and McNeil, 2001 - Modelling Dependence with Copulas and Applications to Risk management. It seems to agree with regard to orientation. My conventions also agree with Kelly (Nov 2007, "Using copulas to model dependence in simulation risk management"). These copulas may also be parametrized by

α=((Subscript[λ, 3])/((Subscript[λ, 1])+(Subscript[λ, 3])));
β=((Subscript[λ, 3])/((Subscript[λ, 2])+(Subscript[λ, 3])));

There is a natural generalization to many dimensions (See ELN, 2001)

8. Some "Archimedean" cases in 2D

We shall define what this means more carefully presently, but for now let's look at some simple examples of 2D simulation.

Frank 2-copula

In this case the simulation of a pair of correlated uniform deviates reduces to the following algorithm. Let ((Subscript[v, 1]), (Subscript[v, 2])) be independent samples from a uniform distribution on [0,1]. Then we set (Subscript[U, 1])=(Subscript[v, 1])and then, as a function of the copula parameter α,

(Subscript[u, 2]) = -1/α Log[1 + (Subscript[v, 2])(1 - Exp[-α])/((Subscript[v, 2])*(Exp[-α (Subscript[u, 1])]-1)- Exp[-α (Subscript[u, 1])])]

(Subscript[u, 2])=-(1/α)log(1+(((1-(e^(-α))) (Subscript[v, 2]))/((Subscript[v, 2]) ((e^(-α (Subscript[u, 1])))-1)-(e^(-α (Subscript[u, 1]))))))

Frankcorrpair[alpha_]:=Module[{uone=RandomReal[],utwo,vtwo=RandomReal[],emau},emau=(E^(-alpha uone));utwo=If[alpha==0,vtwo,-((Log[1+((vtwo (1-(E^(-alpha))))/(vtwo (emau-1)-emau))])/alpha)];{uone,utwo}]

SeedRandom[100];
frankten=Table[Frankcorrpair[10],{10000}];
franktenplot=ListPlot[frankten,AspectRatio->1]

SeedRandom[100];
frankten=Table[Frankcorrpair[-10],{10000}];
ListPlot[frankten,AspectRatio->1]

SeedRandom[100];
frankten=Table[Frankcorrpair[0],{10000}];
ListPlot[frankten,AspectRatio->1]

SeedRandom[100];
frankten=Table[Frankcorrpair[3],{10000}];
ListPlot[frankten,AspectRatio->1]

DensityFunction for Frank

The copula itself is given by:

FrankCop[u_, v_, α_] := Module[{expal = Exp[-α]},
-1/α Log[1 + (expal^u-1)*(expal^v-1)/(expal-1)]]

FrankCop[u,v,a]

To get the copula density for comparison with the scatterplots we differentiate the copula with respect to both u and v:

FrankDenRaw[u_, v_, α_] =Simplify[D[FrankCop[u,v,α], u, v]]

FrankDen[u_,v_,α_]:= FrankDenRaw[u, v, α]

FrankDen[0.5, 0.5, 10]

frankdenten=ContourPlot[FrankDen[u,v,10],{u,0.001`,1},{v,0.001`,1}]

Show[GraphicsRow[{franktenplot,frankdenten}]]

Clayton 2-copula

This is easily done if one has a gamma simulator, as discussed by K. Aas (2004) "Modelling the dependence structure of financial assets: A survey of four copulas." The Clayton case is good when there is one-sided tail dependence.

?GammaDistribution

Random[GammaDistribution[2,1]]

The procedure for this case is as follows, where the copula has a parameter δ:

1. Simulate a Gamma variate X distributed as Gamma(1/δ, 1).
2. Simulate 2 (more generally d) independent standard uniforms (Subscript[V, 1]), (Subscript[V, 2]).
3. Return

(Subscript[U, 1])=(((1-((log((Subscript[V, 1])))/X)))^(-1/δ));
(Subscript[U, 2])=(((1-((log((Subscript[V, 2])))/X)))^(-1/δ))

Claytoncorrpair[δ_]:=Module[{vone=RandomReal[],vtwo=RandomReal[],X=Random[GammaDistribution[(1/δ),1]]},{(((1-((Log[vone])/X)))^(-1/δ)),(((1-((Log[vtwo])/X)))^(-1/δ))}]

Claytoncorrpair[1]

SeedRandom[100];clayten=Table[Claytoncorrpair[0.00001`],{10000}];
clayplotsmall=ListPlot[clayten,AspectRatio->1]

SeedRandom[100];clayten=Table[Claytoncorrpair[1],{10000}];
clayplotone=ListPlot[clayten,AspectRatio->1]

SeedRandom[100];clayten=Table[Claytoncorrpair[3],{10000}];
clayplotthree=ListPlot[clayten,AspectRatio->1]

SeedRandom[100];clayten=Table[Claytoncorrpair[10],{10000}];
clayplotten=ListPlot[clayten,AspectRatio->1]

The Clayton density function

This time I will just write this down directly:

ClayCopDen[u_, v_, δ_] := (1+δ)(u*v)^(-1-δ)(u^(-δ)+v^(-δ))^(-1/δ-2)

ClayCopDen[u, v, δ]

(((u v))^(-δ-1)) ((((u^(-δ))+(v^(-δ))))^(-2-(1/δ))) (δ+1)

ClayCopDen[0.5, 0.5, 3]

claydensmall=ContourPlot[ClayCopDen[u,v,0.0001`],{u,0.001`,1},{v,0.001`,1}]

claydenone=ContourPlot[ClayCopDen[u,v,1],{u,0.001`,1},{v,0.001`,1}]

claydenthree=ContourPlot[ClayCopDen[u,v,3],{u,0.001`,1},{v,0.001`,1}]

claydenten=ContourPlot[ClayCopDen[u,v,10],{u,0.001`,1},{v,0.001`,1}]

Show[GraphicsGrid[{{clayplotsmall,claydensmall,clayplotone,claydenone},{clayplotthree,claydenthree,clayplotten,claydenten}}]]

See Section 6.3.1 of "Copula Methods in Finance" for a discussion of the conditional approach - I will come back to d>3 for Clayton later.

In general, before doing any detailed applications, it is not a bad idea to cross-check
(a) the copula (joint CDF on the square)
(b) the density function on the square
(c) the form of your scatterplots from the proposed simulation method.

9. Digging Deeper into Archimedian Copulas

The recipe approach presented so far gets us simulating quickly in two dimensions but leaves us a little short on insight and generalizations to higher dimensions. We can make more sense of matters by proceeding more formally and introducing a Laplace transform description. The d-dimensional Archimedian copulas may be written as

C((Subscript[u, 1]),…,(Subscript[u, d]))=(ϕ^(-1))(ϕ((Subscript[u, 1]))+…+ϕ((Subscript[u, d])))

The function ϕ is a decreasing function called the generator of the copula. (ϕ^(-1)) denotes the inverse of the generator. The generators for some common cases are:

Frank Generator

(Subscript[ϕ, F])(u,α)=log((((e^(-α))-1)/((e^(-α u))-1)))

Clayton Generator

(Subscript[ϕ, Cl])(u,θ)=(u^(-θ))-1

Gumbel generator

We have not yet discussed this one but it has the simple generator

(Subscript[ϕ, G])(u,θ)=(((-log(u)))^θ)

In general these things make sense provided ϕ is a strictly decreasing convex function mapping [0,1] to [0, ∞] with ϕ(1)=0. See the book by Nelsen (1999) for details.

Laplace transform descriptions

The idea is to regard (ϕ^(-1))(t) as the Laplace transform of an associated density function g(x) with associated CDF G(x). That is, we set

InversePhi[t] = Integrate[Exp[-t x] g[x], {x, 0, Infinity}]

(ϕ^(-1))(t) = (∫_0^∞)(e^(-t x)) g(x)dx

Let's look at these associated densities first, before saying what we do with them.

Density in the Clayton case

In this case the generator and its inverse are just

(Subscript[ϕ, Cl])(u,θ)=(u^(-θ))-1

((Subscript[ϕ, Cl])^(-1))(t,θ)=(((1+t))^(-1/θ))

InverseLaplaceTransform[(1+t)^(-1/θ), t, x]

To recognize this we pull out the gamma density function

Needs["Statistics`ContinuousDistributions`"];
PDF[GammaDistribution[a,1], x]

and we see that the invented density is that of a gamma distribution with a=1/θ.

Density in the Frank case

This is actually discrete! It is easiest to quote the density and check its transform:

Sum[(1-Exp[-α])^k/(k α) Exp[-k t], {k, 1, Infinity}]

and this is the functional inverse of

t = (Subscript[ϕ, F])(u,α)=log((((e^(-α))-1)/((e^(-α u))-1)))

So the associated density is actually discrete with P(V=k) given by, for k=1,2,…,

Density in the Gumbel case

The inverse generator is given by

Exp[-t^(1/θ)]

Mathematica chokes on this one:

InverseLaplaceTransform[Exp[-t^(1/θ)], t, x]

It only comes out nicely for a few values of θ, e.g. 2:

InverseLaplaceTransform[Exp[-t^(1/2)], t, x]

Integrate[%, {x, 0, y}, Assumptions -> y>0]

Plot[%,{y,0,10}]

For general values of θ we obtain one of the so-called stable distributions, which are in fact most easily characterized by their transforms.

The algorithm for d-dimensional Archimedian

This may now be stated simply - for the proof see McNeil, Frey Embrechts (2005). 
1. Generate a random variable V with density g, CDF G;
2. Generate uniform IID (Subscript[X, 1]),(Subscript[X, 2]),(Subscript[…X, d]);
3. Let

(Subscript[U, i])=(ϕ^(-1))[-((log((Subscript[X, i])))/V)]

In practice making V is the hard part unless we are in the Clayton case, where can use Gamma. It is interesting to consider numerical transform methods as well! I do not know, and would be very interested in, any methodology that employs the quantile methodology for this problem. The Frank formula for 2D is a special case that avoids the problem. An algorithm for stable distributions is gievn by K. Aas. 

For non Archimedian in general dimensions one is faced with a difficult iterative process using conditional densities (see Cherubini et al, Copulas in finance).

Clayton in d=3

Claytoncorrtrip[δ_]:=Module[{vone=RandomReal[],vtwo=RandomReal[],vthree=RandomReal[],X=Random[GammaDistribution[(1/δ),1]]},{(((1-((Log[vone])/X)))^(-1/δ)),(((1-((Log[vtwo])/X)))^(-1/δ)),(((1-((Log[vthree])/X)))^(-1/δ))}]

Claytoncorrtrip[1]

simdata = Table[Claytoncorrtrip[5], {i, 1, 1000}];

Show[Graphics3D[Point/@simdata]]

Show[%,ViewPoint->{0,0,100}]

Show[%,ViewPoint->{100,0,0}]

Kendall's tau for Archimedian copulae

There are some nice properties of such copulae that facilitate calibration. In general, to encode dependency, we use concepts like correlation, Spearman's rank correlation and Kendall's τ. (Look them up on Wiki if you do not know the definitions in terms of data - there are nice discussions in terms of data). The Spearman and Kendall quantities are actually functions of the copula. For example, Kendall's τ may be estimated from data. In terms of a general 2D copula it may be written (this is not obvious!)

τ=4 E[C[U,V]]-1

In the 2D Archimedian case this can be shown to be given by

τ=4 (∫_0^1)((ϕ(t))/(ϕ'(t)))dt+1

For Clayton this gives

p[u_]:=u^(-θ)-1;Simplify[4*Integrate[p[u]/p'[u], {u, 0, 1}, Assumptions -> θ>-1]+1]

For Gumbel this gives

p[u_]:=(-Log[u])^θ;
Simplify[4*Integrate[p[u]/p'[u], {u, 0, 1}, Assumptions -> θ>-1]+1]

For Frank this gives

p[u_]:= Log[(E^(-α) - 1)/
    (E^((-α)*u) - 1)];
    Simplify[4*
      Integrate[p[u]/p'[u], {u, 0, 1}, Assumptions -> θ>-1]+1]

Plot[%,{α,0.1,30}]

So this does make sense - the function Li in there is a polylogarithm. Embrechts et al write the answer in terms of Debye functions....

There are many other aspects of calibration, and choosing a copula. Doing some plots is always a good first step, and one may often wish to use measures of tail dependence as well.

References

K. Aas, 2004, Modelling the dependence structure of financial assets: A survey of four copulas. Norwegian Computing Centre, report SAMBA/22/04, December 2004. 
http://www.nr.no/files/samba/bff/SAMBA2204c.pdf

U. Cherubini, E. Luciano, W. Vecchiato, 2004, Copula Methods in Finance, Wiley.

P. Embrechts, F. Lindskog, A.J. McNeil, 2003, Modelling dependence with copulas and applications to risk management. In Handbook of heavy tailed distributions in finance, edited by Rachev ST, published by Elsevier/North-Holland, Amsterdam.  A 2001 preprint version is available here:
http://www.risklab.ch/ftp/papers/DependenceWithCopulas.pdf

D.L. Kelly, 2007, Using Copulas to Model Dependence in Simulation Risk Management,
ASME International Mechanical Engineering Congress and Exposition, INL/CON-07-13037 preprint. 
http://www.inl.gov/technicalpublications/Documents/3752703.pdf

T. Schmidt, 2006, Coping with Copulas, In "Copulas - from Theory to Applications in Finance", RISK Books. Also at:
http://www.math.uni-leipzig.de/~tschmidt/TSchmidt_Copulas.pdf

W.T. Shaw and K.T.A. Lee, 2007, Bivariate Student t distributions with variable marginal degrees of freedom and independence, 
doi:10.1016/j.jmva.2007.08.006