You asked: "Why do you think that there's a solution?"
I am not sure that there's a solution given my conjectures as in equation (1) of the PDF file Problem_stylized.pdf. That is, I do think that equations (2) to (5) (and their translations into Maple) are correct, but equation (5) (i.e., the predictions, as simply derived from a normality assumption) could map to conjectures which are not linear.
Now I ask you:
Given the form of the beta's, alpha's and natlq2 and natlq3, can you think of polynomial conjectures with more reasonable degrees? Maybe quadratic? Square root? Any other power law?
Previosuly in this thread you wrote: "At least all the equations are algebraic functions (rational functions & fractional exponents) with the highest exponent being 4 and the only fractional exponent being 1/2." Is this still true for the latest attached scripts? Did you verify this with some Maple command? Most importantly, how to use this information to think of a more clever form for my conjectures?
I think is important to emphasize that I am trying to solve a two-steps problem. In the first step (not in this thread), I look for the beta's and alpha's as combinations of the lambda's and mu's. In the second step (this whole thread), I try to pin down the lambda's and mu's as expressions of the exogenous parameters and then plug them back into the beta's and alpha's to fully characterize my equilibrium.
That is, (in the first step) I use the linear conjecture for natlq2 and natlq3 (equation (1)) to find the equations for beta's and alpha's in the first place...therefore, if such conjecture for natlq2 and natlq3 (equation (1)) turns out to be non-linear, then also the equations for beta's and alpha's, which are those given at the beginning of the scripts, would need to be revised accordingly.
Finally, as you suggested, I attach here my attempts with numerical approaches.
With this calibration (calibration 1) I obtain some values for my lambda's and mu's but the numerical solutions look a bit strange, i.e. tiny. I run fsolve with the avoid option to check the variation of the solutions across the runs. For each variable, while the signs of the solutions seem to be preserved across runs, their values differ and they are all very tiny. What could it mean?
Here I use a more realistic calibration of my parameters (calibration 2). In this case fsolve fails, but I don't know why. What could it mean?
Finally, I use DirectSearch[SolveEquations]. Contrary to the root-finder fsolve(), this optimiser almost always returns some sort of least-square solution, as shown by the relatively minimal value of the sum of the squares of the residuals (the first two items in the list returned by the command). For calibration 2, I finally obtain some results but the caveat is that the outputs are not necessarily what we would ordinarily call solutions (is my understanding correct?). Moreover, the results I obtain for calibration1 differ from those I got using fsolve as I described in 1). All in all, what could it mean?
In general, how do I know whether there are numerical errors? Is there a way to verify the accuracy of the numerical results?
Again, I thank you for looking into this with such care.