# NSolve for large models

I'm a Mathematica rookie running version 7.0 for Mac OS X x86 (64-bit) (February 19, 2009) on a 2010 MacBook Pro. (I'm also a stackexchange "newbie" who originally posted this on stackoverflow, thus missing the target audience. I hope I'm not breaking any sacred rules by reposting the question here.)

Context: I'm trying to solve an economics problem with a large number of firms, where each firm is given a maximum production quota (i.e., rights to produce). The firms are of different sizes, larger ones being more efficient than smaller ones. There is a market for quotas, thus small, inefficient firms can shut down and sell their quota rights to larger ones; the quota price k must be computed endogenously. By and large, the bigger the number of firms, the bigger the number of equations to solve simultaneously.

The issue is as follows: with a small number of firms n (from 2 to 3), only Solve works, whereas NSolve erroneously yields {}, i.e., no solution. With a medium number of firms (from 4 to 7), only NSolve works. Solve does not seem to work for any problem with 4 or more firms. With 8 and 9 firms, NSolve again erroneously yields {}. With 10 to 16 firms, NSolves works fine. With a larger number of firms (17 and above), the evaluation keeps running without end and I end up aborting (I waited up to 15 minutes).

My hope was to run the model with thousands of firms. My questions are
(1) how can I systematically get Mathematica to converge to the solution for an arbitrary large system and
(2) why does NSolve not work for small systems.

Two remarks:

(1) You will notice some awkwardness in the way I set up the model, starting with a guess k0 for k, which determines which firms are operating and which ones are shutting down. Upon that, for the list of operating firms, I compute the profit maximizing amount of labor, total output, and equilibrium k. If k is different from k0, I restart the computation with k0=k. The question of how to create a more elegant one-step program is a separate issue, though, and I don't think it is connected to what appears to be a model size limitation.

(2) If I solve the model for, say, 11 firms and then jump trying to solve for, say, 13 firms, Mathematica ends up running "endlessly". But if I solve for 12 firms first, then solving for 13 firms is a matter of seconds. So it would seem that despite the ClearAll["Global*"] statement at the beginning of the notebook, it is as if Mathematica seemed to start with a better starting point.

I'm ready to be humbled :)

UPDATE 1: I'm now running version 10.0 and NSolve works fine with all cases involving less than 16 firms, which generate a system of at most 4 equations in 4 unknowns. The case with 17 firms generates a system with 5 equations and 5 unknowns, which appears to be too much for NSolve. I changed the question title accordingly, to focus on NSolve only.

UPDATE 2: If we use the case with 12 firms as an illustration, the generated set of equations is $\left\{ \begin{array}{l} \frac{5154.55 (5-k)}{\text{tl}_1^{0.3}}-100=0 \\ \frac{5727.27 (5-k)}{\text{tl}_2^{0.3}}-100=0 \\ \frac{6300 (5-k)}{\text{tl}_3^{0.3}}-100=0 \\ 7363.64 \text{tl}_1^{0.7}+8181.82 \text{tl}_2^{0.7}+9000 \text{tl}_3^{0.7}=1000000 \end{array} \right.$

which can be rewritten

$\left\{ \begin{array}{l} \text{tl}_1^{0.7} = (5-k)^{7/3} 51.5455^{7/3} \\ \text{tl}_2^{0.7} = (5-k)^{7/3} 57.2727^{7/3} \\ \text{tl}_3^{0.7} = (5-k)^{7/3} 63^{7/3} \\ (5-k)^{7/3} \left( 7363.64 \times 51.5455^{7/3} +8181.82 \times 57.2727^{7/3} +9000 \times 63^{7/3} \right)=1000000 \end{array} \right.$,

a problem simple enough to be solved with a pocket calculator. With more firms, i.e., more equations and more unknowns, the structure should be just as simple; it is therefore a mystery to me why NSolve appears to be overwhelmed.

Here is the code:

\$Version
ClearAll["Global*"]
(* Aggregate Quota Allocation *)
alloc=1000000;
(* Number of firms, variable that I would like to increase arbitrarily; *)
nv=4;
(* Equal quota allocation to each firm; *)
qbar = alloc/nv;
(* Firm sizes vary in amount of capital from smallest to largest;*)
vsmin=0 ;
vsmax=900000 ;
(* Assume discrete uniform distribution of firm sizes;*)
dvs=(vsmax-vsmin)/(nv-1) ;
(* List all firms by size; *)
listvs=Table[vs,{vs,vsmin,vsmax,dvs}]  ;
(* A firm's capital cost per annum is 10 percent of capital size;*)
g=0.1 ;
(* Avoidable fixed cost of operation;*)
f=5000 ;
(* Parameters of production function (same for all firms);*)
a=0.01 ;
alpha=0.7 ;
(* Price of output; *)
p=5;
(* Price of labor (daily wage);*)
wa=100 ;
(* Firm vs is a profit maximizer setting optimal number of days of operation tl[vs] ;
Firms that are too small are better off leasing their quota allocation to larger ones;
The quota leasing price k is computed endogenously;
Firms will exit if "netprofitstay" is less than 0;
We are looking for threshold firm, whose size nvs1 is just above threshold size vs0;
All firms whose size is nvs1 or larger are operating, the others are shutting down;
For now, we are setting the problem up as a manual loop, setting k0 = "k opt"
from previous computation ; *)

k0=4.891714989052472 ;

profitstay0[tl0_,vs0_]:= (p-k0)(a vs0 tl0^alpha) - wa tl0 + k0 qbar - g vs0 - f  ;
netprofitstay0[tl0_,vs0_]:=(p-k0)(a vs0 tl0^alpha) - wa tl0 - f  ;
sol0:=NSolve[{D[profitstay0[tl0,vs0],tl0]==0 ,netprofitstay0[tl0,vs0]==0} ,{tl0,vs0}]  ;
vs1:=vs0/.Last[sol0]
nvs0:=Part[Nearest[listvs,vs1],1] ;
nvs1:=If[nvs0<vs1,nvs0+dvs,nvs0] ;
(* List of firms who operate when k=k0; *)
listvsop:=Table[vs,{vs,nvs1,vsmax,dvs}]  ;

(* Production function;*)
eq[tl[vs_],vs_]:=a vs tl[vs]^alpha ;
(* Profit; *)
profitstay[k_,tl[vs_],vs_]:=(p-k) eq[tl[vs],vs] - wa tl[vs] + k qbar - g vs - f  ;
(* List of first-order conditions for a maximum; *)
foc:=Table[D[profitstay[k,tl[vs],vs],tl[vs]]==0,{vs,listvsop}]  ;
(* List of all endogenous variables; *)
variables:=Flatten[{k,Table[tl[vs],{vs,listvsop}]},1] ;
Print[variables] ;
(* List of all equations, including the one specifying that aggregate quota demand = aggregate
quota supply; *)
equations:=Flatten[{foc,Sum[eq[tl[vs],vs],{vs,listvsop}]==alloc }] ;
Print[equations] ;

(* Find solution;*)
NSolve[equations,variables]
`
-
Update: I upgraded to version 10. As with version 7, the code as posted here runs fine. In fact, NSolve now also works fine for any number of firms between 2 and 16, solving the problem in matters of seconds. For 17 firms or more, however, the evaluation still keeps running without end in sight. –  samdak Jul 26 at 12:12
try working with findroot. your question needs a better title by the way (what does Solve have to do with anything?) –  george2079 Jul 26 at 13:23
Thank you @george2079 for taking a look. I had Solve in the question because, at first, it worked for the small system case (2 or 3 firms), whereas NSolve did not. Now that I'm running version 10, NSolve does indeed also work for the small system case. I'm hesitant to use FindRoot because (1) it requires an initial guess and (2) the generated equations to solve are so simple I can solve them by hand. I therefore suspect that I'm missing something quite basic in the way I'm setting up the notebook. Thanks again. –  samdak Jul 26 at 20:46