fix: risk_metrics.summary_statistics weighted skewness and kurtosis use population central moments without bias correction

[AlexFiliakov]

Problem Statement

risk_metrics.py summary_statistics() (lines 718-725) computes weighted skewness and kurtosis using population central moments (m3, m4) divided by population std cubed/quarted. Both the central moments and the std use population (biased) estimators. This systematically underestimates the magnitude of both statistics.

v1.0 Impact

High. Weighted skewness and kurtosis are used to characterize loss distributions in risk reporting. Underestimated kurtosis (heavy-tail measure) can lead to underestimation of tail risk, which is the core concern of insurance optimization.

Affected Code

ergodic_insurance/risk_metrics.py, lines 716-725:

else:
    mean = np.average(self.losses, weights=self.weights)
    variance = np.average((self.losses - mean) ** 2, weights=self.weights)  # population
    std = np.sqrt(variance)
    # Weighted skewness
    m3 = np.average((self.losses - mean) ** 3, weights=self.weights)  # population 3rd moment
    skew = m3 / (std**3) if std > 0 else 0
    # Weighted kurtosis
    m4 = np.average((self.losses - mean) ** 4, weights=self.weights)  # population 4th moment
    kurt = (m4 / (std**4) - 3) if std > 0 else 0

Current Behavior

Uses population estimators:

• variance = sum(w_i * (x_i - mean)^2) / sum(w_i)
• m3 = sum(w_i * (x_i - mean)^3) / sum(w_i)
• skewness = m3 / variance^(3/2)

Expected Behavior

Apply bias corrections analogous to the unweighted N/(N-1) and N^2/((N-1)(N-2)) corrections, using V1 and V2 (sum of weights, sum of squared weights):

For the variance, the Bessel correction is: var_corrected = var_pop * V1^2 / (V1^2 - V2) (already implemented in risk_adjusted_metrics lines 636-643).

For weighted skewness, the correction factor is:

V1 = np.sum(weights)
V2 = np.sum(weights**2)
V3 = np.sum(weights**3)
skew_corrected = m3 / var_corrected**(3/2) * V1**3 / (V1**3 - 3*V1*V2 + 2*V3)

Alternative Solutions

1. Apply weighted bias corrections for m3 and m4 (most correct)
2. Use the existing corrected variance from risk_adjusted_metrics and derive corrected higher moments
3. Use scipy.stats with frequency weights (if applicable)

Recommended Approach

Option 1: apply the standard weighted cumulant corrections. These are well-documented for reliability weights.

Acceptance Criteria

• Weighted skewness uses bias-corrected formula
• Weighted kurtosis uses bias-corrected formula
• For equal weights, results match scipy.stats.skew/kurtosis with bias=False
• Weighted variance uses same correction as risk_adjusted_metrics (consistency)

References

• Rimoldini, L. (2014). Weighted mean, variance, covariance, and correlation. Astronomy and Computing, 5, 1-8.
• Wikipedia: Weighted arithmetic mean, section on weighted sample variance and higher moments.
• Related: fix: risk_metrics.py parametric VaR uses population variance for weighted case (no Bessel correction) #940/Math: parametric VaR uses population variance for weighted case (no Bessel correction) #685 (weighted VaR population variance), fix: risk_metrics.summary_statistics uses ddof=0 for unweighted std (line 657) #942/fix: RiskMetrics.summary_statistics uses population std (ddof=0) #508 (unweighted ddof=0 in same method)