nested numerical integration: not valid limits?

Question

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\")."

Answer 1

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.

Answer 2

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}} *)