I need to solve, numerically, a very large system of equations (Logit-type). The system has one parameter $d$ which is numerically defined in order to find the solution.
My code works, but it requires 11 minutes to run. I want:
- To increase the efficiency of my code (I have to estimate also a different "pay" function which will generate around 10.000 equations instead of the actual 250).
- I want Mathematica to run my program across many values of $d$ (user-defined).
- I want Mathematica to save solutions for each $d$ in a dataset which can be exported to Stata for a maximum likelihood estimation (my task will be to find the parameter $d$ which better fit my data). Or, maybe, I can import my observations in Mathematica for the maximum likelihood estimation? What's best?
Here is the code:
ClearAll[qmax, qmax2, pay]
qmax[h_] := Floor[8 Sqrt[2] Sqrt[h]];
qmax2[h_, q_] := Min[qmax[h], (q - 1)];
pay[h_, q_] := (32 - q^2/(4 h))*(1/10)*Sum[
Sum[Symbol["t" <> ToString[t] <> "q" <> ToString[s]], {s, 0, qmax2[t, q]}]+
If[q <= qmax[t], 0.5*Symbol["t" <> ToString[t] <> "q" <> ToString[q]], 0],
{t, 1, 10}];
equations = Flatten@Table[
Table[
Symbol["t" <> ToString[h] <> "q" <> ToString[q]] ==
Exp[d*pay[h,q]]/Sum[Exp[d*pay[h, s]], {s, 0, qmax[h]}],
{q, 0, qmax[h]}],
{h, 1, 10}];
vars = Flatten[Table[
{Symbol["t" <> ToString[h] <> "q" <> ToString[q]], 0.01},
{h,1, 10}, {q, 0, qmax[h]}], 1];
d = 1;
FindRoot[equations, vars]
FindRoot
is to start with a good guess. Do you know more about the solution than 0.01 for each component? When you do get to varying the parameterd
, try using the previousd
's answer as an initial guess for the next root -- or get more sophisticated with linear or quadratic extrapolation. $\endgroup$equations = Total[Abs[Flatten@Table[Table[Symbol["t"<>ToString[h]<>"q" <>ToString[q]]*Sum[Exp[d*pay[h, s]], {s,0,qmax[h]}]-Exp[d*pay[h,q]], {q,0,qmax[h]}], {h,1,10}]]] //. {0.5*xxx_+0.5*yyy_->0.5*(xxx+yyy), 2*xxx_+2*yyy_->2*(xxx+yyy)}; vars=Flatten[Table[{Symbol["t"<>ToString[h]<>"q"<>ToString[q]]}, {h,1,10}, {q,0,qmax[h]}]]; d=1; AbsoluteTiming[NMinimize[equations, vars]]
takes about 30 seconds to find an approximate solution. That only times the solve. Can you confirm that this works? Can you show a problem with 1000 equations? Can you give a list of example d values? $\endgroup$