# nested numerical integration: not valid limits?

I encountered an error when I was hoping for some Mathematica (8.0.4) magic sparing me to code up numerical integration, function approximation and root-finding myself. The broader context and the reference functions here come from a previous question:

k=10/3;
H = ParetoDistribution[1.18709*10^6, 0.938482];
Hstar = H;
u[c_,l_]=Log[c]-Log[1+l^(1+k)/(1+k)];
T[z_]=(1-0.84/1.3) * z;
lType[n_]=ArgMax[{u[n l-T[n l],l],l>=0},l];
zType[n_] = n lType[n];


where u is the utility function as in Saez 2001 allowing for income effects, T the actual tax schedule (approximate). Note that lType and zType only makes sense numerically. The original problem with

Type=InverseFunction[zType]


was solved with not SetDelayed but Set everywhere. But the following nested (numerical) integration does not work (for any H distribution you give it). Have I defined it wrong?

g[z_]:=NSolve[T'[z]/(1-T'[z])-(k SurvivalFunction[H,z]/(z PDF[Hstar,z]))
Integrate[(1-g)Exp[Integrate[1-\[Xi]u[zzz]/\[Xi]c[zzz],{zzz,z,zz}]]
PDF[H,zz]/SurvivalFunction[H,z],{zz,z,\[Infinity]}],g,Reals];
N[g[2000000], 10]


Error:

NIntegrate::nlim: zzz = zz is not a valid limit of integration.

EDIT: More interestingly, when I simply specify numeric arguments for the inner integral, the computation starts, though does not finish in an hour:

  tmp[z_?NumericQ, zz_?NumericQ] := NIntegrate[1 - \[Xi]u[zzz]/\[Xi]c[zzz], {zzz, z, zz}]
g[z_] := NSolve[T'[z]/(1 - T'[z]) - (k SurvivalFunction[H,z]/(z PDF[Hstar, z]))
NIntegrate[(1 - g) Exp[tmp[z, zz]] PDF[H, zz]/SurvivalFunction[H, z], {zz,z,
\[Infinity]}], g, Reals]


But the same gives me an error if I compile first. Why? Because the limit at infinity? How should I proceed?

    tmp = Compile[{{z, _Real}, {zz, _Real}},NIntegrate[1 - \[Xi]u[zzz]/\[Xi]c[zzz],
{zzz, z, zz}],Parallelization -> True, CompilationTarget -> "C"]


Then I run g[2000000] and get the error:

CompiledFunction::cfsa: "Argument zz at position 2 should be a !(\"machine-size real number\")."

• If you reformat the code from the link so that people can copy and paste to their notebook, you are more likely to get attention from folks here. As it is, the In's and Out's are mixed in together with comments. Just a suggestion. Commented May 18, 2012 at 21:06
• @MichaelWijaya, thanks, I am happy to improve this, of course. But how do you do what is needed? I tried to all "Copy As" variants of Mathematica, and none were both human-readable and Mathematica-ready. And my preceding declarations are needed for a run anyway, no? So paste many-many lines? Commented May 18, 2012 at 21:10
• @MichaelWijaya: and by the way, I don't think I mixed the outs with the ins. The code there is what it should be. Yes, I pasted then some of the Out as some text of mine. Why is that a problem? Commented May 18, 2012 at 21:11
• Let me try to edit your question to get you started, and you can make changes later. Commented May 18, 2012 at 21:11
• Done with the edits. Please feel free to revert to the original post or make changes. I really do not mean to be difficult. Commented May 18, 2012 at 21:21

I'm putting this in an answer because I need a bit of room.

I don't see the $\xi u$ and $\xi c$ functions defined anywhere. Also you've defined g in terms of itself. g[z_]:= NSolve[FUNC[g],g] which doesn't make sense to me. NSolve will return a rule like g -> something.

The error message is given by NIntegrate which cannot take a variable in the limits, i.e.

  NIntegrate[f[x],{x,0,b}] gives an error.


MMa seems to have conniptions and default to NIntegrate when the function in the Piecewise statement has finite precision. Compare:

Integrate[Piecewise[{{y + 1, y >= 1}}, 0]*Integrate[\[Phi][x], {x, 0, y}], {y, 1, Infinity}]

Integrate[Piecewise[{{y+ 1.1, y >= 1}}, 0]*Integrate[\[Phi][x], {x, 0, y}], {y, 1, Infinity}]  which gives the error.


My suggestion is to try to break up your g[z] function into workable parts. For instance maybe define your inside integral as a function and then integrate that separately. It's easier to check that you're getting the proper behavior from each part if everything isn't mushed into one statement.

• thanks a lot. I'll add the [Xi]'s needed for tmp in an edit. For the definition of g: I think I did not confuse the definition of the function with a variable used within it, but maybe this more than just inelegant. Let me also ask back: Why is NIntegrate a problem, then? And how else would you pick this apart if not the way I tried to pre-define or pre-compile tmp above? (Which I still don't get why it does not work, by the way.) Commented May 19, 2012 at 14:12

With some dummy definitions and small changes, I think I've got a working equivalent of what you might have intended.

There was one major edit: changing the upper limit of your outer NIntegrate expression from Infinity to a large positive number. This seemed necessary because the limit of the corresponding integrand seems to approach 1/0 as the variable of integration approaches Infinity.

I also used FindRoot instead of NSolve. Anyway, please make any changes necessary, and I hope this helps.

k = 10/3;
H = ParetoDistribution[1.18709*10^6, 0.938482];
Hstar = H;
u[c_, l_] := Log[c] - Log[1 + l^(1 + k)/(1 + k)];
T[z_] := (1 - 0.84/1.3)*z;

\[Xi]u[z_] := z
\[Xi]c[z_] := z + 1

int1[l_?NumericQ, u_?NumericQ] :=
NIntegrate[1 - \[Xi]u[zzz]/\[Xi]c[zzz], {zzz, l, u}]

int1[2000000, 10000000]

(* Out: 1.60944 *)

Limit[SurvivalFunction[H, z], z -> Infinity]

(* Out: 0. *)

int2[z_?NumericQ, gg_?NumericQ] :=
NIntegrate[(1 - gg) Exp[int1[z, zz]] PDF[H, zz]/
SurvivalFunction[H, z], {zz, z, 10^9}]

int2[2000000, 1/2]

(* Out: 3.55197 *)

sys[z_?NumericQ, gg_?NumericQ] :=
T'[z]/(1 - T'[z]) -
k SurvivalFunction[H, z]/(z PDF[Hstar, z]) int2[z, gg]

g[z_] := Block[{ggg}, FindRoot[sys[z, ggg], {ggg, -0.1}]]

g[2000000] // AbsoluteTiming

(* Out: {6.5703641, {ggg -> 0.978297}} *)

• William: Thanks a lot. Embarrassingly enough, I stupidly overlooked a crucial part, that the whole function g should be used in the integral, so it is ill-defined to hope for a solution for a single g[z] value anyway. I reposted a question with a much-simplified special case, but the integral equation, hoping to use DSolve: mathematica.stackexchange.com/q/5840/1273 Commented May 20, 2012 at 19:20