To obtain the coefficients requires solving a system of nonlinear equations.

This is the system (below in Mathematica code) I used for third-order polynomials. There are boundary conditions on the skew (gamma1) and kurtosis (gamma2). For example, for symmetric distributions i.e. skew=0, the lower boundary of kurtosis is about -1.15 (I derived the boudary conditions in an article I published a few years ago, if you're interested I can get it for you).

Note: The mean and variance are set to 0 and 1. You can take a distribution and impose a linear tranformation on it and it will not change the values of the skew and kurtosis (or any other higher moments, for that matter).

Code:

```
gamma1 = 2;
gamma2 = 7;
FindRoot[{
a + c == 0,
b^2 + 6*b*d + 2*c^2 + 15*d^2 == 1,
2*c*(b^2 + 24*b*d + 105*d^2 + 2) == gamma1,
24*(b*d + c^2*(1 + b^2 + 28*b*d) + d^2*(12 + 48*b*d +
141*c^2 + 225*d^2)) == gamma2},
{a, -0.230}, {b, 0.88}, {c, 0.23}, {d, 0.019}, AccuracyGoal -> 24]
Solutions:{a -> -0.260023, b -> 0.761585, c -> 0.260023, d -> 0.0530723}
```

Other techniques that I mentioned such as g-and-h transformations and the Generalized Lambda distributions also require the simulatanoeus solutions to a system of nonlinear equations.

These techiques that I mention are especially useful when you need multivariate non-normal distributions with a specified correlation matrix.

For example, suppose you wanted X1, X2, and X3 all with different values of skew and kurtosis and with correlations of r12=0.4; r13=0.5; r23=0.6. These techniques I mention handle this scenario very easily.

Note: when I say non-normal I mean that the values of skew and kurtosis are not equal to zero. A normal distribution has skew of 0 and kurtosis of 0.