Skip to content

FSolve

Description

FSolve solves for the coefficients of Fleishman’s power method, which turns normal random numbers into ones with a requested skewness and excess kurtosis. It uses Newton’s method (up to 50 iterations) and returns the three coefficients, the number of iterations and the residual, so the caller can tell whether it converged. It doesn’t check whether the requested skewness and kurtosis are possible; see Fkmin.

Used by: MCM(_V), MCM_CI(_V), MCM_error(_V), rerror, rerrorV

Source Code

FSolve = LAMBDA(
    skew,       // requested skewness
    kurt,       // requested excess kurtosis (not checked against min required)
    [b], [c], [d], [n],     // for recursion
  LET(
    b0, IF(ISOMITTED(b), 1, b),
    c0, IF(ISOMITTED(c), skew/6, c),
    d0, IF(ISOMITTED(d), 0, d),
    n0, IF(ISOMITTED(n), 0, n),
    F_1, b0^2+6*b0*d0+2*c0^2+15*d0^2-1,
    F_2, 2*c0*(b0^2+24*b0*d0+105*d0^2+2)-skew,
    F_3, 24*(b0*d0+c0^2*(1+b0^2+28*b0*d0)+d0^2*(12+48*b0*d0+141*c0^2+225*d0^2))-kurt,
    resid, SQRT(F_1^2+F_2^2+F_3^2),
    IF(OR(n0 >= 50, resid < 0.0000001),
      HSTACK(b0, c0, d0, n0, resid),
      LET(
        J, VSTACK(
          HSTACK(2*b0+6*d0, 4*c0, 6*b0+30*d0),
          HSTACK(4*c0*b0+48*c0*d0, 2*b0^2+48*b0*d0+210*d0^2+4, 48*c0*b0+420*c0*d0),
          HSTACK(24*(d0+2*b0*c0^2+28*d0*c0^2+48*d0^3),
                 24*(2*c0+2*b0^2*c0+56*b0*d0*c0+282*c0*d0^2),
                 24*(b0+28*b0*c0^2+24*d0+144*b0*d0^2+282*c0^2*d0+900*d0^3))
        ),
        Fvec, VSTACK(F_1,F_2,F_3),
        delta, MMULT(MINVERSE(J), Fvec),
        FSolve(skew, kurt, b0-INDEX(delta,1,1), c0-INDEX(delta,2,1), d0-INDEX(delta,3,1), n0+1)
      )
    )
  )
);