Include parameter-dependence in equilibrium moments calculation

This is not currently implemented in FPP, but I think it would be nice

It is entirely dependent on if the tracking code can give us a stochastic matrix as a function of parameters ( tpsa) instead of scalars, which I think in Julia should be doable if we write the tracking code properly.

Equilibrium moments calculation we would have to:
1. Go to parameter-dependent fixed point to all orders (obtained from factorization of normal form transformation)
2. At the (fully-nonlinear) fixed point get linear normalizing matrix
3. Same calculation as before