Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I am trying to numerically integrate an equation that involves InverseErfc (embedded in the copula defined). The equation looks like the following:

$$ \int_0^T \int_0^\infty \int_0^\infty e^{-Ab(s)x-Ar(s)y}\ f(x)\ c(F(x),G(y))\ g(y)\ dx\ dy\ ds\, $$

The function $\ f(x)$ and $\ g(x)$ are the density functions of $\ F(x)$ and $\ G(x)$ respectively, whereas the function $\ c(F(x),G(y))$ represents the copula chosen (in this case it's student-t copula). When combined, the student-t copula function with exponential marginals are called upon in MMA as GGE[Theta_,Alpha_, Beta_] in the program code.

(** Defining values of parameters & expressions **)
cr = 0.5; ar = -1; br = 0; Sigmar = 0.5; Gammar = 1; r0 = 0.0153; Betar = 100; 
cb = 0.1; ab = -1; bb = 0; Sigmab = 0.05; Gammab = 1; b0 = 0.04; Betab = 200;
Theta = 0.5; T = 1; Rho = 4; DoF = 3; 

compr = Sqrt[(cr ar)^2 + 2 (Sigmar Gammar)^2]; 
compb = Sqrt[(cb ab)^2 + 2 (Sigmab Gammab)^2];

Ar[s_] := (2 r0 (1 - Exp[-s compr]))/((compr - cr ar) + (compr + cr ar) Exp[-s compr]); 
Ab[s_] := (2 b0 (1 - Exp[-s compb]))/((compb - cb ab) + (compb + cb ab Exp[-s compb]);

(** Define the copula - Student-t copula with exponential margins **)
GEE[Theta_, Alpha_, Beta_] := CopulaDistribution[{"MultivariateT", {{1, Theta}, {Theta, 1}}, 
DoF}, {ExponentialDistribution[Alpha], ExponentialDistribution[Beta]}];

(** Obtaining the density function of the above copula and its Laplace transform **)
GumExpExp := PDF[GEE[Theta, Betab, Betar], {x, y}];
P[u_, v_] := LaplaceTransform[GumExpExp, {x, y}, {u, v}]

(** The following is the expression I would like to calculate **)
chat = NIntegrate[P[Ab[s], Ar[s]], {s, 0, T}]

I have been running the code for 8 hours and up till now there is no output. Can someone please point out what is the mistake or give me a different idea on how to obtain the numerical answer? Your help is much appreciated.

share|improve this question
There are a couple of missing parentheses and braces. In the definition of Ab[s_] check the last term. And also the missing ':' in the definition of GEE and also the missing braces. – Sektor Jun 23 '13 at 1:16
Start from the inside: your P[ ] function alone does not evaluate, so of curse you cannot integrate it. The reason P[ ] doesn't evaluate is because GumExpExp does not evaluate properly. The reason GumExpExp doesn'e evaluate is the Gee[ ] is not correct. This is where I gave up. Take things one step at a time, don't try to solve the whole problem at once. – bill s Jun 23 '13 at 4:52
Thanks @bills. Can you please explain on what you meant by "Gee[ ] is not correct". Which part of GEE is wrong (the syntax, or is it because its Laplace transform may not exist)? – SNMRamli Jun 23 '13 at 6:13
Try GumExpExp[x_, y_] := PDF[GEE[Theta, Betab, Betar], {x, y}] and then evaluate GumExpExp[x,y]. Do you get a function of x and y? – bill s Jun 23 '13 at 6:52
Is it necessary to write GumExpExp[x_,y_]? Because I still get a function of x and y without writing the [x_,y_]. – SNMRamli Jun 23 '13 at 7:29

LaplaceTransform attempts a symbolic evaluation of the transform, which Mathematica fails to do in this case. To get an answer, write the transform explicitly as an NIntegrate multiple integral.

It still takes a while, but reducing the accuracy/precision demands will speed things up:

chat = NIntegrate[
  GumExpExp Exp[-Ab[s] x - Ar[s] y],
  {x, 0, Infinity}, {y, 0, Infinity}, {s, 0, T},
  AccuracyGoal -> 5, PrecisionGoal -> 5] // AbsoluteTiming
  {27.454909, 0.999839}

Of course, you can adjust AccuracyGoal and PrecisionGoal to suit your needs and available time.

share|improve this answer

Keep in mind InverseErfc[2]=-infinity and InverseErfc[0]=infinity, so for InverseErfc function domain should be between 0 and 2, if not you did not got numeric answer.

share|improve this answer
This is an answer, but apparently not for this question. – J. M. Jun 23 '13 at 4:23
So should I use Integrate instead of NIntegrate? – SNMRamli Jun 23 '13 at 4:49

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.