This research models the S&P 500’s underlying volatility states using a 2-State Hidden Markov Model (HMM). By assuming the market operates in hidden, non-stationary regimes, we can mathematically isolate periods of low-volatility “bull” markets from high-volatility “panic” phases without relying on lagging moving averages.
Three strategies are evaluated out-of-sample using a walk-forward framework:
Key Finding: The baseline strategy (A) failed out-of-sample, spending only ~20% of the time invested due to crisis-scarred training data. Strategies B and C are designed to correct this regime persistence mismatch.
We fetch 20 years of daily tick data for the SPY ETF to capture multiple major market shocks. We calculate daily logarithmic returns to achieve stationarity for the HMM.
library(quantmod)
library(depmixS4)
library(ggplot2)
library(dplyr)
library(scales)
library(tidyr)
# Fetch data
getSymbols("SPY", from = "2004-01-01", to = Sys.Date(), warnings = FALSE)
spy_returns <- dailyReturn(SPY)
# Format for the model
df_returns <- data.frame(Return = as.numeric(spy_returns),
Date = index(spy_returns))
df_returns$Year <- as.integer(format(df_returns$Date, "%Y"))A 2-state Gaussian HMM is fitted to the full dataset for regime visualization purposes. The Baum-Welch algorithm (EM) estimates the transition matrix and emission parameters; the Viterbi algorithm then decodes the most likely state sequence.
set.seed(42)
hmm_model <- depmix(Return ~ 1,
data = df_returns,
nstates = 2,
family = gaussian())
fitted_hmm <- fit(hmm_model, verbose = FALSE)converged at iteration 36 with logLik: 18252.24 print("Model Successfully Converged: Viterbi Algorithm Applied.")[1] "Model Successfully Converged: Viterbi Algorithm Applied."regime_probabilities <- posterior(fitted_hmm)
df_returns$Regime <- as.factor(regime_probabilities$state)
ggplot(df_returns, aes(x = Date, y = Return)) +
geom_line(color = "gray50", alpha = 0.3, linewidth = 0.5) +
geom_point(aes(color = Regime), alpha = 0.8, size = 1.2) +
scale_color_manual(values = c("1" = "#00b894", "2" = "#ff7675"),
labels = c("Regime 1: Normal", "Regime 2: Panic Phase")) +
scale_y_continuous(labels = percent_format(accuracy = 1)) +
theme_minimal(base_size = 14) +
theme(
legend.position = "bottom",
legend.title = element_blank(),
panel.grid.minor = element_blank(),
panel.grid.major = element_line(color = "gray80", linewidth = 0.2)
) +
labs(title = "S&P 500 Market Regimes",
subtitle = "HMM Diagnosis of Underlying Volatility States",
y = "Daily Return", x = "")
All three strategies share the same walk-forward discipline: the model is only ever fitted on data prior to the test year. A 1-day signal lag is applied throughout to eliminate look-ahead bias. Transaction costs of 0.1% per switch are charged.
# --- Shared parameters ---
first_test_year <- 2010
cost_per_trade <- 0.001
rolling_window <- 3 # years of history for Strategy B
all_years <- sort(unique(df_returns$Year))
test_years <- all_years[all_years >= first_test_year]
# Helper: fit HMM on train_data, predict on test_data
# soft = FALSE -> hard binary position (0 or 1)
# soft = TRUE -> position = posterior probability of normal state (0 to 1)
run_fold <- function(train_data, test_data, soft = FALSE) {
hmm_fold <- tryCatch({
m <- depmix(Return ~ 1, data = train_data, nstates = 2, family = gaussian())
fit(m, verbose = FALSE)
}, error = function(e) NULL)
if (is.null(hmm_fold)) return(NULL)
# Identify normal state = lower volatility state
params <- getpars(hmm_fold)
sd_state1 <- params[7]
sd_state2 <- params[9]
normal_state <- ifelse(sd_state1 < sd_state2, 1, 2)
test_hmm <- tryCatch({
m_test <- depmix(Return ~ 1, data = test_data, nstates = 2, family = gaussian())
setpars(m_test, getpars(hmm_fold))
}, error = function(e) NULL)
if (is.null(test_hmm)) return(NULL)
post <- tryCatch(posterior(test_hmm), error = function(e) NULL)
if (is.null(post)) return(NULL)
if (soft) {
# Use posterior probability of normal state as continuous position size
normal_col <- ifelse(normal_state == 1, "S1", "S2")
test_data$Signal <- post[[normal_col]]
} else {
# Hard: 1 if normal state, 0 if panic state
test_data$Regime <- as.integer(post$state)
test_data$Signal <- ifelse(test_data$Regime == normal_state, 1L, 0L)
}
# 1-day lag: today's signal becomes tomorrow's position
test_data$Position <- c(1, test_data$Signal[-nrow(test_data)])
# Charge transaction cost proportional to position change magnitude
test_data$Turnover <- c(0, abs(diff(test_data$Position)))
test_data$Net_Return <- test_data$Return * test_data$Position -
test_data$Turnover * cost_per_trade
test_data
}The HMM trains on all available history before each test year. The signal is binary: fully invested or fully in cash.
wf_a <- list()
for (test_yr in test_years) {
train_data <- df_returns[df_returns$Year < test_yr, ]
test_data <- df_returns[df_returns$Year == test_yr, ]
if (nrow(train_data) < 200 || nrow(test_data) < 5) next
res <- run_fold(train_data, test_data, soft = FALSE)
if (!is.null(res)) wf_a[[as.character(test_yr)]] <- res
}converged at iteration 24 with logLik: 4778.576 converged at iteration 28 with logLik: 5558.361 converged at iteration 27 with logLik: 6304.838 converged at iteration 29 with logLik: 7152.113 converged at iteration 31 with logLik: 8042.015 converged at iteration 35 with logLik: 8929.911 converged at iteration 34 with logLik: 9753.505 converged at iteration 37 with logLik: 10620.71 converged at iteration 42 with logLik: 11587.4 converged at iteration 42 with logLik: 12410.99 converged at iteration 45 with logLik: 13291.89 converged at iteration 36 with logLik: 13968.26 converged at iteration 36 with logLik: 14820.34 converged at iteration 38 with logLik: 15487.95 converged at iteration 36 with logLik: 16319.24 converged at iteration 36 with logLik: 17184.99 converged at iteration 36 with logLik: 17999.12 wf_a_df <- bind_rows(wf_a) |> arrange(Date)The training window is capped at the most recent 3 years, so the model’s memory of the 2008 crisis fades as markets normalize. The signal remains binary.
wf_b <- list()
for (test_yr in test_years) {
train_data <- df_returns[df_returns$Year >= (test_yr - rolling_window) &
df_returns$Year < test_yr, ]
test_data <- df_returns[df_returns$Year == test_yr, ]
if (nrow(train_data) < 200 || nrow(test_data) < 5) next
res <- run_fold(train_data, test_data, soft = FALSE)
if (!is.null(res)) wf_b[[as.character(test_yr)]] <- res
}converged at iteration 28 with logLik: 2127.165 converged at iteration 29 with logLik: 2105.617 converged at iteration 67 with logLik: 2238.351 converged at iteration 62 with logLik: 2388.751 converged at iteration 43 with logLik: 2490.06 converged at iteration 80 with logLik: 2669.63 converged at iteration 50 with logLik: 2631.139 converged at iteration 53 with logLik: 2614.41 converged at iteration 151 with logLik: 2709.221 converged at iteration 31 with logLik: 2725.246 converged at iteration 32 with logLik: 2728.36 converged at iteration 30 with logLik: 2389.071 converged at iteration 26 with logLik: 2425.116 converged at iteration 36 with logLik: 2215.029 converged at iteration 60 with logLik: 2385.574 converged at iteration 47 with logLik: 2398.047 converged at iteration 56 with logLik: 2521.476 wf_b_df <- bind_rows(wf_b) |> arrange(Date)Instead of switching hard between 0% and 100% invested, the position size at each step equals the posterior probability that the market is in the normal (low-volatility) regime. This produces smoother, lower-turnover exposure and avoids being whipsawed by borderline regime calls.
wf_c <- list()
for (test_yr in test_years) {
train_data <- df_returns[df_returns$Year < test_yr, ]
test_data <- df_returns[df_returns$Year == test_yr, ]
if (nrow(train_data) < 200 || nrow(test_data) < 5) next
res <- run_fold(train_data, test_data, soft = TRUE)
if (!is.null(res)) wf_c[[as.character(test_yr)]] <- res
}converged at iteration 24 with logLik: 4778.576 converged at iteration 29 with logLik: 5558.361 converged at iteration 28 with logLik: 6304.838 converged at iteration 29 with logLik: 7152.113 converged at iteration 30 with logLik: 8042.015 converged at iteration 34 with logLik: 8929.911 converged at iteration 34 with logLik: 9753.505 converged at iteration 37 with logLik: 10620.71 converged at iteration 42 with logLik: 11587.4 converged at iteration 42 with logLik: 12410.99 converged at iteration 45 with logLik: 13291.89 converged at iteration 36 with logLik: 13968.26 converged at iteration 37 with logLik: 14820.34 converged at iteration 38 with logLik: 15487.95 converged at iteration 36 with logLik: 16319.24 converged at iteration 37 with logLik: 17184.99 converged at iteration 37 with logLik: 17999.12 wf_c_df <- bind_rows(wf_c) |> arrange(Date)# Align B and C to Strategy A's dates (they may differ at edges)
b_net <- wf_b_df$Net_Return[match(wf_a_df$Date, wf_b_df$Date)]
c_net <- wf_c_df$Net_Return[match(wf_a_df$Date, wf_c_df$Date)]
plot_df <- data.frame(
Date = wf_a_df$Date,
SPY = cumprod(1 + wf_a_df$Return),
Strat_A = cumprod(1 + wf_a_df$Net_Return),
Strat_B = cumprod(1 + ifelse(is.na(b_net), 0, b_net)),
Strat_C = cumprod(1 + ifelse(is.na(c_net), 0, c_net))
)
plot_long <- pivot_longer(plot_df, -Date, names_to = "Strategy", values_to = "Wealth")
plot_long$Strategy <- factor(plot_long$Strategy,
levels = c("SPY", "Strat_A", "Strat_B", "Strat_C"),
labels = c("Buy & Hold (SPY)",
"A: Expanding + Hard Switch",
"B: Rolling 3yr + Hard Switch",
"C: Expanding + Soft Allocation"))
ggplot(plot_long, aes(x = Date, y = Wealth, color = Strategy, linewidth = Strategy)) +
geom_line(alpha = 0.9) +
scale_color_manual(values = c(
"Buy & Hold (SPY)" = "gray50",
"A: Expanding + Hard Switch" = "#d63031",
"B: Rolling 3yr + Hard Switch" = "#e17055",
"C: Expanding + Soft Allocation" = "#0984e3"
)) +
scale_linewidth_manual(values = c(0.8, 1.0, 1.0, 1.3)) +
scale_y_continuous(labels = dollar_format()) +
theme_minimal(base_size = 14) +
theme(
legend.position = "bottom",
legend.title = element_blank(),
panel.grid.minor = element_blank(),
panel.grid.major = element_line(color = "gray85", linewidth = 0.2),
plot.title = element_text(face = "bold", size = 16)
) +
guides(linewidth = "none") +
labs(title = "Cumulative Wealth: All Strategies vs. Buy & Hold",
subtitle = paste0("Out-of-sample walk-forward (", first_test_year,
"–present) | 0.1% transaction cost applied"),
y = "Portfolio Value ($)", x = "")
trading_days <- 252
ann_return <- function(r) {
r <- r[!is.na(r)]
prod(1 + r)^(trading_days / length(r)) - 1
}
ann_vol <- function(r) sd(r, na.rm = TRUE) * sqrt(trading_days)
sharpe <- function(r, rf = 0.045) (ann_return(r) - rf) / ann_vol(r)
max_dd <- function(r) {
r <- r[!is.na(r)]
wealth <- cumprod(1 + r)
peak <- cummax(wealth)
min((wealth - peak) / peak)
}
sortino <- function(r, rf = 0.045) {
r <- r[!is.na(r)]
excess <- r - rf / trading_days
downside_r <- excess[excess < 0]
downside_vol <- sqrt(mean(downside_r^2)) * sqrt(trading_days)
(ann_return(r) - rf) / downside_vol
}
calmar <- function(r) ann_return(r) / abs(max_dd(r))
time_mkt <- function(pos) scales::percent(mean(pos, na.rm = TRUE), accuracy = 0.1)
fmt_pct <- function(x) scales::percent(x, accuracy = 0.1)
spy_r <- wf_a_df$Return
a_r <- wf_a_df$Net_Return
b_r <- wf_b_df$Net_Return[match(wf_a_df$Date, wf_b_df$Date)]
c_r <- wf_c_df$Net_Return[match(wf_a_df$Date, wf_c_df$Date)]
b_pos <- wf_b_df$Position[match(wf_a_df$Date, wf_b_df$Date)]
c_pos <- wf_c_df$Position[match(wf_a_df$Date, wf_c_df$Date)]
perf_table <- data.frame(
Metric = c(
"Annualized Return",
"Annualized Volatility",
"Sharpe Ratio (rf = 4.5%)",
"Sortino Ratio",
"Calmar Ratio",
"Max Drawdown",
"Avg. Time in Market"
),
`Buy & Hold` = c(
fmt_pct(ann_return(spy_r)), fmt_pct(ann_vol(spy_r)),
round(sharpe(spy_r), 2), round(sortino(spy_r), 2),
round(calmar(spy_r), 2), fmt_pct(max_dd(spy_r)), "100%"
),
`A: Expanding + Hard` = c(
fmt_pct(ann_return(a_r)), fmt_pct(ann_vol(a_r)),
round(sharpe(a_r), 2), round(sortino(a_r), 2),
round(calmar(a_r), 2), fmt_pct(max_dd(a_r)),
time_mkt(wf_a_df$Position)
),
`B: Rolling 3yr + Hard` = c(
fmt_pct(ann_return(b_r)), fmt_pct(ann_vol(b_r)),
round(sharpe(b_r), 2), round(sortino(b_r), 2),
round(calmar(b_r), 2), fmt_pct(max_dd(b_r)),
time_mkt(b_pos)
),
`C: Expanding + Soft` = c(
fmt_pct(ann_return(c_r)), fmt_pct(ann_vol(c_r)),
round(sharpe(c_r), 2), round(sortino(c_r), 2),
round(calmar(c_r), 2), fmt_pct(max_dd(c_r)),
time_mkt(c_pos)
),
check.names = FALSE
)
knitr::kable(perf_table, align = c("l", "r", "r", "r", "r"),
caption = "Out-of-sample performance comparison (0.1% transaction costs, 2010–present)")| Metric | Buy & Hold | A: Expanding + Hard | B: Rolling 3yr + Hard | C: Expanding + Soft |
|---|---|---|---|---|
| Annualized Return | 12.1% | -3.4% | -2.7% | 2.1% |
| Annualized Volatility | 17.2% | 13.7% | 13.7% | 12.3% |
| Sharpe Ratio (rf = 4.5%) | 0.44 | -0.58 | -0.52 | -0.19 |
| Sortino Ratio | 0.42 | -0.75 | -0.67 | -0.24 |
| Calmar Ratio | 0.35 | -0.07 | -0.05 | 0.07 |
| Max Drawdown | -34.1% | -52.4% | -49.1% | -30.5% |
| Avg. Time in Market | 100% | 19.8% | 24.8% | 20.9% |
The baseline strategy (A) failed out-of-sample, spending only ~20% of the period invested. This is a textbook example of regime persistence mismatch: training on a crisis-heavy window (2004–2009) caused the model to systematically over-identify panic regimes during the subsequent bull market, resulting in a -3.4% annualized return and a worse maximum drawdown than buy & hold.
Strategy B corrects the memory problem by restricting training to the most recent 3 years. Without the weight of 2008 dominating the prior, the model recalibrates to prevailing market conditions and spends more time invested in normal regimes.
Strategy C addresses the problem differently: rather than fixing the training window, it removes the discretization step entirely. By using the posterior probability of a normal regime as a continuous position size, it avoids whipsaw exits caused by borderline regime calls and substantially reduces turnover costs.
Research Implication: For regime-switching models applied to trending markets, training window length and signal discretization are first-order design decisions — not hyperparameters to tune after the fact. A model that correctly identifies regimes but carries the wrong historical prior will still fail in production. The comparison across A, B, and C isolates each source of failure independently.