# Numerical Integration with InverseErfc

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.

-
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.

-

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.

-
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