How can I generate lotteries where I can simultaneously determine their expected values, variances and skewness? First of all, thanks for taking the time to look at this. My goal is to be able to generate lotteries, where I can determine the first, second and third moments of their distributions, for a behavioural economics experiment. In other words, I'd like to be able to determine their expected values, variances and skewness.
\begin{align*}
 &\mathrm{Expected\ Value\ :} &&=\frac{\sum_{i}x_i}{n}\\
 &\mathrm{Variance:} &&=\frac{\sum_{i}(x_i-\mathrm{Mean})^2}{n}\\
 &\mathrm{Skewness:} &&=\sum_{i}\left( \frac{x_i-\mathrm{Mean}}{\sqrt{\mathrm{Variance}}}\right) ^3
\end{align*}
The lotteries will have three outcomes that are equally probable, so I believe this amounts to solving for $\{a,b,c\}$ in the following nonlinear system of equations, where $m, v$ and $s$ are just constants.
\begin{align}
 \frac{a+b+c}{3}= m\\
 \frac{
 (a-m)^2+(b-m)^2+(c-m)^2
 }{n}= v\\
 \frac{
 (a-m)^3+(b-m)^3+(c-m)^3
 }{\left( \sqrt{v}\right) ^3} = s
\end{align}
My questions is, how would you find for what values of $m,v$ and $s$ there are non-negative real solutions? I know that I can just input non-negative values for $a, b$ and $c$, but the issue with this is that, for my experiment, I would like to find triples of lotteries that satisfy one of the two following conditions. For three lotteries, say $\{1,2,3\}$, I would like either:
A. $\mathrm{Exp_1}=\mathrm{Exp_2}>\mathrm{Exp_3}$ and $\mathrm{Var_2}=\mathrm{Var_3}<\mathrm{Var_1}$ and $\mathrm{Skew_3}=\mathrm{Skew_1}<\mathrm{Skew_2}$
or
B. $\mathrm{Exp_1}=\mathrm{Exp_2}<\mathrm{Exp_3}$ and $\mathrm{Var_2}=\mathrm{Var_3}>\mathrm{Var_1}$ and $\mathrm{Skew_3}=\mathrm{Skew_1}>\mathrm{Skew_2}$
So, my question to you is, is this possible? And if so, how could I achieve it?
Many thanks and
best wishes,
Jeremy
 A: To me, it's more like a question on zeros of polynomials. It would be nice if someone provides a numerical perspective for it.
Your problem is equivalent to know if there are zeros for the following system of equation.
$$\begin{equation} a+b+c= u \\a^2+b^2+c^2= v \\ a^3+b^3+c^3= w\end{equation}$$
and that system is equivalent to:
$$\begin{equation} a+b+c= u=:x_1 \\ab+bc+ca= \dfrac{u^2-v}{2}=:x_2 \\ abc=\dfrac{2w+u^3-3uv}{6}=:x_3 \end{equation}$$
and this system of equations has three real roots( singular or not) , if the follow equation has three real roots.
$$X^3-x_1X^2+x_2X-x_3=0$$
Any further analysis that you want on $a,b,c$  can be translated to a study on the above equation.
A: The following Python code helps to find $a,b,c$ as a function of $m,v,s$ by transforming the set of equations eqs1 that define mean, variance and skewness into an equivalent set of equations eqs2 that is easy to solve.
In this case, it is enough to look at the last equation (which is printed by the code).
In the example I chose $m=v=s=1$ and the code produces the output
$$
a^3 - 3a^2 + \frac32a + \frac16
$$
which, by wolfram alpha, has the three solutions $-0.093199, 0.76963, 2.3236$.
By symmetry, the unique solution (up to swapping the values of $a,b,c$) is given by:
\begin{align*}
a &= -0.093199\\
b &=0.76963 \\
c &= 2.3236  
\end{align*}
Which means that you cannot realize $m=v=s=1$ if you want $a,b,c >0$.

from sympy import xring, QQ, lex
from sympy.polys import groebnertools as gt

P, variables= xring('c,b,a', QQ, lex)
c,b,a = variables

# choose your values for m,v,s below
m=1 # mean
v=1 # variance
s=1 # skewness

f1 = a+b+c - 3*m   
f2 = (a-m)**2 + (b-m)**2 + (c-m)**2 - 3*v  
f3 = (a-m)**3 + (b-m)**3 + (c-m)**3 - s*v**(3/2)

eqs1 = [f1,f2,f3]
eqs2 = gt.groebner(eqs1, P)

#for f in eqs2:
#    print(, "= 0")
print(eqs2[-1])


Read chapter 3.1 of Cox Little O'Shea if you want to know why this works. Note that they refer to  eqs2  as a Gröbner basis for the Ideal generated by  eqs1.
