Part 1: Seasonal Autoregressive (SAR) Models
Understanding SAR Models
A **SAR(p,P)_s** model with seasonal period s is:
\[ y_t = \phi_1 y_{t-1} + \cdots + \phi_p y_{t-p} + \Phi_1 y_{t-s} + \cdots + \Phi_P y_{t-Ps} + \varepsilon_t \]
where: - \(\phi_i\) are non-seasonal AR coefficients - \(\Phi_j\) are seasonal AR coefficients - \(s\) is the seasonal period (e.g., 12 for monthly data, 4 for quarterly)
Example: SAR(1,1)_12 with Monthly Data
# Simulate SAR(1,1)_12 model
n <- 360 # 30 years of monthly data
phi <- 0.5 # Non-seasonal AR(1)
Phi <- 0.6 # Seasonal AR(1) at lag 12
s <- 12 # Monthly seasonal period
# Generate gamma-distributed innovations (right-skewed)
# This mimics asymmetric economic shocks
innovations <- rgamma(n + 100, shape = 2, scale = 1) - 2
# Generate SAR series
y <- numeric(n + 100)
for (t in (s + 2):(n + 100)) {
y[t] <- phi * y[t-1] + Phi * y[t-s] + innovations[t]
}
y <- y[101:(n + 100)] # Remove burn-in
# Convert to time series object for better plotting
y_ts <- ts(y, frequency = 12, start = c(2000, 1))
# Plot the series
plot(y_ts, main = "SAR(1,1)_12 Process with Skewed Innovations",
ylab = "Value", xlab = "Year", col = "steelblue", lwd = 1.5)
abline(h = 0, col = "gray", lty = 2)
Seasonal Patterns Visualization
# Extract seasonal patterns
months <- rep(1:12, length.out = length(y))
seasonal_df <- data.frame(
Month = factor(months, labels = month.abb),
Value = y
)
# Boxplot by month
boxplot(Value ~ Month, data = seasonal_df,
main = "Seasonal Pattern in SAR(1,1)_12 Data",
xlab = "Month", ylab = "Value",
col = "lightblue", border = "steelblue")
Fitting SAR Models: PMM2 vs Classical Methods
# Fit SAR(1,1)_12 with PMM2
fit_sar_pmm2 <- sar_pmm2(
y,
order = c(1, 1), # (non-seasonal AR order, seasonal AR order)
season = list(period = 12), # Monthly seasonality
method = "pmm2",
verbose = FALSE
)
# Fit with OLS for comparison
fit_sar_ols <- sar_pmm2(
y,
order = c(1, 1),
season = list(period = 12),
method = "ols",
verbose = FALSE
)
# Fit with CSS for comparison
fit_sar_css <- sar_pmm2(
y,
order = c(1, 1),
season = list(period = 12),
method = "css",
verbose = FALSE
)
# Compare estimates
cat("True parameters:\n")
#> True parameters:
cat(" phi (AR1):", phi, "\n")
#> phi (AR1): 0.5
cat(" Phi (SAR1):", Phi, "\n\n")
#> Phi (SAR1): 0.6
cat("OLS Estimates:\n")
#> OLS Estimates:
print(coef(fit_sar_ols))
#> ar1 sar1
#> 0.5121056 0.5860107
cat("\nCSS Estimates:\n")
#>
#> CSS Estimates:
print(coef(fit_sar_css))
#> ar1 sar1
#> 0.5121056 0.5860107
cat("\nPMM2 Estimates:\n")
#>
#> PMM2 Estimates:
print(coef(fit_sar_pmm2))
#> ar1 sar1
#> 0.4937866 0.6072927
# Display PMM2 summary with diagnostics
cat("\n")
summary(fit_sar_pmm2)
#>
#> =======================================================================
#> Seasonal AR Model fitted with PMM2
#> =======================================================================
#>
#> Call:
#> sar_pmm2(x = y, order = c(1, 1), season = list(period = 12),
#> method = "pmm2", verbose = FALSE)
#>
#> Model: SAR(1,1)_12
#> Observations: 360
#> Effective sample size: 360
#>
#> Coefficients:
#> -----------------------------------------
#> ar1 : 0.49379
#> sar1 : 0.60729
#> Mean : 0.09881
#>
#> Residual Distribution Characteristics:
#> -----------------------------------------
#> m2 (variance): 2.44677
#> m3 (skewness ind.): 6.05862
#> m4 (kurtosis ind.): 38.73868
#>
#> Skewness coef. (c3): 1.58301
#> Excess kurtosis(c4): 3.47083
#> Var. reduction (g): 0.54195
#>
#> => Expected 45.8% variance reduction vs OLS [OK]
#>
#> Algorithm Information:
#> -----------------------------------------
#> Converged: TRUE
#> Iterations: 3
#>
#> Residual Statistics:
#> -----------------------------------------
#> Min: -2.0082
#> Q1: -1.1558
#> Median: -0.3302
#> Mean: 0.0013
#> Q3: 0.7034
#> Max: 8.1440
#>
#> =======================================================================Understanding PMM2 Efficiency
The PMM2 method computes the variance correction factor \(g\):
# Extract moments from fit
m2 <- fit_sar_pmm2@m2
m3 <- fit_sar_pmm2@m3
m4 <- fit_sar_pmm2@m4
# Compute g-factor
g_result <- pmm2_variance_factor(m2, m3, m4)
cat("PMM2 variance correction factor:\n")
#> PMM2 variance correction factor:
cat(" g =", round(g_result$g, 4), "\n")
#> g = 0.5419
cat(" Theoretical variance reduction:",
round((1 - g_result$g) * 100, 2), "%\n")
#> Theoretical variance reduction: 45.81 %
# For Gaussian innovations, g ≈ 1 (no improvement)
# For skewed/heavy-tailed innovations, g < 1 (PMM2 more efficient)Residual Diagnostics
# Create diagnostic plots
par(mfrow = c(2, 2))
# 1. Residuals over time
residuals_pmm2 <- fit_sar_pmm2@residuals
plot(residuals_pmm2, type = "l",
main = "PMM2 Residuals",
xlab = "Time", ylab = "Residual",
col = "steelblue")
abline(h = 0, col = "red", lty = 2)
# 2. ACF of residuals
acf(residuals_pmm2, main = "ACF of PMM2 Residuals",
col = "steelblue", lwd = 2)
# 3. Q-Q plot
qqnorm(residuals_pmm2, main = "Q-Q Plot of PMM2 Residuals",
pch = 20, col = "steelblue")
qqline(residuals_pmm2, col = "red", lwd = 2)
# 4. Histogram of residuals
hist(residuals_pmm2, breaks = 30, col = "lightblue",
border = "steelblue", main = "Distribution of PMM2 Residuals",
xlab = "Residual", freq = FALSE)
lines(density(residuals_pmm2), col = "darkblue", lwd = 2)
