Introduction to PMM2: Polynomial Maximization Method

The Problem: When OLS Falls Short

Ordinary Least Squares (OLS) is the workhorse of regression analysis, providing Best Linear Unbiased Estimators (BLUE) under the Gauss-Markov assumptions. However, when error distributions deviate from normality—particularly when they exhibit skewness—OLS estimates, while still unbiased, lose statistical efficiency. This means:

• Higher variance in parameter estimates
• Wider confidence intervals
• Reduced statistical power for hypothesis tests

The Solution: Polynomial Maximization Method (PMM2)

The Polynomial Maximization Method of order 2 (PMM2) is a moment-based estimation technique that leverages higher-order moments of residuals to achieve:

1. Lower variance estimates compared to OLS when errors are non-Gaussian
2. Robustness to skewness and kurtosis in error distributions
3. Asymptotic efficiency gains validated through Monte Carlo studies

Getting Started with EstemPMM

Let’s explore PMM2 with practical examples using the EstemPMM package.

library(EstemPMM)
set.seed(123)

Example 1: Simple Linear Regression with Skewed Errors

We’ll simulate data where errors follow a skewed distribution (using \(\chi^2\)) to demonstrate PMM2’s advantages.

# Sample size
n <- 200

# Generate predictor
x <- rnorm(n, mean = 5, sd = 2)

# True parameters
beta_0 <- 10
beta_1 <- 2.5

# Generate highly skewed errors using chi-squared distribution
# (chi-square with df=3 has skewness ~ 1.63)
errors <- rchisq(n, df = 3) - 3  # Center to have mean ~ 0

# Response variable
y <- beta_0 + beta_1 * x + errors

# Create data frame
dat <- data.frame(y = y, x = x)

# Fit OLS to get residuals for visualization
fit_ols_temp <- lm(y ~ x, data = dat)
residuals_ols <- residuals(fit_ols_temp)

# Plot histogram of residuals
par(mfrow = c(1, 2))
hist(residuals_ols, breaks = 30, main = "Residual Distribution",
     xlab = "Residuals", col = "lightblue", border = "white")
qqnorm(residuals_ols, main = "Q-Q Plot")
qqline(residuals_ols, col = "red", lwd = 2)

par(mfrow = c(1, 1))

# Calculate skewness and kurtosis
cat(sprintf("Skewness: %.3f\n", pmm_skewness(residuals_ols)))
#> Skewness: 1.865
cat(sprintf("Excess Kurtosis: %.3f\n", pmm_kurtosis(residuals_ols)))
#> Excess Kurtosis: 4.710

# Fit PMM2 model
fit_pmm2 <- lm_pmm2(y ~ x, data = dat, verbose = FALSE)

# Fit OLS model for comparison
fit_ols <- lm(y ~ x, data = dat)

# Display coefficients
cat("=== OLS Estimates ===\n")
#> === OLS Estimates ===
print(coef(fit_ols))
#> (Intercept)           x 
#>   10.188412    2.423676

cat("\n=== PMM2 Estimates ===\n")
#> 
#> === PMM2 Estimates ===
print(coef(fit_pmm2))
#> (Intercept)           x 
#>    9.980738    2.465210

cat("\n=== True Values ===\n")
#> 
#> === True Values ===
cat("Intercept:", beta_0, "\n")
#> Intercept: 10
cat("Slope:    ", beta_1, "\n")
#> Slope:     2.5

# Perform bootstrap inference
boot_results <- pmm2_inference(fit_pmm2, formula = y ~ x, data = dat,
                               B = 500)

# Display summary with confidence intervals
summary(boot_results)
#>     Estimate       Std.Error            bias               t.value     
#>  Min.   :2.465   Min.   :0.06404   Min.   :-0.0085748   Min.   :27.25  
#>  1st Qu.:4.344   1st Qu.:0.13959   1st Qu.:-0.0060182   1st Qu.:30.06  
#>  Median :6.223   Median :0.21514   Median :-0.0034615   Median :32.87  
#>  Mean   :6.223   Mean   :0.21514   Mean   :-0.0034615   Mean   :32.87  
#>  3rd Qu.:8.102   3rd Qu.:0.29069   3rd Qu.:-0.0009049   3rd Qu.:35.68  
#>  Max.   :9.981   Max.   :0.36624   Max.   : 0.0016517   Max.   :38.49  
#>     p.value              conf.low       conf.high     
#>  Min.   : 0.000e+00   Min.   :2.340   Min.   : 2.591  
#>  1st Qu.:3.960e-164   1st Qu.:4.070   1st Qu.: 4.618  
#>  Median :7.921e-164   Median :5.801   Median : 6.645  
#>  Mean   :7.921e-164   Mean   :5.801   Mean   : 6.645  
#>  3rd Qu.:1.188e-163   3rd Qu.:7.532   3rd Qu.: 8.672  
#>  Max.   :1.584e-163   Max.   :9.263   Max.   :10.699

Example 2: Multiple Regression

PMM2 extends naturally to multiple predictors.

# Generate multiple predictors
n <- 250
x1 <- rnorm(n, mean = 0, sd = 1)
x2 <- runif(n, min = -2, max = 2)
x3 <- rexp(n, rate = 0.5)

# True model
beta <- c(5, 1.2, -0.8, 2.0)  # Intercept, x1, x2, x3
errors <- rgamma(n, shape = 2) - 2  # skewed gamma errors (gamma3 ~ 1.41)

y <- beta[1] + beta[2]*x1 + beta[3]*x2 + beta[4]*x3 + errors

# Create data frame
dat_multi <- data.frame(y = y, x1 = x1, x2 = x2, x3 = x3)

# Fit models
fit_pmm2_multi <- lm_pmm2(y ~ x1 + x2 + x3, data = dat_multi, verbose = FALSE)
fit_ols_multi <- lm(y ~ x1 + x2 + x3, data = dat_multi)

# Compare coefficients
comparison <- data.frame(
  Parameter = c("Intercept", "x1", "x2", "x3"),
  True = beta,
  OLS = coef(fit_ols_multi),
  PMM2 = coef(fit_pmm2_multi)
)

print(comparison, row.names = FALSE)
#>  Parameter True        OLS       PMM2
#>  Intercept  5.0  4.8525846  4.8556232
#>         x1  1.2  1.1399811  1.1067664
#>         x2 -0.8 -0.8521322 -0.8109516
#>         x3  2.0  2.0548042  2.0529738

Example 3: Polynomial Regression

PMM2 handles polynomial terms seamlessly.

# Generate data with quadratic relationship
n <- 200
x <- seq(-3, 3, length.out = n)
errors <- rchisq(n, df = 5) - 5

# True quadratic model: y = 2 + 3x - 0.5x^2
y <- 2 + 3*x - 0.5*x^2 + errors

dat_poly <- data.frame(y = y, x = x)

# Fit polynomial model
fit_pmm2_poly <- lm_pmm2(y ~ x + I(x^2), data = dat_poly, verbose = FALSE)

# Summary
summary(fit_pmm2_poly, formula = y ~ x + I(x^2), data = dat_poly)
#> PMM2 estimation results
#> Call:
#> lm_pmm2(formula = y ~ x + I(x^2), data = dat_poly, verbose = FALSE)
#> 
#> Coefficients (asymptotic standard errors, see ?vcov.PMM2fit):
#>              Estimate Std. Error z value  Pr(>|z|)    
#> (Intercept)  1.783387   0.296107  6.0228 1.715e-09 ***
#> x            3.014761   0.103356 29.1687 < 2.2e-16 ***
#> I(x^2)      -0.424796   0.066386 -6.3989 1.565e-10 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Central moments of initial residuals:
#>   m2 = 9.442913 
#>   m3 = 26.18281 
#>   m4 = 320.0677 
#> 
#> Theoretical characteristics of PMM2 (S = 2):
#>   c3 = 0.9023129 
#>   c4 = 0.5894661 
#>   g  = 0.6855844  (expected ratio Var[PMM2]/Var[OLS])
#> 
#> Algorithm information:
#>   Convergence status: TRUE 
#>   Iterations: 3 
#> 
#> Approximate statistical inference via bootstrap (B=100):
#>               Estimate  Std.Error         bias   t.value       p.value
#> (Intercept)  1.7833866 0.30186673  0.027088085  5.907861  3.465791e-09
#> x            3.0147610 0.09463962  0.011280959 31.855169 1.116315e-222
#> I(x^2)      -0.4247964 0.06292167 -0.003420012 -6.751194  1.466333e-11
#>               conf.low  conf.high
#> (Intercept)  1.1917386  2.3750345
#> x            2.8292707  3.2002512
#> I(x^2)      -0.5481206 -0.3014722

Diagnostic Plots

PMM2 provides standard diagnostic methods:

# Generate diagnostic plots
plot(fit_pmm2_multi)
title("PMM2 Diagnostics")

The diagnostic plots include:

1. Residuals vs Fitted: Check for non-linearity and heteroscedasticity
2. Q-Q Plot: Assess normality of residuals
3. Scale-Location: Check homoscedasticity assumption
4. Residuals vs Leverage: Identify influential observations

Comparing PMM2 with OLS: Monte Carlo Evidence

Based on simulation studies documented in the package:

These ratios less than 1.0 confirm that PMM2 achieves lower Mean Squared Error compared to classical methods (CSS/OLS) when errors deviate from normality.

When to Use PMM2

Use PMM2 when:

• Error distributions exhibit skewness (PMM2 draws its gain from \(\gamma_3\))
• You need more precise estimates than OLS provides
• Working with financial, environmental, or biomedical data where non-normality is common
• Sample sizes are moderate to large (n ≥ 100 recommended)

Do not use PMM2 when:

• Errors are symmetric, however heavy their tails: then \(g_2 = 1\) (see below). Use PMM3 (lm_pmm3()) or, for outliers, lm_wpmm3().

Stick with OLS when:

• Errors are approximately Gaussian
• Sample sizes are very small relative to the number of regressors: at n = 30 with five regressors, lm_pmm2() returns the OLS start in 14-28% of samples with skewed errors
• Computational simplicity is paramount
• You need unbiased estimates regardless of efficiency

Theoretical Limits

Symmetric errors: degree 2 gains nothing. With \(\gamma_3 = 0\) the factor \(g_2 = 1 - \gamma_3^2 / (2 + \gamma_4)\) equals 1 exactly, whatever \(\gamma_4\) is. The same holds for every estimator built on linear and quadratic functions of the residuals (second-order least squares, the linear-plus-quadratic optimal estimating function). On symmetric laws all of the reserve lies at degree 3: lm_pmm3(). If PMM2 beats OLS on symmetric errors in a small sample, the gain comes from noise in the sample skewness. It shrinks as \(n\) grows: in the calibration study of pmm_dispatch(), the PMM2/OLS MSE ratio for Laplace errors stays between 0.94 and 0.99, and it approaches 1 as \(n\) grows.

Symmetric heavy tails: expect a tie with robust estimators. On symmetric laws, robust estimators (Huber) and rank estimators (Wilcoxon / Hodges-Lehmann) are already close to efficient. PMM3 and the weak-moment estimators lm_wpmm3() can therefore match them but rarely beat them. PMM’s clear advantages are skewed errors and platykurtic or bimodal symmetric errors, where those estimators do not help.

PMM and GMM. PMM2 solves the Godambe-optimal estimating equation built from the moment conditions \(E[\varepsilon^k] = \alpha_k\), \(k \le S\). Its asymptotic variance equals that of GMM on the same conditions with the optimal weight matrix. The two are the same estimator only asymptotically. In a finite sample PMM takes a root of the estimating equation, GMM takes a minimum of a quadratic form, and the PMM root need not minimise the GMM criterion. The estimating equation can also have several roots. The consistent one is the root reached by Newton iteration from a \(\sqrt{n}\)-consistent start. lm_pmm2() starts from OLS, damps steps that do not reduce the estimating function, and warns when the Jacobian at the root is not positive definite.

Practical Recommendations

1. Always fit OLS first as a baseline comparison
2. Check residual diagnostics to assess non-normality
3. Use bootstrap inference (pmm2_inference()) for robust confidence intervals
4. Compare standard errors between OLS and PMM2 to quantify efficiency gains
5. Validate results with out-of-sample prediction when possible

Next Steps

To learn more about EstemPMM:

• Time Series Applications: See vignette("pmm2_time_series") for AR, MA, ARMA, and ARIMA models
• Bootstrap Inference: See vignette("bootstrap_inference") for detailed statistical inference procedures
• Symmetric Errors: See vignette("pmm3_symmetric_errors") for PMM3
• Theoretical Background: Refer to Zabolotnii et al. (2018) cited in package documentation

Summary

The PMM2 method provides a powerful alternative to OLS when error distributions deviate from normality. By leveraging higher-order moments, PMM2 achieves:

• 10-50% reduction in parameter estimate variance
• Improved statistical power for hypothesis testing
• Robust performance across various non-Gaussian error distributions

The EstemPMM package makes PMM2 accessible through an intuitive interface that mirrors standard R modeling functions like lm().

References

1. Zabolotnii, S., Warsza, Z.L., Tkachenko, O. (2018). “Polynomial Estimation of Linear Regression Parameters for the Asymmetric PDF of Errors”. In: Automation 2018. Advances in Intelligent Systems and Computing, vol 743. Springer, Cham. https://doi.org/10.1007/978-3-319-77179-3_75

2. Package documentation: ?lm_pmm2, ?pmm2_inference, ?compare_with_ols

3. GitHub repository: https://github.com/SZabolotnii/EstemPMM