Moments and Skewness¶

Moments describe the shape of a distribution, with mean (1st moment), variance (2nd), skewness (3rd), and kurtosis (4th) quantifying center, spread, asymmetry, and tail heaviness respectively.

Moments Definition¶

1. Raw Moments¶

E[X^n]:

from scipy import stats
import numpy as np

data = np.random.normal(0, 1, 1000)

# First raw moment = mean
m1 = np.mean(data)
print(f"1st moment (mean): {m1:.3f}")

# Second raw moment
m2 = np.mean(data**2)
print(f"2nd moment: {m2:.3f}")

# Third raw moment
m3 = np.mean(data**3)
print(f"3rd moment: {m3:.3f}")

2. Central Moments¶

E[(X - μ)^n]:

mean = np.mean(data)

# Second central moment = variance
c2 = np.mean((data - mean)**2)
print(f"Variance: {c2:.3f}")

# Third central moment
c3 = np.mean((data - mean)**3)
print(f"3rd central moment: {c3:.3f}")

# Using scipy
print(stats.moment(data, moment=2))  # Variance
print(stats.moment(data, moment=3))  # 3rd central

3. Standardized Moments¶

E[(X - μ)^n / σ^n]:

std = np.std(data, ddof=1)

# Skewness (3rd standardized moment)
skew = np.mean(((data - mean) / std)**3)
print(f"Skewness: {skew:.3f}")

# Using scipy
print(stats.skew(data, bias=False))

4. Relationship¶

# Raw → Central
# μ₂' = E[X²] = Var(X) + (E[X])²
# μ₃' = E[X³] = μ₃ + 3μ(Var) + μ³

# Central → Standardized
# γ₃ = μ₃ / σ³ (skewness)
# γ₄ = μ₄ / σ⁴ - 3 (excess kurtosis)

5. Moment Generating Function¶

# MGF: M(t) = E[e^(tX)]
# Derivatives at t=0 give moments

# For normal distribution
t = 0.5
norm_mgf = np.exp(mean * t + 0.5 * std**2 * t**2)
print(f"Normal MGF at t={t}: {norm_mgf:.3f}")

6. Cumulants¶

# Cumulant generating function
# K(t) = log(M(t))

# First cumulant = mean
# Second cumulant = variance
# Third cumulant = third central moment

7. Sample vs Population¶

# Sample moments (unbiased)
sample_var = np.var(data, ddof=1)  # Divide by n-1

# Population moments (biased)
pop_var = np.var(data, ddof=0)     # Divide by n

print(f"Sample variance: {sample_var:.3f}")
print(f"Population variance: {pop_var:.3f}")

Skewness¶

1. Definition¶

Measure of asymmetry:

# Skewness = E[(X - μ)³] / σ³
skewness = stats.skew(data, bias=False)
print(f"Skewness: {skewness:.3f}")

# Interpretation:
# 0: Symmetric (normal)
# >0: Right-skewed (positive, long right tail)
# <0: Left-skewed (negative, long left tail)

2. Positive Skew¶

# Generate right-skewed data
right_skew = np.random.exponential(scale=1, size=1000)
print(f"Skewness: {stats.skew(right_skew):.3f}")  # ~2.0

# Visualization
plt.hist(right_skew, bins=50, density=True, alpha=0.7)
plt.axvline(np.mean(right_skew), color='r', label='Mean')
plt.axvline(np.median(right_skew), color='g', label='Median')
plt.legend()
plt.title('Positive Skew: Mean > Median')
plt.show()

# Mean pulled toward tail
print(f"Mean: {np.mean(right_skew):.3f}")
print(f"Median: {np.median(right_skew):.3f}")  # Mean > Median

3. Negative Skew¶

# Generate left-skewed data
left_skew = -np.random.exponential(scale=1, size=1000)
print(f"Skewness: {stats.skew(left_skew):.3f}")  # ~-2.0

# Mean < Median for left skew

4. Skewness Test¶

# Test if skewness significantly differs from 0
skewtest_stat, p_value = stats.skewtest(data)
print(f"Skew test: statistic={skewtest_stat:.3f}, p={p_value:.4f}")

# p < 0.05: Significant skewness (reject symmetry)

5. Pearson Skewness Coefficients¶

# Mode skewness
mean = np.mean(data)
mode = stats.mode(data, keepdims=True).mode[0]
std = np.std(data, ddof=1)
pearson_mode = (mean - mode) / std

# Median skewness (more robust)
median = np.median(data)
pearson_median = 3 * (mean - median) / std

print(f"Pearson mode skewness: {pearson_mode:.3f}")
print(f"Pearson median skewness: {pearson_median:.3f}")

6. Quantile Skewness¶

# Based on quartiles (very robust)
q1, q2, q3 = np.percentile(data, [25, 50, 75])
bowley_skew = (q3 + q1 - 2*q2) / (q3 - q1)
print(f"Bowley skewness: {bowley_skew:.3f}")

# Ranges from -1 to +1

7. Transformations¶

# Box-Cox transformation to reduce skewness
from scipy import stats

data_positive = np.abs(data) + 1  # Must be positive
transformed, lambda_param = stats.boxcox(data_positive)

print(f"Original skewness: {stats.skew(data_positive):.3f}")
print(f"Transformed skewness: {stats.skew(transformed):.3f}")
print(f"Lambda parameter: {lambda_param:.3f}")

Kurtosis¶

1. Definition¶

Measure of tail heaviness:

# Excess kurtosis = E[(X - μ)⁴] / σ⁴ - 3
kurtosis = stats.kurtosis(data, bias=False)
print(f"Excess kurtosis: {kurtosis:.3f}")

# Interpretation:
# 0: Normal (mesokurtic)
# >0: Heavy tails (leptokurtic)
# <0: Light tails (platykurtic)

2. Leptokurtic (Heavy Tails)¶

# t-distribution has heavy tails
heavy_tail = stats.t.rvs(df=5, size=1000)
print(f"Kurtosis: {stats.kurtosis(heavy_tail):.3f}")  # Positive

# More extreme values than normal

3. Platykurtic (Light Tails)¶

# Uniform distribution has light tails
light_tail = np.random.uniform(-1, 1, 1000)
print(f"Kurtosis: {stats.kurtosis(light_tail):.3f}")  # ~-1.2

# Fewer extreme values than normal

4. Kurtosis Test¶

# Test if kurtosis significantly differs from 0 (normal)
kurttest_stat, p_value = stats.kurtosistest(data)
print(f"Kurt test: statistic={kurttest_stat:.3f}, p={p_value:.4f}")

# p < 0.05: Significant deviation from normal kurtosis

5. Sample Kurtosis Bias¶

# Unbiased vs biased estimator
kurt_unbiased = stats.kurtosis(data, bias=False)  # Recommended
kurt_biased = stats.kurtosis(data, bias=True)

print(f"Unbiased: {kurt_unbiased:.3f}")
print(f"Biased: {kurt_biased:.3f}")

6. Outlier Sensitivity¶

# Kurtosis very sensitive to outliers
data_clean = np.random.normal(0, 1, 100)
data_outlier = np.concatenate([data_clean, [10, -10]])

print(f"Clean kurtosis: {stats.kurtosis(data_clean):.3f}")
print(f"With outliers: {stats.kurtosis(data_outlier):.3f}")  # Much higher

7. Robust Alternatives¶

# L-kurtosis (based on L-moments, more robust)
# Not in scipy, but conceptually important

# Or use quantile-based measure
q = np.percentile(data, [10, 25, 75, 90])
robust_kurt = (q[3] - q[0]) / (q[2] - q[1])
print(f"Robust kurtosis proxy: {robust_kurt:.3f}")

Normal Test¶

1. Jarque-Bera Test¶

# Tests both skewness and kurtosis
jb_stat, p_value = stats.jarque_bera(data)
print(f"Jarque-Bera: statistic={jb_stat:.3f}, p={p_value:.4f}")

# H0: Data is normally distributed
# p < 0.05: Reject normality

2. Anderson-Darling Test¶

# More powerful than Jarque-Bera
result = stats.anderson(data, dist='norm')
print(f"A-D statistic: {result.statistic:.3f}")

# Compare to critical values
for i, sl in enumerate(result.significance_level):
    print(f"{sl}%: {result.critical_values[i]:.3f}")

3. Shapiro-Wilk Test¶

# Most powerful test for normality (small samples)
sw_stat, p_value = stats.shapiro(data)
print(f"Shapiro-Wilk: statistic={sw_stat:.4f}, p={p_value:.4f}")

# W close to 1 indicates normality

4. Kolmogorov-Smirnov Test¶

# General distribution test
ks_stat, p_value = stats.kstest(data, 'norm', args=(np.mean(data), np.std(data)))
print(f"K-S test: statistic={ks_stat:.4f}, p={p_value:.4f}")

5. D'Agostino-Pearson Test¶

# Combines skew and kurtosis tests
k2, p_value = stats.normaltest(data)
print(f"D'Agostino: statistic={k2:.3f}, p={p_value:.4f}")

# Equivalent to combining skewtest and kurtosistest

6. Q-Q Plot¶

# Visual normality check
stats.probplot(data, dist="norm", plot=plt)
plt.title("Q-Q Plot")
plt.show()

# Points on line → Normal
# S-curve → Heavy tails
# Reverse S → Light tails

7. Omnibus Test¶

# Not directly in scipy, but can be constructed
skew = stats.skew(data)
kurt = stats.kurtosis(data)
n = len(data)

# Jarque-Bera statistic
jb = n/6 * (skew**2 + kurt**2/4)
p_value = 1 - stats.chi2.cdf(jb, df=2)
print(f"Omnibus: JB={jb:.3f}, p={p_value:.4f}")

Practical Applications¶

1. Financial Returns¶

# Stock returns often have fat tails
returns = np.random.standard_t(df=5, size=1000) * 0.02  # 2% daily vol

print(f"Skewness: {stats.skew(returns):.3f}")
print(f"Kurtosis: {stats.kurtosis(returns):.3f}")  # Positive → Fat tails

# Implications: Normal models underestimate risk

2. Income Distribution¶

# Income typically right-skewed
income = np.random.lognormal(mean=10.5, sigma=0.5, size=1000)

print(f"Skewness: {stats.skew(income):.3f}")  # Strongly positive
print(f"Mean: ${np.mean(income):.0f}")
print(f"Median: ${np.median(income):.0f}")  # Median < Mean

3. Outlier Detection¶

# High kurtosis suggests outliers
if stats.kurtosis(data) > 3:
    print("Warning: Heavy tails detected")
    # Consider robust statistics or outlier removal

4. Transformation Selection¶

# Choose transformation based on skewness
if stats.skew(data) > 1:
    # Strong right skew → log transform
    transformed = np.log(data + 1)
elif stats.skew(data) < -1:
    # Strong left skew → exponential transform
    transformed = np.exp(data)
else:
    # Mild skew → square root or Box-Cox
    transformed, _ = stats.boxcox(data + abs(data.min()) + 1)

5. Model Selection¶

# Normal vs t-distribution based on kurtosis
if abs(stats.kurtosis(data)) < 1:
    model = stats.norm
    print("Using normal distribution")
else:
    # Fit t-distribution for heavy tails
    df, loc, scale = stats.t.fit(data)
    model = stats.t(df=df, loc=loc, scale=scale)
    print(f"Using t-distribution (df={df:.1f})")

6. Risk Measures¶

# VaR underestimated if kurtosis ignored
mean, std = np.mean(returns), np.std(returns)

# Normal assumption
var_normal = stats.norm.ppf(0.05, loc=mean, scale=std)

# Empirical (accounts for actual distribution)
var_empirical = np.percentile(returns, 5)

print(f"VaR (normal): {var_normal:.4f}")
print(f"VaR (empirical): {var_empirical:.4f}")

7. Power Analysis¶

# Sample size for detecting skewness/kurtosis
from statsmodels.stats.power import normal_power_het

# Example: Detect skewness of 0.5 with power 0.8
# Requires large sample size

Summary¶

Moments: - 1st: Mean (location) - 2nd: Variance (spread) - 3rd: Skewness (asymmetry) - 4th: Kurtosis (tail heaviness)

Skewness: - 0: Symmetric - >0: Right-skewed (mean > median) - <0: Left-skewed (mean < median)

Kurtosis (excess): - 0: Normal tails - >0: Heavy tails (more extremes) - <0: Light tails (fewer extremes)

Key insight: Higher moments reveal distribution shape beyond mean and variance, with skewness indicating asymmetry (important for risk assessment) and kurtosis measuring tail thickness (critical for extreme event modeling).