PMM3 for Time Series: AR, MA, ARMA, and ARIMA Models

Introduction

Classical time series estimation methods (CSS, MLE) assume Gaussian innovations. The PMM2 method addresses asymmetric non-Gaussian innovations, but many real-world time series exhibit symmetric platykurtic errors — lighter tails than the normal distribution. Examples include:

Bounded measurement errors (sensor readings within a fixed range)
Uniformly distributed shocks (e.g., rounding errors)
Truncated normal processes (censored observations)
Beta-symmetric disturbances (bounded economic indicators)

The PMM3 method (S=3) extends the Polynomial Maximization Method to these cases, using a cubic stochastic polynomial and Newton-Raphson solver based on the score equation \(Z_{1\nu} = \varepsilon_\nu(\kappa - \varepsilon_\nu^2)\).

This vignette demonstrates PMM3 for time series models:

AR — Autoregressive models
MA — Moving Average models
ARMA — Combined AR + MA models
ARIMA — Models with differencing for non-stationary series
Real-data examples using datasets included in the package
Comparison with MLE and PMM2

When PMM3 vs PMM2 vs OLS?

Innovation property	Recommended method
Skewed (\(\|\gamma_3\| > 0.3\))	PMM2 (`ar_pmm2`, `ma_pmm2`, …)
Symmetric, platykurtic (\(\gamma_4 < -0.7\))	PMM3 (`ar_pmm3`, `ma_pmm3`, …)
Near-Gaussian	MLE / CSS

Use pmm_dispatch() on residuals from a classical fit to automate this choice.

Setup

library(EstemPMM)
set.seed(42)

Part 1: AR Models with PMM3

AR(1) with Uniform Innovations

Uniform innovations are the canonical platykurtic case: \(\gamma_4 \approx -1.2\), symmetric, bounded.

n <- 300
phi_true <- 0.7

# Uniform innovations: symmetric, platykurtic (gamma4 ≈ -1.2)
innov <- runif(n, -sqrt(3), sqrt(3))  # variance = 1

# Generate AR(1) series
ar_data <- arima.sim(n = n, list(ar = phi_true), innov = innov)

# Plot
plot(ar_data, main = "AR(1) Process with Uniform Innovations",
     ylab = "Value", col = "steelblue", lwd = 1.2)

Fit AR(1): PMM3 vs MLE

# PMM3 fit
fit_pmm3_ar1 <- ar_pmm3(ar_data, order = 1, include.mean = FALSE)

# MLE fit (classical)
fit_mle_ar1 <- arima(ar_data, order = c(1, 0, 0), include.mean = FALSE)

# Compare
cat("True AR(1) coefficient:", phi_true, "\n")
#> True AR(1) coefficient: 0.7
cat("MLE estimate:          ", coef(fit_mle_ar1)["ar1"], "\n")
#> MLE estimate:           0.7431419
cat("PMM3 estimate:         ", coef(fit_pmm3_ar1)["ar1"], "\n")
#> PMM3 estimate:          0.7151052

# PMM3 diagnostics
summary(fit_pmm3_ar1)
#> PMM3 time series estimation results (S = 3, symmetric errors)
#> Model type: AR(1)
#> Call:
#> ts_pmm3(x = x, order = order, model_type = "ar", max_iter = max_iter, 
#>     tol = tol, adaptive = adaptive, step_max = step_max, include.mean = include.mean, 
#>     initial = initial, na.action = na.action, verbose = verbose)
#> 
#> Coefficients:
#> AR:       ar1 
#> 0.7151052 
#> 
#> Central moments of initial residuals:
#>   m2 = 1.013623 
#>   m4 = 1.837523 
#>   m6 = 3.978298 
#> 
#> Theoretical characteristics of PMM3 (S = 3):
#>   gamma4 = -1.211538 
#>   gamma6 = 6.993114 
#>   g3     = 0.297445  (expected ratio Var[PMM3]/Var[OLS])
#>   kappa  = 1.292902 
#> 
#> Algorithm information:
#>   Convergence status: TRUE 
#>   Iterations: 3

The g3 value in the summary shows the theoretical variance reduction factor. For uniform innovations, \(g_3 \approx 0.64\), meaning PMM3 achieves about 36% variance reduction over OLS.

AR(2) Model

phi_true <- c(0.5, -0.3)

# Beta-symmetric innovations: symmetric, platykurtic
innov2 <- (rbeta(400, 2, 2) - 0.5) * sqrt(12 * 2)  # scaled for var ≈ 1
ar2_data <- arima.sim(n = 400, list(ar = phi_true), innov = innov2)

fit_pmm3_ar2 <- ar_pmm3(ar2_data, order = 2, include.mean = FALSE)
fit_mle_ar2  <- arima(ar2_data, order = c(2, 0, 0), include.mean = FALSE)

comparison_ar2 <- data.frame(
  Parameter = c("ar1", "ar2"),
  True   = phi_true,
  MLE    = c(coef(fit_mle_ar2)["ar1"], coef(fit_mle_ar2)["ar2"]),
  PMM3   = coef(fit_pmm3_ar2)[c("ar1", "ar2")]
)

print(comparison_ar2, row.names = FALSE)
#>  Parameter True        MLE       PMM3
#>        ar1  0.5  0.4502433  0.4468475
#>        ar2 -0.3 -0.3468338 -0.3250832

cat("\nPMM3 g3 =", fit_pmm3_ar2@g_coefficient,
    "(expected variance ratio PMM3/OLS)\n")
#> 
#> PMM3 g3 = 0.6997661 (expected variance ratio PMM3/OLS)

Part 2: MA Models with PMM3

MA models are more challenging because innovations are not directly observable. PMM3 uses CSS residuals as initial estimates and refines them via Newton-Raphson.

MA(1) Model

n <- 300
theta_true <- 0.6

# Truncated normal innovations (|x| ≤ 2): symmetric, platykurtic
innov_raw <- rnorm(n * 2)
innov_tn <- innov_raw[abs(innov_raw) <= 2]
innov_tn <- innov_tn[1:n]
innov_tn <- innov_tn / sd(innov_tn)  # standardize

# Generate MA(1) series
ma_data <- arima.sim(n = n, list(ma = theta_true), innov = innov_tn)

# Fit
fit_pmm3_ma1 <- ma_pmm3(ma_data, order = 1, include.mean = FALSE)
fit_mle_ma1  <- arima(ma_data, order = c(0, 0, 1), include.mean = FALSE)

cat("True MA(1) coefficient:", theta_true, "\n")
#> True MA(1) coefficient: 0.6
cat("MLE estimate:          ", coef(fit_mle_ma1)["ma1"], "\n")
#> MLE estimate:           0.6569534
cat("PMM3 estimate:         ", coef(fit_pmm3_ma1)["ma1"], "\n")
#> PMM3 estimate:          0.6323384
cat("\nConvergence:", fit_pmm3_ma1@convergence, "\n")
#> 
#> Convergence: TRUE
cat("Iterations: ", fit_pmm3_ma1@iterations, "\n")
#> Iterations:  3

MA(2) Model

theta_true <- c(0.5, 0.3)

innov_u <- runif(400, -sqrt(3), sqrt(3))
ma2_data <- arima.sim(n = 400, list(ma = theta_true), innov = innov_u)

fit_pmm3_ma2 <- ma_pmm3(ma2_data, order = 2, include.mean = FALSE)
fit_mle_ma2  <- arima(ma2_data, order = c(0, 0, 2), include.mean = FALSE)

comparison_ma2 <- data.frame(
  Parameter = c("ma1", "ma2"),
  True   = theta_true,
  MLE    = c(coef(fit_mle_ma2)["ma1"], coef(fit_mle_ma2)["ma2"]),
  PMM3   = coef(fit_pmm3_ma2)[c("ma1", "ma2")]
)

print(comparison_ma2, row.names = FALSE)
#>  Parameter True       MLE      PMM3
#>        ma1  0.5 0.5262435 0.5171439
#>        ma2  0.3 0.3078768 0.2931746

Part 3: ARMA Models

ARMA models combine autoregressive and moving average components. PMM3 handles both components simultaneously.

ARMA(1,1) Model

n <- 400
phi_true   <- 0.5
theta_true <- 0.3

# Uniform innovations
innov_arma <- runif(n, -sqrt(3), sqrt(3))

arma_data <- arima.sim(n = n,
                       list(ar = phi_true, ma = theta_true),
                       innov = innov_arma)

# Fit
fit_pmm3_arma <- arma_pmm3(arma_data, order = c(1, 1), include.mean = FALSE)
fit_mle_arma  <- arima(arma_data, order = c(1, 0, 1), include.mean = FALSE)

cat("True: AR =", phi_true, ", MA =", theta_true, "\n\n")
#> True: AR = 0.5 , MA = 0.3

cat("MLE estimates:\n")
#> MLE estimates:
print(coef(fit_mle_arma)[c("ar1", "ma1")])
#>       ar1       ma1 
#> 0.4258325 0.3115166

cat("\nPMM3 estimates:\n")
#> 
#> PMM3 estimates:
print(coef(fit_pmm3_arma))
#>       ar1       ma1 
#> 0.4254178 0.3708970

summary(fit_pmm3_arma)
#> PMM3 time series estimation results (S = 3, symmetric errors)
#> Model type: ARMA(1,1)
#> Call:
#> ts_pmm3(x = x, order = order, model_type = "arma", max_iter = max_iter, 
#>     tol = tol, adaptive = adaptive, step_max = step_max, include.mean = include.mean, 
#>     initial = initial, na.action = na.action, verbose = verbose)
#> 
#> Coefficients:
#> AR:       ar1 
#> 0.4254178 
#> MA:      ma1 
#> 0.370897 
#> 
#> Central moments of initial residuals:
#>   m2 = 1.007571 
#>   m4 = 1.773758 
#>   m6 = 3.745631 
#> 
#> Theoretical characteristics of PMM3 (S = 3):
#>   gamma4 = -1.252799 
#>   gamma6 = 7.453816 
#>   g3     = 0.2795877  (expected ratio Var[PMM3]/Var[OLS])
#>   kappa  = 1.270544 
#> 
#> Algorithm information:
#>   Convergence status: TRUE 
#>   Iterations: 3

Diagnostic Plots

plot(fit_pmm3_arma)

The four diagnostic panels:

Residuals vs Time — should show no patterns
Q-Q Plot — deviations from the line indicate non-normality (expected for platykurtic errors: lighter tails than normal)
ACF — should fall within confidence bands (white noise)
Histogram — platykurtic residuals appear flatter than a bell curve

Part 4: ARIMA Models with Differencing

For non-stationary series, ARIMA models apply differencing before estimation. PMM3 uses a nonlinear solver with stats::arima(fixed=...) for proper residual computation and numDeriv::jacobian for derivatives. This approach correctly captures the full ARIMA dynamics after differencing.

ARIMA(1,1,0) Model

n <- 300
phi_true <- 0.6

# Generate stationary AR(1) with uniform innovations
innov_arima <- runif(n, -sqrt(3), sqrt(3))
x_stat <- arima.sim(n = n, list(ar = phi_true), innov = innov_arima)

# Integrate (cumulative sum) to create non-stationary series
x_integrated <- cumsum(x_stat)

plot(x_integrated, type = "l",
     main = "ARIMA(1,1,0) Process (Integrated AR(1))",
     ylab = "Value", xlab = "Time", col = "darkgreen", lwd = 1.2)

Fit ARIMA(1,1,0)

fit_pmm3_arima <- arima_pmm3(x_integrated, order = c(1, 1, 0),
                             include.mean = FALSE)
fit_mle_arima  <- arima(x_integrated, order = c(1, 1, 0),
                        include.mean = FALSE)

cat("True AR(1) coefficient:", phi_true, "\n")
#> True AR(1) coefficient: 0.6
cat("MLE estimate:          ", coef(fit_mle_arima)["ar1"], "\n")
#> MLE estimate:           0.5788831
cat("PMM3 estimate:         ", coef(fit_pmm3_arima)["ar1"], "\n")
#> PMM3 estimate:          0.6318376

summary(fit_pmm3_arima)
#> PMM3 time series estimation results (S = 3, symmetric errors)
#> Model type: ARIMA(1,1,0)
#> Call:
#> ts_pmm3(x = x, order = order, model_type = "arima", max_iter = max_iter, 
#>     tol = tol, adaptive = adaptive, step_max = step_max, include.mean = include.mean, 
#>     initial = initial, na.action = na.action, verbose = verbose)
#> 
#> Coefficients:
#> AR:       ar1 
#> 0.6318376 
#> 
#> Central moments of initial residuals:
#>   m2 = 1.034877 
#>   m4 = 1.910823 
#>   m6 = 4.259607 
#> 
#> Theoretical characteristics of PMM3 (S = 3):
#>   gamma4 = -1.215801 
#>   gamma6 = 7.080315 
#>   g3     = 0.3086522  (expected ratio Var[PMM3]/Var[OLS])
#>   kappa  = 1.2847 
#> 
#> Algorithm information:
#>   Convergence status: TRUE 
#>   Iterations: 5

ARIMA(1,1,1) Model

phi_true   <- 0.4
theta_true <- 0.3

innov_111 <- runif(350, -sqrt(3), sqrt(3))
x_stat2 <- arima.sim(n = 350, list(ar = phi_true, ma = theta_true),
                     innov = innov_111)
x_int2 <- cumsum(x_stat2)

fit_pmm3_111 <- arima_pmm3(x_int2, order = c(1, 1, 1), include.mean = FALSE)
fit_mle_111  <- arima(x_int2, order = c(1, 1, 1), include.mean = FALSE)

comparison_arima <- data.frame(
  Parameter = c("ar1", "ma1"),
  True   = c(phi_true, theta_true),
  MLE    = c(coef(fit_mle_111)["ar1"], coef(fit_mle_111)["ma1"]),
  PMM3   = coef(fit_pmm3_111)[c("ar1", "ma1")]
)

print(comparison_arima, row.names = FALSE)
#>  Parameter True       MLE      PMM3
#>        ar1  0.4 0.3234249 0.3745386
#>        ma1  0.3 0.2558448 0.3151259

Part 5: Forecasting with PMM3 Models

All PMM3 time series objects support the predict() method.

# Forecast 10 steps ahead from the AR(1) model
preds_ar <- predict(fit_pmm3_ar1, n.ahead = 10)

cat("AR(1) PMM3 forecasts (10 steps):\n")
#> AR(1) PMM3 forecasts (10 steps):
print(round(as.numeric(preds_ar), 4))
#>  [1] -1.8096 -1.2941 -0.9254 -0.6618 -0.4732 -0.3384 -0.2420 -0.1731 -0.1238
#> [10] -0.0885

# Plot original series with forecasts
n_plot <- 80
plot_data <- tail(as.numeric(ar_data), n_plot)
pred_vals <- as.numeric(preds_ar)

plot(seq_along(plot_data), plot_data, type = "l",
     xlim = c(1, n_plot + 10),
     ylim = range(c(plot_data, pred_vals)),
     main = "AR(1) PMM3: Observed and Forecast",
     xlab = "Time", ylab = "Value",
     col = "steelblue", lwd = 1.5)
lines(n_plot + 1:10, pred_vals, col = "red", lwd = 2, lty = 2)
legend("topleft", legend = c("Observed", "Forecast"),
       col = c("steelblue", "red"), lty = c(1, 2), lwd = c(1.5, 2))

Forecasting ARMA Models

# Forecast from the ARMA(1,1) model
preds_arma <- predict(fit_pmm3_arma, n.ahead = 10)

cat("ARMA(1,1) PMM3 forecasts (10 steps):\n")
#> ARMA(1,1) PMM3 forecasts (10 steps):
print(round(as.numeric(preds_arma), 4))
#>  [1] 1.4471 0.6156 0.2619 0.1114 0.0474 0.0202 0.0086 0.0036 0.0016 0.0007

Part 6: Method Selection with pmm_dispatch()

The pmm_dispatch() function automates the choice between OLS, PMM2, and PMM3 by analyzing innovation properties from an initial classical fit.

# Fit classical model first
fit_classical <- arima(ar_data, order = c(1, 0, 0), include.mean = FALSE)
res_classical <- residuals(fit_classical)

# Dispatch
recommendation <- pmm_dispatch(res_classical)
#>   n = 300 | gamma3 = +0.014 | gamma4 = -1.216
#>   g2(PMM2) = 0.9998 | g3(PMM3) = 0.2948
#>   >>> |gamma3| = 0.014 < 0.3 (symmetric) and gamma4 = -1.216 < -0.7 (platykurtic). PMM3 expected to reduce variance by 70.5%.

# Compare all three methods on uniform innovations
cat("Recommended method:", recommendation$method, "\n")
#> Recommended method: PMM3
cat("gamma3 =", round(recommendation$gamma3, 3),
    " gamma4 =", round(recommendation$gamma4, 3), "\n")
#> gamma3 = 0.014  gamma4 = -1.216

if (!is.null(recommendation$g3)) {
  cat("g3 =", round(recommendation$g3, 3),
      " (PMM3 variance reduction factor)\n")
}
#> g3 = 0.295  (PMM3 variance reduction factor)

Part 7: Real-Data Example — DJIA 2002

The djia2002 dataset contains daily closing prices and changes of the Dow Jones Industrial Average for July–December 2002.

data(djia2002)
changes <- na.omit(djia2002$change)
n_obs <- length(changes)

plot(changes, type = "l",
     main = "DJIA Daily Changes (Jul-Dec 2002)",
     ylab = "Change (USD)", xlab = "Trading Day",
     col = "steelblue", lwd = 1.2)
abline(h = 0, col = "red", lty = 2)

Analyze Innovation Properties

# Fit a classical AR(1)
fit_djia_mle <- arima(changes, order = c(1, 0, 0))
res_djia <- residuals(fit_djia_mle)

cat("Innovation diagnostics:\n")
#> Innovation diagnostics:
cat("  Skewness (gamma3):", round(pmm_skewness(res_djia), 3), "\n")
#>   Skewness (gamma3): 0.445
cat("  Kurtosis (gamma4):", round(pmm_kurtosis(res_djia), 3), "\n")
#>   Kurtosis (gamma4): 0.263
cat("  Symmetry test:\n")
#>   Symmetry test:
sym <- test_symmetry(res_djia)
cat("    symmetric:", sym$is_symmetric, "\n")
#>     symmetric: FALSE
cat("    message:  ", sym$message, "\n")
#>     message:   |gamma3| = 0.445 > 0.3: residuals are asymmetric. Consider PMM2 instead.

# Let pmm_dispatch decide
dispatch <- pmm_dispatch(res_djia)
#>   n = 126 | gamma3 = +0.445 | gamma4 = +0.263
#>   g2(PMM2) = 0.9126 | g3(PMM3) = 0.9900
#>   >>> |gamma3| = 0.445 > 0.3 and g2 = 0.913 < 0.95: moderate asymmetry, PMM2 worthwhile (8.7% reduction).
cat("\nRecommended method:", dispatch$method, "\n")
#> 
#> Recommended method: PMM2

Fit AR(1) with All Three Methods

fit_djia_pmm2 <- ar_pmm2(changes, order = 1)
fit_djia_pmm3 <- ar_pmm3(changes, order = 1)

comparison_djia <- data.frame(
  Method = c("MLE", "PMM2", "PMM3"),
  ar1 = c(coef(fit_djia_mle)["ar1"],
          coef(fit_djia_pmm2)["ar1"],
          coef(fit_djia_pmm3)["ar1"])
)
print(comparison_djia, row.names = FALSE)
#>  Method         ar1
#>     MLE -0.03956444
#>    PMM2 -0.02144724
#>    PMM3 -0.04472054

The DJIA data is known to have positive skewness, so pmm_dispatch() typically recommends PMM2 as the appropriate method for this series. PMM3 may still converge but is not expected to improve upon PMM2 when the innovations are asymmetric.

Part 8: Real-Data Example — WTI Crude Oil Prices

The DCOILWTICO dataset contains WTI crude oil daily spot prices.

data(DCOILWTICO)

# Compute log-returns (remove NAs)
prices <- na.omit(DCOILWTICO$DCOILWTICO)
log_returns <- diff(log(prices))
log_returns <- log_returns[is.finite(log_returns)]

plot(log_returns, type = "l",
     main = "WTI Crude Oil: Daily Log-Returns",
     ylab = "Log-return", xlab = "Trading Day",
     col = "darkgreen", lwd = 0.8)
abline(h = 0, col = "red", lty = 2)

Innovation Analysis

# Fit AR(1) to log-returns
fit_oil_mle <- arima(log_returns, order = c(1, 0, 0))
res_oil <- residuals(fit_oil_mle)

cat("Log-return innovation diagnostics:\n")
#> Log-return innovation diagnostics:
cat("  Skewness (gamma3):", round(pmm_skewness(res_oil), 3), "\n")
#>   Skewness (gamma3): -0.498
cat("  Kurtosis (gamma4):", round(pmm_kurtosis(res_oil), 3), "\n")
#>   Kurtosis (gamma4): 2.093

# Dispatch
dispatch_oil <- pmm_dispatch(res_oil)
#>   n = 1248 | gamma3 = -0.498 | gamma4 = +2.093
#>   g2(PMM2) = 0.9394 | g3(PMM3) = 0.9010
#>   >>> |gamma3| = 0.498 > 0.3 and g2 = 0.939 < 0.95: moderate asymmetry, PMM2 worthwhile (6.1% reduction).
cat("\nRecommended method:", dispatch_oil$method, "\n")
#> 
#> Recommended method: PMM2

Fit and Compare

fit_oil_pmm2 <- ar_pmm2(log_returns, order = 1)
fit_oil_pmm3 <- ar_pmm3(log_returns, order = 1)

comparison_oil <- data.frame(
  Method = c("MLE", "PMM2", "PMM3"),
  ar1 = c(coef(fit_oil_mle)["ar1"],
          coef(fit_oil_pmm2)["ar1"],
          coef(fit_oil_pmm3)["ar1"])
)
print(comparison_oil, row.names = FALSE)
#>  Method        ar1
#>     MLE 0.02971020
#>    PMM2 0.02250418
#>    PMM3 0.04288416

vf2_oil <- pmm2_variance_factor(fit_oil_pmm2@m2, fit_oil_pmm2@m3, fit_oil_pmm2@m4)
cat("\nPMM2 g2 =", vf2_oil$g2, "\n")
#> 
#> PMM2 g2 =
cat("PMM3 g3 =", fit_oil_pmm3@g_coefficient, "\n")
#> PMM3 g3 = 0.9002976

Crude oil log-returns typically exhibit heavy tails (positive excess kurtosis), so PMM2 is usually the better choice. This example demonstrates how pmm_dispatch() guides users to the correct method.

Part 9: Adaptive Mode

By default, PMM3 computes \(\kappa\) from the initial residuals and holds it fixed. In adaptive mode, \(\kappa\) is re-estimated from the current residuals at each Newton-Raphson iteration:

# Fixed kappa (default)
fit_fixed <- ar_pmm3(ar_data, order = 1, include.mean = FALSE, adaptive = FALSE)

# Adaptive kappa
fit_adapt <- ar_pmm3(ar_data, order = 1, include.mean = FALSE, adaptive = TRUE)

comparison_adaptive <- data.frame(
  Mode   = c("Fixed kappa", "Adaptive kappa"),
  ar1    = c(coef(fit_fixed)["ar1"], coef(fit_adapt)["ar1"]),
  g3     = c(fit_fixed@g_coefficient, fit_adapt@g_coefficient),
  iters  = c(fit_fixed@iterations, fit_adapt@iterations)
)
print(comparison_adaptive, row.names = FALSE)
#>            Mode       ar1       g3 iters
#>     Fixed kappa 0.7151052 0.297445     3
#>  Adaptive kappa 0.7147911 0.297445     4

Adaptive mode may improve estimates when the distribution departs significantly from the initial OLS residuals, at the cost of additional iterations.

Part 10: Model Comparison via AIC

The AIC() method is available for all PMM3 time series objects.

# Fit several AR orders with PMM3
aic_vals <- sapply(1:5, function(p) {
  fit <- ar_pmm3(ar_data, order = p, include.mean = FALSE)
  AIC(fit)
})

aic_df <- data.frame(Order = 1:5, AIC = round(aic_vals, 2))
print(aic_df, row.names = FALSE)
#>  Order    AIC
#>      1 858.54
#>      2 861.11
#>      3 863.25
#>      4 865.24
#>      5 867.37

cat("\nBest AR order by AIC:", which.min(aic_vals), "\n")
#> 
#> Best AR order by AIC: 1

Practical Guidelines

When to Use PMM3 for Time Series

Recommended when innovations are:

Symmetric (\(|\gamma_3| \leq 0.3\))
Platykurtic (\(\gamma_4 < -0.7\))
Example distributions: uniform, beta-symmetric, truncated normal
Sample size \(n \geq 200\) for reliable moment estimation

Use PMM2 instead when:

Innovations are skewed (\(|\gamma_3| > 0.3\))
Financial returns with asymmetric tails

Stick with MLE/CSS when:

Innovations are approximately Gaussian
Very small samples (\(n < 100\))

Diagnostic Workflow

Fit a classical model (MLE/CSS) as baseline
Examine residuals for non-normality:
- pmm_skewness(residuals) — check symmetry
- pmm_kurtosis(residuals) — check tail behavior
- test_symmetry(residuals) — formal symmetry check
Use pmm_dispatch(residuals) for automatic method selection
Refit with PMM3 if platykurtic and symmetric
Compare AIC across model orders
Validate with diagnostic plots and out-of-sample forecasts

Computational Notes

AR models: Fast convergence (3–10 iterations), linearized solver
MA models: Moderate (10–30 iterations), linearized solver
ARMA models: Moderate (10–30 iterations), linearized solver
ARIMA models: Nonlinear solver with numerical Jacobian — slower but accurately captures differencing dynamics

Efficiency Summary

Based on Monte Carlo simulations (R = 1000, n = 200, uniform innovations):

Model	MSE Ratio (PMM3/MLE)	Variance Reduction
AR(1)	0.55–0.70	30–45%
AR(2)	0.60–0.75	25–40%
MA(1)	0.90–1.00	0–10%
ARMA(1,1)	0.70–0.85	15–30%
ARIMA(1,1,0)	0.32–0.70	30–68%

Key insight: PMM3 provides the largest gains for AR and ARIMA models with platykurtic innovations. MA models show smaller improvements because CSS residuals are less precise starting points, but PMM3 remains consistent.

Available Functions

Function	Model	Class Returned
`ar_pmm3()`	AR(p)	`ARPMM3`
`ma_pmm3()`	MA(q)	`MAPMM3`
`arma_pmm3()`	ARMA(p,q)	`ARMAPMM3`
`arima_pmm3()`	ARIMA(p,d,q)	`ARIMAPMM3`
`ts_pmm3()`	Any of above	`TS3fit` subclass

All classes support: coef(), residuals(), fitted(), predict(), summary(), plot(), AIC().

Next Steps

PMM3 for linear regression: See vignette("pmm3_symmetric_errors")
PMM2 time series: See vignette("pmm2_time_series") for asymmetric innovations
Bootstrap inference: See vignette("bootstrap_inference")
PMM2 basics: See vignette("pmm2_introduction")

References

Zabolotnii S., Warsza Z.L., Tkachenko O. (2018). Polynomial Estimation of Linear Regression Parameters for the Asymmetric PDF of Errors. Springer AISC, vol 743. doi:10.1007/978-3-319-77179-3_75
Zabolotnii S., Tkachenko O., Warsza Z.L. (2022). Application of the Polynomial Maximization Method for Estimation Parameters of Autoregressive Models with Asymmetric Innovations. Springer AISC, vol 1427. doi:10.1007/978-3-031-03502-9_37
Package functions: ?ar_pmm3, ?ma_pmm3, ?arma_pmm3, ?arima_pmm3
Method selection: ?pmm_dispatch, ?test_symmetry
GitHub: https://github.com/SZabolotnii/EstemPMM