Numerically stable option for the stationary distribution and its gradient

[jlperla]

The current "regularization" hack to deal with lower-rank stationary distributions could be replaced with something more stable. That has two parts:

• [ ] Add an option which will do some sort of eigenvalue decomposition with the output of the lyapdinstead of regularization directly or taking the cholesky https://github.com/HighDimensionalEconLab/DifferentiableStateSpaceModels.jl/blob/v0.4.14/src/generate_perturbation.jl#L250-L251
  • This would probably work by setting all small negative values to be 0, which makes it positive semi-deifnite and low-rank
  • At that point, the it could be multiplied back to create the matrix, and then take the cholesky, or maybe the eigenvalue decomposition or QR could be manipulated directly.
• [ ] That probably work work reasonably well. The bigger concern is what happens in the derivative of that matrix wrt to the parameters when the matrix itself is lower-rank? Is https://github.com/HighDimensionalEconLab/DifferentiableStateSpaceModels.jl/blob/v0.4.14/src/generate_perturbation_derivatives.jl#L160 still the right calculation or does it become numerically unstable/etc.

[jlperla]

Doing this with a pivoted cholesky instead