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.

  • $\begingroup$ 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. $\endgroup$ – Sektor Jun 23 '13 at 1:16
  • $\begingroup$ 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. $\endgroup$ – bill s Jun 23 '13 at 4:52
  • $\begingroup$ 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)? $\endgroup$ – SNMRamli Jun 23 '13 at 6:13
  • $\begingroup$ 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? $\endgroup$ – bill s Jun 23 '13 at 6:52
  • $\begingroup$ Is it necessary to write GumExpExp[x_,y_]? Because I still get a function of x and y without writing the [x_,y_]. $\endgroup$ – 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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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