2
$\begingroup$

I'd like to evaluate the following numerical integration using Mathematica

$$ \ \int_0^T\int_0^\infty xe^{-0.04 s}g(x,s) dxds\ $$

where g(x,s) is a Gaussian copula function with say, marginal exponential distributions with parameter A and parameter B (i.e. X~exp(A) and S~exp(B)). In the following MMA attempt, the g(x,s) is defined by GE[Theta,A,B].

This is what I've tried so far:

GE[Theta_, A_, B_] := CopulaDistribution[{"Binormal", Theta}, {ExponentialDistribution[A],ExponentialDistribution[B]}];

Delta = 0.04; A = 0.10; B = 10; Theta = 0.9; T=5;
GExpExp[x_,s_] := PDF[GE[Theta, A, B], {x, s}]

NIntegrate[Exp[-Delta*s]*x GExpExp, {x, 0, Infinity}, {s, 0, T}]

and the following is the error I obtain:

*NIntegrate::inumr: The integrand 1.1547 E^(-10 s-0.1 x-s [Delta]+2 InverseErfc[2 \ Plus[<<2>>]]^2-1.33333 (0.707107 InverseErfc[<<1>>]-Power[<<2>>] \ InverseErfc[<<1>>])^2) x has evaluated to non-numerical values for all sampling points in the region with boundaries {{0.05,36},{0.5,1}}. >>*

I then try to integrate the expression with respect to only one variable, in this case I chose variable x, so that I'd end up with a univariate expression (and hence I might be able to work on a simpler problem from then on)

  Integrate[Exp[-Delta*s]*GExpExp1[x, s], {x, 0, Infinity}]

  (* 11.547 E^(-10.04 s + 1.33333 InverseErfc[4 - 4 E^(-10 s)]^2) *)

It seems that eventually I still have to deal with the Inverse[NormalDistribution] function, which I think is the source of this issue. I then continued with the following: 

  T = 5; 
  Assuming[s > 0, LaplaceTransform[NIntegrate[11.547005383792516` E^(-10.04` s + 1.3333333333333337` InverseErfc[4 - 4 E^(-10 s)]^2), {s, 0, T}], s, T]]

  (* NIntegrate::inumr: The integrand 11.547 E^(-10.04 s+1.33333 InverseErfc[4+Times[<<2>>]]^2) has evaluated to non-numerical values for all sampling points in the region with boundaries {{1.,2}}. >>  *)

It could be because the function does not converge, but even after reducing the limit interval from [0,T] to [2,3] (and of course changing T to 5 instead of 2) I still can't obtain a number. Does anyone know the correct syntax that I should use so that I can obtain a numerical answer eventually? Your assistance is much appreciated.

Improve this question
$\endgroup$
8
  • $\begingroup$ Certainly, T is not a number, and you do not seem to have assigned a value to it... $\endgroup$
    – J. M.'s eventual burnout ♦
    Jun 7, 2013 at 19:31
  • $\begingroup$ Instead of GExpExp = PDF[GE[\[Theta], A, B, {x, s}] .. use: GExpExp = PDF[GE[Theta, A, B], {x, s}] ... note the different placement of the closing brackets ] ... and then you will need to choose a value for T too, or set up K[T] := NIntegrate[ blah] .... not .... K[T] = NIntegrate[ ...] $\endgroup$
    – wolfies
    Jun 7, 2013 at 19:34
  • $\begingroup$ You still haven't addressed the thing @wolfies was pointing out. In any event, K is a built-in reserved symbol, so please do not use it to denote your function. $\endgroup$
    – J. M.'s eventual burnout ♦
    Jun 8, 2013 at 4:17
  • $\begingroup$ Thanks for pointing out the errors guys. I have edited the question by including the value of parameter T as well as have written the correct expression of PDF[GE[Theta,A,B]] like how it appears in my code. $\endgroup$
    – SNMRamli
    Jun 8, 2013 at 4:24
  • 1
    $\begingroup$ There does indeed appear to be a problem here now (after everything has been fixed up) ... mma 9 hangs or kernel crashes on my Mac. You can force it to produce an answer by adding an option such as NIntegrate[blah, Method -> "MonteCarlo"] ... but, in any event, it would appear that the integral does not converge. $\endgroup$
    – wolfies
    Jun 8, 2013 at 8:58

1 Answer 1

Reset to default
1
$\begingroup$

If I use exact coefficients, I get an exact answer with Integrate after a couple of minutes:

GE[Theta_, A_, B_] := 
  CopulaDistribution[{"Binormal", Theta}, {ExponentialDistribution[A],
     ExponentialDistribution[B]}];

Delta = 4/100; A = 10/100; B = 10; Theta = 90/100; T = 5;
GExpExp[x_, s_] := PDF[GE[Theta, A, B], {x, s}]

Integrate[Exp[-Delta*s]*GExpExp[x, s], {x, 0, Infinity}, {s, 0, T}]
(*
  (250 (-1 + E^(251/5)))/(251 E^(251/5))
*)
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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