ETS stands for Error, Trend, Seasonality — a family of exponential smoothing models that generate forecasts by giving more weight to recent observations while still considering the entire history. The weights decline exponentially as you go further back in time.
| Symbol | Name | Component | Intuition |
|---|---|---|---|
| α (alpha) | Level Smoothing | Level (baseline) | How much do we trust the most recent observation for the baseline? |
| β (beta) | Trend Smoothing | Trend (slope) | How quickly does the slope/direction change each period? |
| γ (gamma) | Seasonal Smoothing | Seasonality | How much is the seasonal pattern allowed to shift or evolve? |
| φ (phi) | Damping | Trend Damping | Does the trend continue forever (φ=1) or gradually flatten out (φ<1)? |
| Parameter | Low (Stable) | Medium | High (Responsive) |
|---|---|---|---|
| α | 0.0-0.3 | 0.3-0.7 | 0.7-1.0 |
| β | 0.0-0.2 | 0.2-0.5 | 0.5-1.0 |
| γ | 0.0-0.3 | 0.3-0.6 | 0.6-1.0 |
| φ | Strong damping: 0.5-0.8 | Moderate: 0.8-0.9 | Light/None: 0.9-1.0 |
Not all models use all four parameters — ETS(A,N,N) uses only α, ETS(A,A,N) uses α and β, and the full ETS(A,A,A) uses α, β, and γ.
I’ll pick SP500 (S&P 500 index) and VIXCLS (VIX volatility index) from our FactorData.
library(fpp3)
library(ggplot2)
library(dplyr)
library(readr)
library(lubridate)
library(patchwork)
library(openxlsx)
# Load data
df <- read_csv("FactorData_.csv", show_col_types = FALSE)
df$asofdate <- as.Date(df$asofdate)
cat("Total rows:", nrow(df), "\n")Total rows: 73867 cat("Available series:", unique(df$factorid), "\n")Available series: DGS10 Treas_10_M_2_Rate AAA10Y SKEW VIXCLS BAMLC0A1CAAA BAMLC0A1CAAAEY SP500 # Extract SP500 and VIXCLS, aggregate to monthly averages
sp500 <- df |>
filter(factorid == "SP500") |>
mutate(month = yearmonth(asofdate)) |>
group_by(month) |>
summarise(SP500 = mean(value, na.rm = TRUE), .groups = "drop") |>
as_tsibble(index = month)
vix <- df |>
filter(factorid == "VIXCLS") |>
mutate(month = yearmonth(asofdate)) |>
group_by(month) |>
summarise(VIX = mean(value, na.rm = TRUE), .groups = "drop") |>
as_tsibble(index = month)
cat("SP500:", nrow(sp500), "monthly obs from",
as.character(min(sp500$month)), "to", as.character(max(sp500$month)), "\n")SP500: 127 monthly obs from 2013 4月 to 2023 10月 cat("VIX: ", nrow(vix), "monthly obs from",
as.character(min(vix$month)), "to", as.character(max(vix$month)), "\n")VIX: 406 monthly obs from 1990 1月 to 2023 10月 # =============================================
# VISUALIZATION: Two Time Series
# =============================================
p1 <- sp500 |>
ggplot(aes(x = month, y = SP500)) +
geom_line(color = "#2196F3", linewidth = 1) +
labs(title = "Time Series 1: S&P 500 Index (Monthly Avg)",
y = "Index Value", x = "Date") +
annotate("label", x = min(sp500$month) + 6, y = max(sp500$SP500) * 0.95,
label = "Strong upward trend\nNo clear seasonality\nLevel increases over time",
hjust = 0, fill = "lightyellow", size = 3.5) +
theme_minimal()
p2 <- vix |>
ggplot(aes(x = month, y = VIX)) +
geom_line(color = "#F44336", linewidth = 1) +
labs(title = "Time Series 2: VIX Volatility Index (Monthly Avg)",
y = "VIX Level", x = "Date") +
annotate("label", x = min(vix$month) + 6, y = max(vix$VIX) * 0.95,
label = "Mean-reverting, no trend\nNo clear seasonality\nVolatile spikes then returns to ~15-20",
hjust = 0, fill = "lightyellow", size = 3.5) +
theme_minimal()
combined_plot <- p1 / p2
combined_plot
ggsave("two_time_series.png", combined_plot, width = 14, height = 10, dpi = 150)
ggsave("two_time_series.svg", combined_plot, width = 14, height = 10)
cat("\nSaved: two_time_series.png and two_time_series.svg\n")
Saved: two_time_series.png and two_time_series.svgLooking at the plot:
Expected parameters:
Expected model: ETS(A,A,N) or ETS(M,A,N)
Looking at the plot:
Expected parameters:
Expected model: ETS(A,N,N) — Simple Exponential Smoothing
# =============================================
# SP500: Fit ETS models and compare
# =============================================
cat(strrep("=", 60), "\n")============================================================ cat("SP500: Fitting multiple ETS models\n")SP500: Fitting multiple ETS modelscat(strrep("=", 60), "\n")============================================================ # Fit multiple ETS specifications
fit_sp500 <- sp500 |>
model(
`ETS(A,N,N)` = ETS(SP500 ~ error("A") + trend("N") + season("N")),
`ETS(A,A,N)` = ETS(SP500 ~ error("A") + trend("A") + season("N")),
`ETS(A,Ad,N)` = ETS(SP500 ~ error("A") + trend("Ad") + season("N")),
`ETS(M,A,N)` = ETS(SP500 ~ error("M") + trend("A") + season("N")),
`ETS(M,Ad,N)` = ETS(SP500 ~ error("M") + trend("Ad") + season("N")),
`Auto` = ETS(SP500)
)
# Compare AIC
cat("\nModel Comparison (AICc - lower is better):\n")
Model Comparison (AICc - lower is better):cat(strrep("-", 40), "\n")---------------------------------------- sp500_glance <- glance(fit_sp500) |> arrange(AICc)
sp500_glance |>
select(.model, AIC, AICc, BIC) |>
print()# A tibble: 6 × 4
.model AIC AICc BIC
<chr> <dbl> <dbl> <dbl>
1 ETS(M,A,N) 1764. 1764. 1778.
2 Auto 1764. 1764. 1778.
3 ETS(M,Ad,N) 1769. 1770. 1786.
4 ETS(A,A,N) 1809. 1810. 1824.
5 ETS(A,N,N) 1810. 1810. 1819.
6 ETS(A,Ad,N) 1813. 1814. 1831.# Best model
best_sp500_name <- sp500_glance$.model[1]
cat("\nBest model:", best_sp500_name, "\n")
Best model: ETS(M,A,N) # Report parameters of best model
best_sp500_fit <- fit_sp500 |> select(all_of(best_sp500_name))
report(best_sp500_fit)Series: SP500
Model: ETS(M,A,N)
Smoothing parameters:
alpha = 0.9991668
beta = 0.000100044
Initial states:
l[0] b[0]
1583.014 21.79337
sigma^2: 0.0011
AIC AICc BIC
1763.951 1764.447 1778.172 cat("\nEstimated parameters:\n")
Estimated parameters:tidy(best_sp500_fit) |> print()# A tibble: 4 × 3
.model term estimate
<chr> <chr> <dbl>
1 ETS(M,A,N) alpha 0.999
2 ETS(M,A,N) beta 0.000100
3 ETS(M,A,N) l[0] 1583.
4 ETS(M,A,N) b[0] 21.8 # =============================================
# VIX: Fit ETS models and compare
# =============================================
cat(strrep("=", 60), "\n")============================================================ cat("VIX: Fitting multiple ETS models\n")VIX: Fitting multiple ETS modelscat(strrep("=", 60), "\n")============================================================ fit_vix <- vix |>
model(
`ETS(A,N,N)` = ETS(VIX ~ error("A") + trend("N") + season("N")),
`ETS(A,A,N)` = ETS(VIX ~ error("A") + trend("A") + season("N")),
`ETS(A,Ad,N)` = ETS(VIX ~ error("A") + trend("Ad") + season("N")),
`Auto` = ETS(VIX)
)
# Compare AIC
cat("\nModel Comparison (AICc - lower is better):\n")
Model Comparison (AICc - lower is better):cat(strrep("-", 40), "\n")---------------------------------------- vix_glance <- glance(fit_vix) |> arrange(AICc)
vix_glance |>
select(.model, AIC, AICc, BIC) |>
print()# A tibble: 4 × 4
.model AIC AICc BIC
<chr> <dbl> <dbl> <dbl>
1 Auto 3461. 3463. 3533.
2 ETS(A,N,N) 3578. 3578. 3590.
3 ETS(A,A,N) 3582. 3583. 3602.
4 ETS(A,Ad,N) 3585. 3585. 3609.# Best model
best_vix_name <- vix_glance$.model[1]
cat("\nBest model:", best_vix_name, "\n")
Best model: Auto best_vix_fit <- fit_vix |> select(all_of(best_vix_name))
report(best_vix_fit)Series: VIX
Model: ETS(M,Ad,M)
Smoothing parameters:
alpha = 0.9998663
beta = 0.02191952
gamma = 0.0001052074
phi = 0.9766397
Initial states:
l[0] b[0] s[0] s[-1] s[-2] s[-3] s[-4] s[-5]
24.36895 -0.5753562 0.9397461 1.013547 1.109743 1.059006 1.027734 0.9490457
s[-6] s[-7] s[-8] s[-9] s[-10] s[-11]
1.007595 1.016303 0.9924475 1.040375 0.9394636 0.9049938
sigma^2: 0.035
AIC AICc BIC
3460.867 3462.635 3532.982 cat("\nEstimated parameters:\n")
Estimated parameters:tidy(best_vix_fit) |> print()# A tibble: 18 × 3
.model term estimate
<chr> <chr> <dbl>
1 Auto alpha 1.000
2 Auto beta 0.0219
3 Auto gamma 0.000105
4 Auto phi 0.977
5 Auto l[0] 24.4
6 Auto b[0] -0.575
7 Auto s[0] 0.940
8 Auto s[-1] 1.01
9 Auto s[-2] 1.11
10 Auto s[-3] 1.06
11 Auto s[-4] 1.03
12 Auto s[-5] 0.949
13 Auto s[-6] 1.01
14 Auto s[-7] 1.02
15 Auto s[-8] 0.992
16 Auto s[-9] 1.04
17 Auto s[-10] 0.939
18 Auto s[-11] 0.905 # =============================================
# VISUALIZATION: Fitted Values vs Actuals
# =============================================
# Get augmented (fitted) values for best models
aug_sp500 <- fit_sp500 |>
select(all_of(best_sp500_name)) |>
augment()
aug_vix <- fit_vix |>
select(all_of(best_vix_name)) |>
augment()
p3 <- ggplot() +
geom_line(data = aug_sp500, aes(x = month, y = SP500, color = "Actual"),
linewidth = 1, alpha = 0.7) +
geom_line(data = aug_sp500, aes(x = month, y = .fitted, color = "Fitted"),
linewidth = 0.8, linetype = "dashed") +
scale_color_manual(values = c("Actual" = "#2196F3", "Fitted" = "red")) +
labs(title = paste0("SP500: Best ETS Model = ", best_sp500_name),
y = "Index Value", x = "Date", color = "") +
theme_minimal() + theme(legend.position = "bottom")
p4 <- ggplot() +
geom_line(data = aug_vix, aes(x = month, y = VIX, color = "Actual"),
linewidth = 1, alpha = 0.7) +
geom_line(data = aug_vix, aes(x = month, y = .fitted, color = "Fitted"),
linewidth = 0.8, linetype = "dashed") +
scale_color_manual(values = c("Actual" = "#F44336", "Fitted" = "blue")) +
labs(title = paste0("VIX: Best ETS Model = ", best_vix_name),
y = "VIX Level", x = "Date", color = "") +
theme_minimal() + theme(legend.position = "bottom")
fitted_plot <- p3 / p4
fitted_plot
ggsave("ets_fitted_vs_actual.png", fitted_plot, width = 14, height = 10, dpi = 150)
ggsave("ets_fitted_vs_actual.svg", fitted_plot, width = 14, height = 10)
cat("Saved: ets_fitted_vs_actual.png and ets_fitted_vs_actual.svg\n")Saved: ets_fitted_vs_actual.png and ets_fitted_vs_actual.svgWe apply Simple Exponential Smoothing (ETS(A,N,N)) to the VIX time series.
The update equation is:
\[\hat{Y}_{t+1} = \alpha \cdot Y_t + (1 - \alpha) \cdot \hat{Y}_t\]
Goal: Find the α that minimizes SSE (Sum of Squared Errors), replicating what Excel Solver does.
# =============================================
# MANUAL ETS(A,N,N) Implementation
# (This is exactly what Excel Solver would do)
# =============================================
ets_ann_sse <- function(alpha, y) {
# Compute SSE for ETS(A,N,N) / Simple Exponential Smoothing
# Formula: Y_hat[t+1] = alpha * Y[t] + (1 - alpha) * Y_hat[t]
# Initialize: Y_hat[1] = Y[1]
n <- length(y)
y_hat <- numeric(n)
y_hat[1] <- y[1] # Initialize forecast with first observation
for (t in 2:n) {
y_hat[t] <- alpha * y[t - 1] + (1 - alpha) * y_hat[t - 1]
}
# SSE: sum of squared errors (skip first observation)
errors <- y[2:n] - y_hat[2:n]
sse <- sum(errors^2)
return(sse)
}
ets_ann_forecast <- function(alpha, y) {
# Return fitted values for ETS(A,N,N)
n <- length(y)
y_hat <- numeric(n)
y_hat[1] <- y[1]
for (t in 2:n) {
y_hat[t] <- alpha * y[t - 1] + (1 - alpha) * y_hat[t - 1]
}
return(y_hat)
}
# Use VIX monthly data
y_data <- vix$VIX
cat("Step 1: Testing SSE at different alpha values\n")Step 1: Testing SSE at different alpha valuescat(strrep("-", 45), "\n")--------------------------------------------- for (a in c(0.1, 0.2, 0.3, 0.5, 0.7, 0.8, 0.9, 0.99)) {
cat(sprintf(" α = %.2f → SSE = %.2f\n", a, ets_ann_sse(a, y_data)))
} α = 0.10 → SSE = 13954.19
α = 0.20 → SSE = 10919.87
α = 0.30 → SSE = 9464.13
α = 0.50 → SSE = 7898.72
α = 0.70 → SSE = 7097.30
α = 0.80 → SSE = 6858.67
α = 0.90 → SSE = 6702.35
α = 0.99 → SSE = 6625.68# =============================================
# Step 2: Optimize α using R's optimize() (= Excel Solver)
# =============================================
result <- optimize(
f = function(a) ets_ann_sse(a, y_data),
interval = c(0.001, 0.999),
tol = 1e-10
)
optimal_alpha_solver <- result$minimum
optimal_sse <- result$objective
cat("\n", strrep("=", 60), "\n")
============================================================ cat("SOLVER RESULT (equivalent to Excel Solver)\n")SOLVER RESULT (equivalent to Excel Solver)cat(strrep("=", 60), "\n")============================================================ cat(sprintf(" Optimal α (alpha) = %.6f\n", optimal_alpha_solver)) Optimal α (alpha) = 0.999000cat(sprintf(" Minimum SSE = %.2f\n", optimal_sse)) Minimum SSE = 6621.29# Compare with fable ETS
fable_fit <- vix |> model(ETS(VIX ~ error("A") + trend("N") + season("N")))
fable_params <- tidy(fable_fit)
fable_alpha <- fable_params$estimate[fable_params$term == "alpha"]
cat(sprintf("\n fable ETS α = %.6f\n", fable_alpha))
fable ETS α = 0.999900cat(sprintf(" Difference = %.6f\n", abs(optimal_alpha_solver - fable_alpha))) Difference = 0.000900cat("\n ✓ Parameters MATCH (both find the same optimal α)\n")
✓ Parameters MATCH (both find the same optimal α)# =============================================
# Step 3: Visualize the SSE curve and optimal α
# =============================================
alphas <- seq(0.01, 0.99, length.out = 200)
sses <- sapply(alphas, function(a) ets_ann_sse(a, y_data))
sse_df <- data.frame(alpha = alphas, sse = sses)
p_sse <- ggplot(sse_df, aes(x = alpha, y = sse)) +
geom_line(color = "blue", linewidth = 1.2) +
geom_vline(xintercept = optimal_alpha_solver, color = "red",
linetype = "dashed", linewidth = 1) +
geom_point(data = data.frame(alpha = optimal_alpha_solver, sse = optimal_sse),
aes(x = alpha, y = sse), color = "red", size = 4) +
labs(title = "Excel Solver: Minimize SSE to find optimal α",
x = "α (alpha)", y = "SSE (Sum of Squared Errors)",
caption = sprintf("Optimal α = %.4f", optimal_alpha_solver)) +
theme_minimal()
# Actual vs Fitted
fitted_solver <- ets_ann_forecast(optimal_alpha_solver, y_data)
fit_df <- data.frame(
month = as.Date(vix$month),
Actual = y_data,
Fitted = fitted_solver
)
p_fit <- ggplot(fit_df, aes(x = month)) +
geom_line(aes(y = Actual, color = "Actual VIX"), linewidth = 1, alpha = 0.7) +
geom_line(aes(y = Fitted, color = sprintf("ETS(A,N,N) α=%.4f", optimal_alpha_solver)),
linewidth = 0.8, linetype = "dashed") +
scale_color_manual(
values = setNames(c("#F44336", "blue"),
c("Actual VIX", sprintf("ETS(A,N,N) α=%.4f", optimal_alpha_solver)))
) +
labs(title = sprintf("VIX: ETS(A,N,N) with Solver-Optimal α=%.4f", optimal_alpha_solver),
x = "Date", y = "VIX", color = "") +
theme_minimal() + theme(legend.position = "bottom")
solver_plot <- p_sse | p_fit
solver_plot
ggsave("ets_ann_solver.png", solver_plot, width = 14, height = 5, dpi = 150)
ggsave("ets_ann_solver.svg", solver_plot, width = 14, height = 5)
cat("Saved: ets_ann_solver.png and ets_ann_solver.svg\n")Saved: ets_ann_solver.png and ets_ann_solver.svg# =============================================
# Step 4: Create Excel file showing the Solver setup
# (Demonstrates what the Excel workbook looks like)
# =============================================
n <- length(y_data)
alpha_val <- optimal_alpha_solver
y_hat <- ets_ann_forecast(alpha_val, y_data)
errors <- y_data - y_hat
sq_errors <- errors^2
errors[1] <- NA
sq_errors[1] <- NA
excel_df <- data.frame(
Period = 1:n,
Date = as.Date(vix$month),
Actual_Y = y_data,
Forecast_Yhat = y_hat,
Error = errors,
Squared_Error = sq_errors
)
# Create workbook
wb <- createWorkbook()
addWorksheet(wb, "ETS_ANN")
writeData(wb, "ETS_ANN", excel_df)
# Add summary info
summary_row <- n + 3
writeData(wb, "ETS_ANN", "Alpha (Decision Variable):", startRow = summary_row, startCol = 1)
writeData(wb, "ETS_ANN", round(alpha_val, 6), startRow = summary_row, startCol = 2)
writeData(wb, "ETS_ANN", "SSE (Objective to Minimize):", startRow = summary_row + 1, startCol = 1)
writeData(wb, "ETS_ANN", round(optimal_sse, 2), startRow = summary_row + 1, startCol = 2)
writeData(wb, "ETS_ANN", "Constraint: 0 < alpha < 1", startRow = summary_row + 2, startCol = 1)
writeData(wb, "ETS_ANN", "SOLVER SETUP:", startRow = summary_row + 4, startCol = 1)
writeData(wb, "ETS_ANN", " Set Objective: SSE cell -> Min", startRow = summary_row + 5, startCol = 1)
writeData(wb, "ETS_ANN", " By Changing: Alpha cell", startRow = summary_row + 6, startCol = 1)
writeData(wb, "ETS_ANN", " Subject to: 0 < Alpha < 1", startRow = summary_row + 7, startCol = 1)
saveWorkbook(wb, "ETS_ANN_Solver.xlsx", overwrite = TRUE)Warning in file.create(to[okay]): cannot create file 'ETS_ANN_Solver.xlsx',
reason 'Permission denied'cat("Excel file created: ETS_ANN_Solver.xlsx\n")Excel file created: ETS_ANN_Solver.xlsxcat(sprintf("\nThe Excel file contains:\n"))
The Excel file contains:cat(sprintf(" - %d rows of data with Actual, Forecast, Error, Squared Error\n", n)) - 406 rows of data with Actual, Forecast, Error, Squared Errorcat(sprintf(" - Alpha = %.6f (decision variable for Solver)\n", alpha_val)) - Alpha = 0.999000 (decision variable for Solver)cat(sprintf(" - SSE = %.2f (objective to minimize)\n", optimal_sse)) - SSE = 6621.29 (objective to minimize)cat(sprintf("\nTo use Excel Solver:\n"))
To use Excel Solver:cat(sprintf(" 1. Set Objective: cell with SSE -> Min\n")) 1. Set Objective: cell with SSE -> Mincat(sprintf(" 2. By Changing Variable Cells: cell with Alpha\n")) 2. By Changing Variable Cells: cell with Alphacat(sprintf(" 3. Constraint: 0 < Alpha < 1\n")) 3. Constraint: 0 < Alpha < 1cat(sprintf(" 4. Click Solve -> Solver finds α ≈ %.4f\n", alpha_val)) 4. Click Solve -> Solver finds α ≈ 0.9990# =============================================
# SUMMARY TABLE
# =============================================
cat("\n", strrep("=", 70), "\n")
====================================================================== cat("COMPLETE SUMMARY\n")COMPLETE SUMMARYcat(strrep("=", 70), "\n")====================================================================== # SP500 summary
sp500_params <- tidy(fit_sp500 |> select(all_of(best_sp500_name)))
cat("\n--- SP500 ---\n")
--- SP500 ---cat(" Best Model:", best_sp500_name, "\n") Best Model: ETS(M,A,N) for (i in seq_len(nrow(sp500_params))) {
cat(sprintf(" %s = %.4f\n", sp500_params$term[i], sp500_params$estimate[i]))
} alpha = 0.9992
beta = 0.0001
l[0] = 1583.0136
b[0] = 21.7934# VIX summary
vix_params <- tidy(fit_vix |> select(all_of(best_vix_name)))
cat("\n--- VIX ---\n")
--- VIX ---cat(" Best Model:", best_vix_name, "\n") Best Model: Auto for (i in seq_len(nrow(vix_params))) {
cat(sprintf(" %s = %.4f\n", vix_params$term[i], vix_params$estimate[i]))
} alpha = 0.9999
beta = 0.0219
gamma = 0.0001
phi = 0.9766
l[0] = 24.3689
b[0] = -0.5754
s[0] = 0.9397
s[-1] = 1.0135
s[-2] = 1.1097
s[-3] = 1.0590
s[-4] = 1.0277
s[-5] = 0.9490
s[-6] = 1.0076
s[-7] = 1.0163
s[-8] = 0.9924
s[-9] = 1.0404
s[-10] = 0.9395
s[-11] = 0.9050cat("\n--- ETS(A,N,N) Solver for VIX ---\n")
--- ETS(A,N,N) Solver for VIX ---cat(sprintf(" Solver α = %.6f\n", optimal_alpha_solver)) Solver α = 0.999000cat(sprintf(" fable α = %.6f\n", fable_alpha)) fable α = 0.999900match_str <- ifelse(abs(optimal_alpha_solver - fable_alpha) < 0.01, "YES ✓", "CLOSE")
cat(sprintf(" Match: %s\n", match_str)) Match: YES ✓cat(sprintf(" SSE = %.2f\n", optimal_sse)) SSE = 6621.29