Exact Simulation of reciprocal Archimedean copulas
Jan-Frederik Mai






plain
	

		Exact simulation of reciprocal Archimedean copulas
		

		
		Jan-Frederik Mai

		

		XAIA Investment

		Sonnenstr. 19, 80331 München

		email: jan-frederik.mai@xaia.com,

		phone: +49 89 589275-131.

	


The decreasing enumeration of the points of a Poisson random measure whose mean measure is Radon on (0,∞] can be represented as a non-increasing function of the jump times of a standard Poisson process. This observation allows to generalize the essential idea from a well-known exact simulation algorithm for arbitrary extreme-value copulas to copulas of a more general family of max-infinitely divisible distributions, with reciprocal Archimedean copulas being a particular example. 



§ INTRODUCTION

A copula C:[0,1]^d → [0,1] is a multivariate distribution function of a random vector whose components are all uniformly distributed on [0,1], cf. <cit.> for background.
The family of reciprocal Archimedean copulas has been introduced and analyzed in <cit.>. With a parameterizing univariate distribution function F:[0,∞) → [0,1], a copula in this class has the analytical form

    C_F(u_1,…,u_d) = ∏_A ∈𝒫_d,oF( ∑_k ∈ AF^-1(u_k))/ ∏_A ∈𝒫_d,eF( ∑_k ∈ AF^-1(u_k)),

where 𝒫_d,o (resp. 𝒫_d,e) is the set of all non-empty subsets of {1,…,d} with odd (resp. even) cardinality. The nomenclature of this copula family is justified by some striking analogies with the well-understood family of Archimedean copulas, cf. <cit.> for background on the latter. For instance, <cit.> show that C_F is a proper copula in dimension d if and only if the parameterizing function F:[0,∞) → [0,∞) has the form 

    F(t) = exp( -∫_t^∞( 1-t/x)^d-1 ν(dx))=:exp(-Λ(t)),

with ν a non-finite Radon measure on (0,∞] such that ν({∞})=0, called the radial measure. In this case, a random vector 𝐘=(Y_1,…,Y_d) with distribution function C_F(F(t_1),…,F(t_d)) has stochastic representation

    𝐘 = (max_k≥ 1{R_k Q_1^(k)},…,max_k≥ 1{R_k Q_d^(k)}),

where 𝐐^(k)=(Q_1^(k),…,Q_d^(k)), k ∈, are independent and uniformly distributed on the d-dimensional simplex S_d:={𝐪=(q_1,…,q_d)≥ 0 : q_1+…+q_d=1}, and, independently, {R_k}_k ≥ 1 is an enumeration of the points of a Poisson random measure on (0,∞] with mean measure ν. It is important to notice that the index k in this maximum runs through an enumeration of the points {R_k}_k≥ 1 of the Poisson random measure. This collection of points is almost surely countably infinite since the measure ν is non-finite[As an educational side remark, <cit.> add in each component of (<ref>) a zero in the set of which the maximum is taken over. In the context of general max-infinitely divisible distributions, cf. Remark <ref> below, this is usual and necessary, since for finite mean measure one might otherwise encounter an empty set. However, in the present situation of non-finite ν this is unnecessary.]. The most prominent member of the family of reciprocal Archimedean copulas is the Galambos copula, the terminology dating back to <cit.>, which arises for the choice ν(dx)=Γ(d+1/θ)/(Γ(d) Γ(1/θ)) x^-1/θ-1 dx, cf. <cit.>. Further background on the Galambos copula can be found in <cit.>. 
The stochastic representation (<ref>) is difficult to simulate from due to the infinite maximum, which is why <cit.> only propose an approximative simulation strategy. For extreme-value copulas, a family whose intersection with reciprocal Archimedean copulas equals the Galambos copula, an alternative and exact simulation strategy is developed in <cit.>, based on an idea originally due to <cit.>. Section <ref> shows how the idea of this algorithm can be generalized to include arbitrary reciprocal Archimedean copulas. Section <ref> concludes.



§ EXACT SIMULATION OF RECIPROCAL ARCHIMEDEAN COPULAS

Let ν be a Radon measure on (0,∞], i.e. S(t):=ν((t,∞])<∞ for all t>0, with the property that ν({∞})=0. We denote u_ν:=ν( (0,∞]) and define a pseudo-inverse via[Notice that both S and S^-1 are non-increasing and right-continuous by definition.] 

    S^-1(t):=inf{x>0 : S(x) ≤ t},   t ∈ [0,u_ν ].

For later reference we remark that right-continuity of S implies

    y < S(x)  ⇔ S^-1(y) > x,   x ∈ (0,∞), y ∈ [0,u_ν ].

We explicitly allow ν to be finite here, in which case S^-1(t) is only defined for t ≤ u_ν<∞, but for the application to reciprocal Archimedean copulas only the case when u_ν=∞ is relevant. For x>0 we denote by δ_x the Dirac measure at x. We denote by P=∑_k = 1^Nδ_R_k a Poisson random measure on (0,∞] with mean measure ν, the random variable N ∈_0 ∪{∞} representing the number of points {R_k} of P. <cit.> is an excellent textbook for background on Poisson random measures. The most important, and characterizing, property of a Poisson random measure on a measurable space E with mean measure ν is the Laplace functional formula

    [ e^-∫_Ef(x) P(dx)] = exp(-∫_E ( 1-e^-f(x)) ν(dx) ),

where ν is a Radon measure on E, and f is a non-negative, Borel-measurable function on E. Recall that if ν is non-finite, P has countably many points R_1,R_2,…. But if ν is finite, the number of points N has a Poisson distribution with parameter u_ν. Hence, regarding notation it is convenient for us to treat both cases jointly by denoting the points of P by R_1,…,R_N, possibly allowing for the value N=∞ in case of non-finite ν. Without loss of generality, we further enumerate the points R_k such that R_1≥ R_2≥… almost surely.
The following auxiliary result follows from (<ref>) by a change of variables from (the possibly complicated measure) ν(dx) to the Lebesgue measure dx, resulting in a stochastic representation of {R_k} that is convenient for our purpose of simulating reciprocal Archimedean copulas. Even though this computation is presumably standard in the literature on Poisson random measures, we state it as a separate lemma and provide a proof here, because it is educational, one key ingredient for the derived simulation algorithm, and apparently lesser known in the literature on copulas and dependence modeling.


Let {ϵ_k}_k ∈ be a sequence of independent and identically distributed exponential random variables with unit mean. Introducing the random variable 

    N_ϵ:=∞   u_ν = ∞,
    ∑_k ≥ 11_{ϵ_1+…+ϵ_k ≤ u_ν}   
    ,

we have the distributional equality

    {R_k}_k = 1,…,Nd={S^-1(ϵ_1+…+ϵ_k)}_k = 1,…,N_ϵ.



Define the point measure P̂:=∑_k =1^N_ϵδ_S^-1(ϵ_1+…+ϵ_k). We notice that P̃=∑_k ≥ 1δ_ϵ_1+…+ϵ_k equals a Poisson random measure on [0,∞) with mean measure the Lebesgue measure dx, hence

    [ e^-∫_(0,∞]f(x) P̂(dx)] = [ e^-∑_k = 1^N_ϵf( S^-1(ϵ_1+…+ϵ_k))]  
          = [ e^-∫_[0,∞)f(S^-1(x)) 1_{x ≤ u_ν} P̃(dx)]
         (<ref>)=exp( -∫_[0,∞)( 1-e^-f(S^-1( x)) 1_{x ≤ u_ν}) dx)  
          =exp( -∫_0^u_ν( 1-e^-f(S^-1( x))) dx) (∗)=exp( -∫_0^∞( 1-e^-f(x)) ν(dx))
         (<ref>)=[ e^-∫_(0,∞]f(x) P(dx)].

To verify equation (∗), denote by λ the Lebesgue measure on (0,∞), and observe that the map G:=S^-1:(0,u_ν)→ (0,∞) is measurable. Consider the measure G_λ defined by G_λ(E):=λ(G^-1(E)), E a Borel set in (0,∞) and G^-1(E) its pre-image under G in (0,u_ν). Then we observe for x ∈ (0,∞) that

    G_λ( (x,∞])    = λ( {y ∈ (0,u_ν) :  S^-1(y) > x})(<ref>)=λ( {y ∈ (0,u_ν) : y < S(x) })
        = S(x) = ν((x,∞]).

Consequently, G_λ=ν and we have the measure-theoretic change of variable formula

    ∫_(0,∞)g(x) ν(dx)   =∫_(0,∞)g(x) G_λ(dx) = ∫_(0,u_ν)g(G(x) ) λ(dx)
        =∫_0^u_νg(S^-1(x) ) dx.

Applying it to the function g(x)=1-e^-f(x) implies (∗). The claim now follows from uniqueness of the Laplace functional of Poisson random measure, since f was an arbitrary non-negative, Borel-measurable function.



It is not the first time that the change of variables technique of Lemma <ref> is found useful for an application to simulation. To provide another example, <cit.> uses essentially the same technique to represent a non-negative[<cit.> in his article also considers infinitely divisible random variables on .] infinitely divisible random variable X with associated Lévy measure[In comparison with the measure ν of the present article, a Lévy measure satisfies the additional integrability condition ∫_0^1x ν(dx)<1.] ν as

    X = ∑_k =1^NR_k d=∑_k=1^N_ϵS^-1(ϵ_1+…+ϵ_k)
 
and discusses the possibility to simulate X based on this stochastic representation.


If the radial measure ν is absolutely continuous with positive density on (0,∞), the function S is continuous and strictly decreasing and S^-1 is the regular inverse. The following example sheds some light on the situation in the case of discrete measures ν.

[Discrete radial measures]
Let ∞>a_1>a_2>…>0 with lim_k →∞a_k=0 and b_k ≥ 0 with ∑_k ≥ 1b_k = ∞ (to guarantee that ν is non-finite), and consider ν = ∑_k ≥ 1b_k δ_a_k on (0,∞), obviously Radon on (0,∞]. The functions S and S^-1 in this case are given by

    S(t)    = ∑_k≥ 1( ∑_i=1^kb_i) 1_[a_k+1,a_k)(t),   t ≥ 0,
    
    S^-1(t)   =∑_k≥ 1a_k 1_[∑_i=1^k-1b_i,∑_i=1^kb_i )(t),   t ≥ 0.

As a concrete one-parametric example, let θ>0 and a_k=1/k, b_k=θ, k ≥ 1. It is not difficult to observe that in this case the formulas above boil down to

    S(t) = lim_u ↓ tθ ⌊1/u⌋,   S^-1(t)=1/⌈t/θ⌉,

with ⌊ . ⌋ and ⌈ . ⌉ denoting the usual floor and ceiling functions mapping [0,∞) to _0. The associated distribution function F=exp(-Λ) generating the d-dimensional reciprocal Archimedean copula associated with ν is determined by

    Λ(t) = θ ∑_k =1^⌊ 1/t ⌋( 1-k t)^d-1(d=2)=θ ⌊1/t⌋ ( 1-t {⌊1/t⌋+1}/2),

where the last equation in the bivariate case is stated explicitly for later reference.


Next, we turn to the simulation of reciprocal Archimedean copulas and introduce the notation

    𝐘_n   :=(max_k = 1,…,n{R_k Q^(k)_1},…,max_k = 1,…,n{R_k Q^(k)_d}),
    
     M_n   :=𝐘_n,   n ≥ 1.

It is observed that every single component of R_n+1 𝐐^(n+1) is smaller or equal than R_n, since the sequence {R_n}_n ≥ 1 is non-increasing by our enumeration. This implies 

    𝐘 = 𝐘_∞ = 𝐘_M, M:=min{n ≥ 1 : R_n+1≤ M_n}.
 
Since R_n is almost surely decreasing to zero and M_n is almost surely non-decreasing, M is almost surely finite. Consequently, in order to simulate 𝐘, it is sufficient to simulate iteratively 𝐘_n for n=1,2,…, until the stopping criterion R_n+1≤ M_n takes place, i.e. until n=M. This is precisely the simulation idea of Algorithm 1 in <cit.> for extreme-value copulas, see also <cit.>, enhanced to fit the scope of reciprocal Archimedean copulas as well with the help of Lemma <ref>. Algorithm <ref> summarizes this strategy in pseudo code. It requires evaluation of S^-1 and of F, which are the sole numerical obstacles. Notice further that we propose to simulate the random vectors 𝐐^(k) according to the well-known stochastic representation

    𝐐^(k)d=( E_1/E_1+…+E_d,…,E_d/E_1+…+E_d),

where E_1,…,E_d are independent exponential random variables with unit mean, cf. <cit.>.




If the random vectors 𝐐^(k) in (<ref>) are not uniformly distributed on S_d, but instead follow some other distribution on S_d, Algorithm <ref> can still be used for simulation, provided one has at hand a simulation algorithm for 𝐐^(k). In this case, one breaks out of the cosmos of reciprocal Archimedean copulas. In the particular case ν(dx)=d/x^2 dx this generalized algorithm equals precisely <cit.> for arbitrary extreme-value copulas. However, for other radial measures ν one also breaks out of the cosmos of extreme-value copulas. The resulting distribution function of 𝐘 in the general case is

    (𝐘≤𝐲) = exp(-[S( min{y_1/Q_1^(1),… ,y_d/Q_d^(1)}) ]),

for 𝐲=(y_1,…,y_d) ∈ [0,∞)^d. For x>0 the function t ↦ x S(t)=x ν((t,∞]) is the survival function of the Radon measure x ν, showing that the d-variate function (𝐘≤𝐲)^x is again a distribution function. Multivariate distribution functions with this property are called max-infinitely divisible, see <cit.> for background. 


Concerning the implementation of Algorithm <ref>, the biggest numerical difficulty is the evaluation of the inverse S^-1. One might be lucky to have a closed form of S^-1 available. For example, in case of the Galambos copula the function S^-1 is given by S^-1(t)=c_θ t^-θ with constant c_θ:=(Γ(d) Γ(1/θ)/Γ(d+1/θ)/θ)^-θ, cf. <cit.>. Interestingly, Algorithm <ref> in this particular case is different than the one derived in <cit.>, which is designed for extreme-value copulas rather than reciprocal Archimedean copulas, but which also includes the Galambos copula. The algorithm of <cit.> is always based on the decreasing sequence d/(ϵ_1+…+ϵ_k), but ours on R_k=S^-1(ϵ_1+…+ϵ_k) instead. For θ=1, these two sequences agree, so the simulation algorithms coincide[Except for a different simulation strategy of the uniform law on the simplex S_d.]. For θ≠ 1, however, they are truly different, since Algorithm <ref> always simulates the uniform law on the simplex and varies the sequence R_k, while <cit.> stick with the sequence d/(ϵ_1+…+ϵ_k) and instead vary the measure on the simplex.
Example <ref> shows that discrete measures ν give rise to S^-1 having (piecewise constant) closed form. As an example, scatter plots for the bivariate reciprocal Archimedean copula associated with the generator F(t)=exp(-Λ(t)), for Λ in (<ref>), are depicted in Figure <ref>. From our simulations, it appears as though the resulting family converges to the independence copula for θ→∞, and to some limiting copula (but not the upper Fréchet bound) for θ→ 0. Furthermore, since S^-1(t)=1/⌈ t/θ⌉ maps to the discrete set {1/k : k ∈}, the copula assigns positive mass to the one-dimensional subsets

    A_k:={(F(u/k),F(1-u/k)) : u ∈ (0,1)}⊂ [0,1]^2,   k ≥ 1,

of the unit square. By construction, this mass is decreasing in k, and the first sets A_1,A_2,… are clearly visible in the scatter plots.
  
 


In case S^-1 is not given in closed form, it is  convenient to recall that S is given in terms of Λ by the Williamson transform inversion formula 

    S(t) = ∑_k=0^d-2(-1)^k Λ^(k)(t)/k!t^k+(-1)^d-1 Λ^(d-1)_+(t)/(d-1)!(d-1)^k,

with Λ^(d-1)_+ denoting the right-hand derivative of Λ^(d-2), cf. <cit.>. In particular, Λ is d-monotone, which implies that the first d-2 derivatives exist, and Λ^(d-2) is convex (so that Λ^(d-1)_+ exists). If S has discontinuities, such as in the case of a discrete measure ν like in Example <ref>, Λ is not (d+1)-monotone, i.e. Λ^(d-1) does not exist. Besides traditional Newton-Raphson inversion, here are two more ideas for dealing with S^-1 in Algorithm <ref>: 


  (i) One idea to evaluate S^-1 approximatively is to approximate S by a piecewise constant function, which amounts to approximation of ν by a discrete measure. Example <ref> then gives the closed form of the resulting piecewise constant approximation to S^-1. This implies an approximative simulation algorithm for the reciprocal Archimedean copula in concern. In light of the relation between S and Λ in (<ref>), this means that the exponent Λ in concern is approximated by an exponent that is proper d-monotone, i.e. not (d+1)-monotone. Similarly, other invertible approximations might be feasible as well, for example a piecewise linear approximation of S. Notice that a piecewise constant approximation of a continuous survival function S results in an approximating reciprocal Archimedean copula with singular component, while a piecewise linear approximation would not have this drawback.

  (ii) We assume that S is differentiable. In step n+1 of the while-loop in Algorithm <ref>, the values x_n:=ϵ_1+…+ϵ_n and S^-1(x_n) have already been computed (in the previous step), and one seeks to simulate the random variable S^-1(x_n+ϵ_n+1). In the initial step we have x_0=0 and S^-1(0)=∞. The random variable S^-1(x_n+ϵ_n+1) has density

    f(x):=e^x_n e^-S(x) S^'(x),   x ∈(0,S^-1(x_n)),   n ≥ 0,

which is often given in closed form, e.g. by virtue of formula (<ref>), so can be evaluated efficiently. If one can find a random variable X, which one can simulate from and whose density g satisfies f ≤ c g for some c ≥ 1, traditional rejection acceptance sampling can be applied, cf. <cit.>, which results in an exact simulation algorithm. 



Using the stochastic representation (<ref>) for the uniform distribution on the simplex, the simulation of the random vector 𝐐^(k) in the k-th while-loop requires to simulate d independent exponential random variables, hence has complexity order linear in d. However, the number of required while-loops is random itself. In the special case ν(dx)=d x^-2 dx, which corresponds to the Galambos copula with parameter θ=1 and S^-1(x)=d/x, Algorithm <ref> coincides with <cit.>. <cit.>, which in turn refers to <cit.>, shows that the expected number of required while-loops in this case equals [M]=d [max{X_1,…,X_d}], where (X_1,…,X_d) is a random vector with survival copula the Galambos copula in concern and all one-dimensional margins unit exponentially distributed. It follows from a computation in <cit.> that

    1 ≤[max{X_1,…,X_d}] = ∑_i=1^ddi (-1)^i+1/H_i≤ d,

where H_n:=1+1/2+…+1/n, n=1,2,…, denotes the harmonic series. Hence, in this special case the complexity order of the algorithm is known explicitly as a function of the dimension d and lies somewhere between d^2 and d^3. Unfortunately, the proof of this result relies heavily on the fact that d/R_k=ϵ_1+…+ϵ_k. In the case of general radial measure, we have

    [M]     = 1+∑_m=1^∞(M>m)=1+∑_m=1^∞(R_k+1 > M_k, k=1,…,m)
        = 1+∑_m=1^∞( S^-1( ∑_i=1^m+1ϵ_i)/S^-1(∑_i=1^kϵ_i)> max{Q^(k)_1,…,Q^(k)_d }, k=1,…,m).
 
From this formula for S^-1(x)=d/x, <cit.> make use of the fact that ∑_i=1^kϵ_i/∑_i=1^m+1ϵ_i, for k=1,…,m, correspond to order statistics of independent samples from the uniform law on [0,1], which we cannot for general ν (hence S^-1). However, using this known Galambos case as benchmark, from the last formula we observe that [M] is larger than in the known benchmark case if the function x ↦ S^-1(x) x is increasing, and that it is smaller if x ↦ S^-1(x) x is decreasing. Although we cannot say a lot about the case when x ↦ S^-1(x) x is neither increasing nor decreasing, this at least provides a feeling for the effect of the choice of radial measure on the expected runtime of the algorithm. In the case when ν is absolutely continuous with density f_ν one may check whether x ↦ S^-1(x) x is increasing (resp. decreasing) by checking whether f_ν(x) x ≥ S(x)  (resp. f_ν(x) x≤ S(x)) for all x>0, which is a quite intuitive condition in terms of the density. Heuristically, it says that heavier tails of the radial measure make the algorithm faster, and vice versa.





§ CONCLUSION

An exact simulation algorithm for reciprocal Archimedean copulas has been presented. It was based on the concatenation of two ideas. On the one hand, via a change of variables transformation the points of a Poisson random measure on (0,∞], whose mean measure equals the radial measure of the reciprocal Archimedean copula, have been represented as a decreasing function of the jump times of a standard Poisson process. On the other hand, an idea of <cit.> has been enhanced from a simulation algorithm for extreme-value copulas to copulas of more general max-infinitely divisible distributions, of which reciprocal Archimedean copulas are a particular representative. 


[Bondesson (1982)]bondesson82 L. Bondesson, On simulation from infinitely divisible distributions, Advances in Applied Probability 14 (1982) pp. 855–869.
[Dombry et al. (2016)]dombry16 C. Dombry, S. Engelke, M. Oesting, Exact simulation of max-stable processes, Biometrika 103:2 (2016) pp. 303–317.
[Fang et al. (1990)]fang90 K.-T. Fang, S. Kotz, K.-W. Ng, Symmetric multivariate and related distributions, Chapman and Hall, London (1990).
[Galambos (1975)]galambos75 J. Galambos, Order statistics of samples from multivariate distributions, Journal of the American Statistical Association 70:351 (1975) pp. 674–680.
[Genest et al. (2018)]genest18 C. Genest, J. Nešlehová, L.-P. Rivest, The class of multivariate max-id copulas with L_1-norm symmetric exponent measures, Bernoulli, in press (2018).
[Mai (2014)]mai14 J.-F. Mai, A note on the Galambos copula and its associated Bernstein function, Dependence Modeling 2 (2014) pp. 22–29.
[Mai, Scherer (2017)]maischerer17 J.-F. Mai, M. Scherer, Simulating Copulas, 2nd edition, World Scientific Press (2017).
[McNeil, Nešlehová (2009)]neilnev09 A.J. McNeil, J. Nešlehová, Multivariate Archimedean copulas, d-monotone functions and l_1-norm symmetric distributions, Annals of Statistics 37:5B (2009) pp. 3059–3097.
[Nelsen (2006)]nelsen06 R.B. Nelsen, An introduction to copulas, 2nd edition, Springer (2006).
[Oesting et al. (2018)]oesting18 M. Oesting, M. Schlather, C. Zhou, Exact and fast simulation of max-stable processes on a compact set using the normalized spectral representation, Bernoulli 24:1 (2018), pp. 1497–1530.
[Resnick (1987)]resnick87 S. Resnick, Extreme values, regular variation and point processes, Springer-Verlag (1987).
[Schlather (2002)]schlather02 M. Schlather, Models for stationary max-stable random fields, Extremes 5:1 (2002), pp. 33–44.


