quantiles of sums using copula distributions too slow

Question

Trying to create a table for quantiles of the sum of two dependent random variables using built-in copula distributions (Clayton, Frank, Gumbel) with Beta marginals. Tried NProbability and FindRoot with various methods -- not fast enough. An example of the copula-marginal combinations I need to explore is the following:

nProbClayton[t_?NumericQ, c_?NumericQ] := 
        NProbability[  x + y <= t, {x, y}  \[Distributed]    
               CopulaDistribution[{"Clayton", c}, {BetaDistribution[8, 2], 
                                                   BetaDistribution[8, 2]}]]

For a single evaluation of the numeric probability using

nProbClayton[1.9, 1/10] // Timing // Quiet

I get

{4.914, 0.939718}

on a Vista 64bit Core2 Duo T9600 2.80GHz machine (MMA 8.0.4)

To get a quantile of the sum, using

FindRoot[nProbClayton[q, 1/10] == 1/100, {q, 1, 0, 2}// Timing // Quiet

with various methods

( `Method -> Automatic`, `Method -> "Brent"`, `Method -> "Secant"` ) 

takes about a minute to find a single quantile: Timings are

{48.781, {q -> 0.918646}}
{50.045, {q -> 0.918646}}
{65.396, {q -> 0.918646}}

For other copula-marginal combinations timings are marginally better.

Need: any tricks/methods to improve timings.

Albatross

Answer 1

The CDF of a Clayton-Pareto copula with parameter c can be calculated according to

cdf[c_] := Module[{c1 = CDF[BetaDistribution[8, 2]]}, 
   (c1[#1]^(-1/c) + c1[#2]^(-1/c) - 1)^(-c) &]

Then, cdf[c][t1,t2] is the probability that x<=t1 and y<=t2. This means that you can calculate the probability that x+y<=t according to

prob[t_?NumericQ, c_?NumericQ] := 
   NIntegrate[Derivative[1, 0][cdf[c]][x, t - x], {x, 0, t}]

The timings I get on my machine are

prob[1.9, .1] // Timing

(* ==> {0.087518, 0.939825} *)

Note that I get a different value for the probability from the one in the original post. However, running nProbClayton[1.9,0.1] produces a warning about slow convergence which could mean that the result in the original post is off. Also, if I change x+y<=t to x+y>t in the original definition of nProbClayton and calculate 1-nProbClayton[1.9,0.1] I get 0.939825 (without warnings) which is the same result as above.

For the quantile of the sum I get

FindRoot[prob[q, .1] == .01, {q, 1, 0, 2}] // Timing

(* ==> {1.19123, {q -> 0.912486}} *)

Again, I get a different result from the one in the original post but similar to before, changing x+y<=t to x+y>t and calculating FindRoot[nProbClayton[q, 1/10] == 1-1/100, {q, 1, 0, 2}] returns the same value for q as above.

Defraud

Thank you so much! This is the kind of method I was hoping for. I get similar timings and the same final result for q, but prob[.1, 1.9] // Timing gives (* ==> {0.1089999,1.6476002080383803*^-10}*)` on my machine. – Albatross

@kguler that should have been prob[1.9, 0.1]. I had c and t reversed in a previous version but I'd forgotten to change it in my answer. I've fixed it now. – Defraud

A minor modification of cdf to cdf[c_] := Module[{c1 = CDF[BetaDistribution[8, 2]],c2 = CDF[BetaDistribution[8, 2]]}, (c1[#1]^(-1/c) + c2[#2]^(-1/c) - 1)^(-c) &] works as a template for arbitrary bivariate distributions. Regarding <= in original post, yes, your suggested change or simplychanging <= to < in nProbClayton gives 0.939825 and no slwcon warnings. – Albatross