--- title: "PMM2 for Time Series: AR, MA, ARMA, ARIMA, and Seasonal Models" author: "Serhii Zabolotnii" date: "`r Sys.Date()`" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{PMM2 for Time Series: AR, MA, ARMA, ARIMA, and Seasonal Models} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 5, warning = FALSE, message = FALSE ) ``` ## Introduction Time series models are fundamental tools in economics, finance, engineering, and many other fields. Classical estimation methods like **Conditional Sum of Squares (CSS)** and **Maximum Likelihood (ML)** assume Gaussian innovations. However, real-world time series often exhibit: - **Heavy tails** (e.g., financial returns) - **Skewness** (e.g., environmental data, commodity prices) - **Outliers** (e.g., measurement errors) The **PMM2 method** extends to time series models, providing **more efficient parameter estimates** when innovations deviate from normality. This vignette demonstrates PMM2 applications to: 1. **Autoregressive (AR)** models 2. **Moving Average (MA)** models 3. **ARMA** models 4. **ARIMA** models with differencing 5. **Seasonal models (SAR, SMA)** for data with seasonal patterns --- ## Setup ```{r load_package} library(EstemPMM) set.seed(42) ``` --- ## Part 1: Autoregressive (AR) Models ### AR(1) Model with Skewed Innovations We'll simulate an AR(1) process with **gamma-distributed innovations** to demonstrate PMM2's robustness. ```{r ar1_simulation} # Model parameters n <- 300 phi <- 0.6 # AR coefficient intercept <- 5 # Generate skewed innovations using gamma distribution # (Gamma(shape=2) has skewness = sqrt(2) ≈ 1.41) innovations <- rgamma(n, shape = 2, rate = 2) - 1 # Center around 0 # Generate AR(1) series x <- numeric(n) x[1] <- intercept + innovations[1] for(i in 2:n) { x[i] <- intercept + phi * (x[i-1] - intercept) + innovations[i] } # Plot the time series plot(x, type = "l", main = "AR(1) Process with Skewed Innovations", ylab = "Value", xlab = "Time", col = "steelblue", lwd = 1.5) abline(h = intercept, col = "red", lty = 2) ``` ### Estimate AR(1) Model: PMM2 vs CSS ```{r ar1_estimation} # Fit using PMM2 fit_pmm2_ar1 <- ar_pmm2(x, order = 1, method = "pmm2", include.mean = TRUE) # Fit using CSS (classical method) fit_css_ar1 <- ar_pmm2(x, order = 1, method = "css", include.mean = TRUE) # Compare coefficients cat("True AR coefficient:", phi, "\n\n") cat("CSS Estimate:\n") print(coef(fit_css_ar1)) cat("\nPMM2 Estimate:\n") print(coef(fit_pmm2_ar1)) # Summary with diagnostics summary(fit_pmm2_ar1) ``` **Observation:** PMM2 typically provides estimates closer to the true value when innovations are non-Gaussian. --- ### AR(2) Model For higher-order AR models, PMM2 maintains its efficiency advantages. ```{r ar2_example} # AR(2) model: x_t = 0.5*x_{t-1} - 0.3*x_{t-2} + e_t n <- 400 phi <- c(0.5, -0.3) # Generate t-distributed innovations (heavy tails) innovations <- rt(n, df = 4) # Generate AR(2) series x <- arima.sim(n = n, model = list(ar = phi), innov = innovations) # Fit models fit_pmm2_ar2 <- ar_pmm2(x, order = 2, method = "pmm2", include.mean = FALSE) fit_css_ar2 <- ar_pmm2(x, order = 2, method = "css", include.mean = FALSE) # Compare comparison_ar2 <- data.frame( Parameter = c("phi1", "phi2"), True = phi, CSS = coef(fit_css_ar2), PMM2 = coef(fit_pmm2_ar2) ) print(comparison_ar2, row.names = FALSE) ``` --- ## Part 2: Moving Average (MA) Models MA models are challenging because innovations are not directly observable. PMM2 uses **CSS residuals** as starting estimates, then refines parameters. ### MA(1) Model ```{r ma1_simulation} # MA(1) model: x_t = e_t + theta*e_{t-1} n <- 300 theta <- 0.7 # Generate chi-squared innovations (right-skewed) innovations <- rchisq(n, df = 3) - 3 # Generate MA(1) series x <- numeric(n) x[1] <- innovations[1] for(i in 2:n) { x[i] <- innovations[i] + theta * innovations[i-1] } # Fit MA(1) models fit_pmm2_ma1 <- ma_pmm2(x, order = 1, method = "pmm2", include.mean = FALSE) fit_css_ma1 <- ma_pmm2(x, order = 1, method = "css", include.mean = FALSE) # Compare estimates cat("True MA coefficient:", theta, "\n\n") cat("CSS Estimate:\n") print(coef(fit_css_ma1)) cat("\nPMM2 Estimate:\n") print(coef(fit_pmm2_ma1)) # Check convergence cat("\nPMM2 Convergence:", fit_pmm2_ma1@convergence, "\n") cat("Iterations:", fit_pmm2_ma1@iterations, "\n") ``` --- ### MA(2) Model ```{r ma2_example} # MA(2) model n <- 400 theta <- c(0.6, 0.3) # Mixed distribution: 80% normal + 20% contamination innovations <- ifelse(runif(n) < 0.8, rnorm(n), rnorm(n, mean = 0, sd = 3)) # Generate MA(2) series x <- arima.sim(n = n, model = list(ma = theta), innov = innovations) # Fit models fit_pmm2_ma2 <- ma_pmm2(x, order = 2, method = "pmm2", include.mean = FALSE) fit_css_ma2 <- ma_pmm2(x, order = 2, method = "css", include.mean = FALSE) # Comparison comparison_ma2 <- data.frame( Parameter = c("theta1", "theta2"), True = theta, CSS = coef(fit_css_ma2), PMM2 = coef(fit_pmm2_ma2) ) print(comparison_ma2, row.names = FALSE) ``` **Monte Carlo Evidence:** For MA(2) models with skewed errors, PMM2 achieves **26-32% variance reduction** compared to CSS (see roadmap documentation). --- ## Part 3: ARMA Models ARMA models combine AR and MA components, providing flexible modeling for stationary series. ### ARMA(1,1) Model ```{r arma11_simulation} # ARMA(1,1) model n <- 500 phi <- 0.5 # AR coefficient theta <- 0.4 # MA coefficient # Generate exponentially distributed innovations (right-skewed) innovations <- rexp(n, rate = 1) - 1 # Generate ARMA(1,1) series x <- arima.sim(n = n, model = list(ar = phi, ma = theta), innov = innovations) # Fit models fit_pmm2_arma <- arma_pmm2(x, order = c(1, 1), method = "pmm2", include.mean = FALSE) fit_css_arma <- arma_pmm2(x, order = c(1, 1), method = "css", include.mean = FALSE) # Extract coefficients cat("True parameters: AR =", phi, ", MA =", theta, "\n\n") cat("CSS Estimates:\n") print(coef(fit_css_arma)) cat("\nPMM2 Estimates:\n") print(coef(fit_pmm2_arma)) # Summary summary(fit_pmm2_arma) ``` --- ### Diagnostic Plots for ARMA Models ```{r arma_diagnostics, fig.width=7, fig.height=6} # Plot residuals and ACF plot(fit_pmm2_arma) title("ARMA(1,1) Model Diagnostics") ``` The diagnostic plots help verify: 1. **Residuals vs Time**: Should show no patterns 2. **ACF of Residuals**: Should be within confidence bands (white noise) 3. **Q-Q Plot**: Assess residual distribution --- ## Part 4: ARIMA Models with Differencing For non-stationary series, ARIMA models incorporate differencing to achieve stationarity. ### ARIMA(1,1,1) Model ```{r arima111_simulation} # Generate non-stationary series requiring differencing n <- 400 phi <- 0.4 theta <- 0.3 # Skewed innovations innovations <- rgamma(n, shape = 3, rate = 3) - 1 # Generate integrated series: ARIMA(1,1,1) # First generate stationary ARMA(1,1) x_stationary <- arima.sim(n = n, model = list(ar = phi, ma = theta), innov = innovations) # Integrate once (cumulative sum) x <- cumsum(x_stationary) # Plot the non-stationary series plot(x, type = "l", main = "ARIMA(1,1,1) Process", ylab = "Value", xlab = "Time", col = "darkgreen", lwd = 1.5) ``` ### Fitting ARIMA Models ```{r arima_estimation} # Fit ARIMA(1,1,1) using PMM2 fit_pmm2_arima <- arima_pmm2(x, order = c(1, 1, 1), method = "pmm2", include.mean = FALSE) # Fit using classical CSS method fit_css_arima <- arima_pmm2(x, order = c(1, 1, 1), method = "css", include.mean = FALSE) # Compare coefficients cat("True parameters: AR =", phi, ", MA =", theta, "\n\n") cat("CSS Estimates:\n") print(coef(fit_css_arima)) cat("\nPMM2 Estimates:\n") print(coef(fit_pmm2_arima)) # Summary summary(fit_pmm2_arima) ``` **Monte Carlo Evidence:** ARIMA(1,1,1) models show **24-52% variance reduction** with PMM2 vs CSS when innovations are non-Gaussian. --- ## Part 5: Bootstrap Inference for Time Series Bootstrap methods provide robust confidence intervals for time series parameters. ```{r ts_bootstrap_inference} # Perform bootstrap inference for the ARMA(1,1) model boot_results <- ts_pmm2_inference(fit_pmm2_arma, B = 500, seed = 123, method = "block", block_length = 20, parallel = FALSE) # Display summary with confidence intervals summary(boot_results) ``` The bootstrap results provide: - **Standard errors** for each parameter - **95% confidence intervals** (percentile and bias-corrected) - **Convergence diagnostics** across bootstrap samples --- ## Part 6: Comparing Methods with Monte Carlo The package provides tools for systematic comparison through Monte Carlo simulation. ```{r monte_carlo_comparison, eval=FALSE} # Define model specifications specs <- list( list( model = "ma", order = 1, theta = 0.7, label = "MA(1)", innovations = list(type = "gamma", shape = 2) ), list( model = "arma", order = c(1, 1), theta = list(ar = 0.5, ma = 0.4), label = "ARMA(1,1)", innovations = list(type = "student_t", df = 4) ) ) # Run Monte Carlo comparison mc_results <- pmm2_monte_carlo_compare( model_specs = specs, methods = c("css", "pmm2"), n = 300, n_sim = 1000, seed = 456, progress = FALSE ) # Display results print(mc_results$summary) print(mc_results$gain) ``` **Interpretation:** - `MSE_ratio < 1.0` indicates PMM2 has **lower variance** than CSS - Efficiency gains typically range from **10-50%** depending on innovation distribution --- ## Part 7: Prediction with Time Series Models PMM2-fitted models support forecasting: ```{r prediction} # Forecast 10 steps ahead for the ARMA(1,1) model forecast_horizon <- 10 pred_raw <- predict(fit_pmm2_arma, n.ahead = forecast_horizon) predictions <- if (is.list(pred_raw)) as.numeric(pred_raw$pred) else as.numeric(pred_raw) # Display forecasts cat("Forecasts for next", forecast_horizon, "periods:\n") print(predictions) # Plot original series with forecasts n_plot <- 100 plot_data <- tail(x, n_plot) plot(seq_along(plot_data), plot_data, type = "l", xlim = c(1, n_plot + forecast_horizon), ylim = range(c(plot_data, predictions)), main = "ARMA(1,1) Forecasts", xlab = "Time", ylab = "Value", col = "steelblue", lwd = 1.5) # Add forecasts lines(n_plot + seq_len(forecast_horizon), predictions, 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)) ``` --- ## Practical Guidelines ### When to Use PMM2 for Time Series **Recommended when:** - Innovations exhibit **skewness or heavy tails** (check residuals from classical fit) - Working with **financial data** (returns often non-Gaussian) - **Environmental time series** (precipitation, pollution levels) - **Engineering measurements** with asymmetric errors - Sample size **n ≥ 200** for reliable moment estimation **Stick with classical methods when:** - Innovations are approximately **Gaussian** (check Q-Q plot) - Very **small sample sizes** (n < 100) - **Computational speed** is critical - Model is for **pedagogical purposes** only --- ### Model Selection Strategy 1. **Visualize** the time series and check for stationarity 2. **Fit classical model** (CSS/ML) as baseline 3. **Examine residuals** for non-normality: - Skewness: `pmm_skewness(residuals)` - Kurtosis: `pmm_kurtosis(residuals)` - Visual: Q-Q plots, histograms 4. **If residuals are non-Gaussian**, refit with PMM2 5. **Compare standard errors** using bootstrap 6. **Validate** with out-of-sample forecasting --- ### Computational Considerations PMM2 is iterative and may take longer than classical methods: - **AR models**: Typically fast (5-15 iterations) - **MA models**: Moderate (10-30 iterations depending on order) - **ARMA/ARIMA**: Can be slower (20-50 iterations) For large datasets (n > 5000), consider: - Using `verbose = FALSE` to reduce output - Starting with classical estimates as `initial` parameter - Reducing `max_iter` if convergence is slow --- ## Summary Statistics: Efficiency Gains Based on extensive Monte Carlo studies documented in the package: | Model Type | Innovation Distribution | Variance Reduction | Sample Size | |------------|------------------------|-------------------|-------------| | AR(1) | Gamma (skew = 1.4) | 15-20% | 200-1000 | | MA(1) | Chi-squared (df=3) | 15-23% | 200-1000 | | MA(2) | Mixed contamination | 26-32% | 200-1000 | | ARMA(1,1) | Student-t (df=4) | 30-45% | 200-1000 | | ARIMA(1,1,1) | Exponential (shifted) | 24-52% | 300-1000 | | SAR(1,1)_12 | Gamma (shape=2) | 20-30% | 200-500 | | SMA(1)_4 | Exponential | 25-34% | 200-500 | **Key Insight:** The more **non-Gaussian** the innovations, the **larger the efficiency gain** from PMM2. --- ## Advanced Topics ### Custom Innovation Distributions You can use `pmm2_monte_carlo_compare()` with custom distributions: ```{r custom_innovations, eval=FALSE} # Custom innovation generator my_innovations <- function(n) { # Mixture: 90% Gaussian + 10% extreme outliers base <- rnorm(n) outliers <- sample(c(-5, 5), n, replace = TRUE, prob = c(0.05, 0.05)) ifelse(runif(n) < 0.9, base, outliers) } # Use in Monte Carlo specs <- list( list( model = "ar", order = 1, theta = 0.6, innovations = list(generator = my_innovations) ) ) mc_results <- pmm2_monte_carlo_compare(specs, methods = c("css", "pmm2"), n = 300, n_sim = 500) ``` --- ## Part 8: Seasonal Models (SAR and SMA) Many real-world time series exhibit seasonal patterns - monthly data with yearly cycles, quarterly economic data, etc. The EstemPMM package provides PMM2 estimation for seasonal autoregressive (SAR) and seasonal moving average (SMA) models. ### Seasonal Autoregressive (SAR) Models SAR models capture dependencies at seasonal lags. For example, SAR(1,1)_12 for monthly data relates the current value to both the previous month and the same month last year. ```{r sar_example} # Simulate SAR(1,1)_12 with seasonal period = 12 n <- 300 phi <- 0.6 # Non-seasonal AR coefficient Phi <- 0.4 # Seasonal AR coefficient (lag 12) s <- 12 # Seasonal period # Generate gamma-distributed innovations (right-skewed) innovations <- rgamma(n, shape = 2, rate = 2) - 1 # Generate SAR series y <- numeric(n) for (t in (s+2):n) { y[t] <- phi * y[t-1] + Phi * y[t-s] + innovations[t] } # Plot the series plot(y, type = "l", main = "SAR(1,1)_12 Process with Skewed Innovations", ylab = "Value", xlab = "Time", col = "steelblue", lwd = 1.5) # Fit SAR model with PMM2 fit_sar_pmm2 <- sar_pmm2( y, order = c(1, 1), season = list(period = 12), method = "pmm2" ) # Compare with OLS fit_sar_ols <- sar_pmm2( y, order = c(1, 1), season = list(period = 12), method = "ols" ) # Display results cat("True parameters: phi =", phi, ", Phi =", Phi, "\n\n") cat("OLS Estimates:\n") print(coef(fit_sar_ols)) cat("\nPMM2 Estimates:\n") print(coef(fit_sar_pmm2)) # Summary with diagnostics summary(fit_sar_pmm2) ``` **Performance:** SAR models with PMM2 show **20-30% variance reduction** compared to OLS when innovations are asymmetric. ### Seasonal Moving Average (SMA) Models SMA models use seasonal lags in the moving average component. ```{r sma_example} # Simulate SMA(1)_4 (quarterly seasonal pattern) n <- 200 Theta <- 0.6 # Seasonal MA coefficient s <- 4 # Quarterly data # Generate exponentially distributed innovations innovations <- rexp(n, rate = 1) - 1 # Generate SMA series y <- numeric(n) for (t in 1:s) { y[t] <- innovations[t] } for (t in (s+1):n) { y[t] <- innovations[t] + Theta * innovations[t-s] } # Plot the series plot(y, type = "l", main = "SMA(1)_4 Process", ylab = "Value", xlab = "Time", col = "darkgreen", lwd = 1.5) # Fit SMA model with PMM2 fit_sma_pmm2 <- sma_pmm2( y, order = 1, season = list(period = 4), method = "pmm2" ) # Fit with CSS for comparison fit_sma_css <- sma_pmm2( y, order = 1, season = list(period = 4), method = "css" ) # Display results cat("True parameter: Theta =", Theta, "\n\n") cat("CSS Estimate:\n") print(coef(fit_sma_css)) cat("\nPMM2 Estimate:\n") print(coef(fit_sma_pmm2)) cat("\nConvergence:", fit_sma_pmm2@convergence, "\n") cat("Iterations:", fit_sma_pmm2@iterations, "\n") ``` **Performance:** Extensive Monte Carlo validation (500 replications) shows SMA models achieve up to **34.1% variance reduction** with PMM2, often exceeding theoretical predictions for certain parameter combinations. ### Comparing Seasonal Methods ```{r compare_seasonal, eval=FALSE} # Systematic comparison of SAR estimation methods compare_sar_methods(y, order = c(1, 1), period = 12) # Universal comparison for non-seasonal models compare_ts_methods( y, model_type = "arima", order = c(1, 0, 1) ) ``` ### Practical Applications of Seasonal Models **Recommended for:** - **Monthly economic data** (sales, production, employment) with annual patterns - **Quarterly financial data** (GDP, earnings reports) - **Energy consumption** with daily, weekly, or seasonal cycles - **Climate data** (temperature, rainfall) with annual seasonality - **Retail sales** with holiday and seasonal effects --- ## Conclusion The **PMM2 method for time series** provides: 1. **Robust estimation** for AR, MA, ARMA, ARIMA, SAR, and SMA models 2. **10-50% variance reduction** when innovations deviate from normality 3. **Seasonal pattern modeling** with SAR/SMA for data with periodic cycles 4. **Seamless integration** with standard R workflows 5. **Comprehensive diagnostics** and bootstrap inference For practitioners working with **real-world time series** exhibiting non-Gaussian behavior, PMM2 offers a principled approach to achieving **more precise parameter estimates** without sacrificing interpretability. --- ## Next Steps - **Bootstrap Inference**: See `vignette("03-bootstrap-inference")` for detailed statistical inference - **Linear Regression**: See `vignette("01-pmm2-introduction")` for PMM2 basics - **Advanced Applications**: Explore package demos with `demo(package = "EstemPMM")` --- ## References 1. Zabolotnii, S., Warsza, Z.L., Tkachenko, O. (2018). "Polynomial Estimation of Linear Regression Parameters for the Asymmetric PDF of Errors". Springer, AUTOMATION 2018. 2. Package functions: `?ar_pmm2`, `?ma_pmm2`, `?arma_pmm2`, `?arima_pmm2`, `?ts_pmm2_inference` 3. Monte Carlo tools: `?pmm2_monte_carlo_compare` 4. GitHub: https://github.com/SZabolotnii/EstemPMM