Pricing options from live data with quantmod, compared to Black-Scholes

Overview

This vignette prices European and Asian options on a real, live-pulled underlying using only empiricaloptions and its own price history – no option-quote market is used to calibrate anything. We then compare the European prices to Black-Scholes prices computed at the underlying’s own trailing historical volatility, which is the natural baseline to check against when (as empiricaloptions assumes) no implied-volatility surface is available: it isolates what the AR(1) + stationary block bootstrap machinery buys you over the simplest possible historical alternative, a plain GBM assumption fed the same estimation window.

library(empiricaloptions)
suppressPackageStartupMessages({
  library(quantmod)
  library(zoo)
})

1. Pull the underlying

We use SPY (S&P 500 ETF) as a liquid, easily-verified stand-in for the kind of underlying empiricaloptions targets – the method itself does not care whether an option market exists for it.

end_date   <- Sys.Date()
start_date <- end_date - 3 * 365   # trailing 3y of daily data

getSymbols("SPY", src = "yahoo", from = start_date, to = end_date,
           auto.assign = TRUE)
#> [1] "SPY"

spy_close <- as.numeric(Cl(SPY))
spy_dates <- index(SPY)

na_idx <- which(is.na(spy_close))
if (length(na_idx) > 0) {
  z <- na.approx(zoo(spy_close, order.by = spy_dates), maxgap = 3, na.rm = FALSE)
  spy_close <- as.numeric(z)
  spy_dates <- spy_dates[!is.na(spy_close)]
  spy_close <- spy_close[!is.na(spy_close)]
}

S0 <- tail(spy_close, 1)
cat(sprintf("Series: %d obs, %s to %s. S0 = %.2f\n",
            length(spy_close), spy_dates[1], tail(spy_dates, 1), S0))
#> Series: 750 obs, 2023-10-02 to 2026-09-28. S0 = 765.61

2. Build the empirical pool

A single eq_pool() call does the AR(1) filtering and stationary block bootstrap described in the companion paper. We use a flat short-term discount rate as a stand-in for the funding rate; substitute a real curve if you have one.

r_flat <- 0.045   # illustrative flat discount rate

pool <- eq_pool(prices = spy_close, r = r_flat, n_sims = 10000)
pool
#> <empirical residual pool>
#>   observations: 749 | paths: 10000 | rate: 0.0450
#>   AR(1) phi: -0.0642 | ann. vol: 0.1527
#>   Ljung-Box p: 0.9660 | ARCH-LM p: 2.65e-23
#>   skew: 0.5427 | excess kurtosis: 18.0155

The printed diagnostics matter here: a Ljung-Box $p$-value near 1 says the AR(1) step has removed the linear structure it was meant to remove, while a very small ARCH-LM $p$-value confirms the residuals still carry volatility clustering – exactly the dependence the block (rather than i.i.d.) bootstrap is designed to preserve.

3. Price European options and compare to Black-Scholes

For the Black-Scholes comparison we use the pool’s own trailing annualized volatility (pool$ann_vol) as the input $\sigma$ – i.e. we ask: relative to the simplest possible historical-volatility GBM price, what does preserving the higher-order dependence structure (fat tails, volatility clustering) change?

bs_call <- function(S, K, r, sigma, T) {
  d1 <- (log(S / K) + (r + 0.5 * sigma^2) * T) / (sigma * sqrt(T))
  S * pnorm(d1) - K * exp(-r * T) * pnorm(d1 - sigma * sqrt(T))
}
bs_put <- function(S, K, r, sigma, T) {
  bs_call(S, K, r, sigma, T) - S + K * exp(-r * T)
}

T_grid <- c(0.0833, 0.25, 0.5, 1)         # 1m, 3m, 6m, 1y
moneyness <- seq(0.85, 1.15, by = 0.025)
K_grid <- round(S0 * moneyness)

results <- do.call(rbind, lapply(T_grid, function(Tj) {
  call_res <- eq_price(pool, S0 = S0, T = Tj, K = K_grid, type = "call")
  put_res  <- eq_price(pool, S0 = S0, T = Tj, K = K_grid, type = "put")
  bs_c <- bs_call(S0, K_grid, r_flat, pool$ann_vol, Tj)
  bs_p <- bs_put(S0, K_grid, r_flat, pool$ann_vol, Tj)

  data.frame(
    T = Tj, K = K_grid, moneyness = K_grid / S0,
    model_call = call_res$prices, bs_call = bs_c,
    model_put = put_res$prices, bs_put = bs_p,
    call_diff_pctS0 = 100 * (call_res$prices - bs_c) / S0
  )
}))

head(results, 10)
#>          T   K moneyness  model_call     bs_call   model_put       bs_put
#> 651 0.0833 651 0.8503024 117.0589105 117.0463048  0.01321958 6.139657e-04
#> 670 0.0833 670 0.8751192  98.1839134  98.1248425  0.06713435 8.063460e-03
#> 689 0.0833 689 0.8999360  79.3873002  79.2565216  0.19943301 6.865434e-02
#> 708 0.0833 708 0.9247528  60.8845909  60.6557233  0.62563548 3.967678e-01
#> 727 0.0833 727 0.9495696  43.0248635  42.9582520  1.69481986 1.628208e+00
#> 746 0.0833 746 0.9743865  26.6858217  27.3645694  4.28468992 4.963438e+00
#> 766 0.0833 766 1.0005094  12.8015368  14.7122544 10.32557536 1.223629e+01
#> 785 0.0833 785 1.0253262   4.8811160   6.9308703 21.33406635 2.338382e+01
#> 804 0.0833 804 1.0501430   1.7595853   2.7400882 37.14144746 3.812195e+01
#> 823 0.0833 823 1.0749599   0.6295614   0.9013219 54.94033538 5.521210e+01
#>     call_diff_pctS0
#> 651     0.001646480
#> 670     0.007715533
#> 689     0.017081631
#> 708     0.029893503
#> 727     0.008700447
#> 746    -0.088654506
#> 766    -0.249568015
#> 785    -0.267728264
#> 804    -0.128068195
#> 823    -0.035495948

4. Visualize

par(mfrow = c(2, 2), mar = c(4, 4, 2, 1))
for (Tj in T_grid) {
  d <- results[results$T == Tj, ]
  plot(d$K, d$bs_call, type = "b", pch = 1, lty = 2, col = "steelblue",
       xlab = "Strike", ylab = "Call price",
       main = sprintf("T = %.2f", Tj),
       ylim = range(c(d$bs_call, d$model_call)))
  lines(d$K, d$model_call, type = "b", pch = 19, col = "coral")
  abline(v = S0, col = "gray70", lty = 3)
}

Model call price (solid) vs Black-Scholes at historical volatility (dashed), by maturity.

par(mfrow = c(1, 1))

plot(results$moneyness, results$call_diff_pctS0, pch = 19,
     col = adjustcolor(as.numeric(as.factor(results$T)), 0.6),
     xlab = "Moneyness K/S0", ylab = "Model - BS (% of S0)",
     main = "Call price difference by moneyness")
abline(h = 0, col = "gray50", lty = 2)
legend("topleft", legend = sprintf("T=%.2f", sort(unique(results$T))),
       col = seq_along(unique(results$T)), pch = 19, cex = 0.8, bty = "n")

Price difference (model minus Black-Scholes), as a percentage of S0, by moneyness.

As in the DAX example in the companion paper, expect the largest deviations away from at-the-money strikes: the block bootstrap carries whatever skew and excess kurtosis are actually present in SPY’s realized residuals (see pool$skew and pool$excess_kurt below), while the Black-Scholes reference curve is, by construction, symmetric in log-moneyness.

data.frame(
  ann_vol = pool$ann_vol,
  skew = pool$skew,
  excess_kurt = pool$excess_kurt,
  phi = pool$phi
)
#>     ann_vol      skew excess_kurt         phi
#> 1 0.1526977 0.5427204    18.01553 -0.06417146

5. Internal consistency: put-call parity

Because calls and puts are priced from the same simulated terminal draws, put-call parity should hold to floating-point precision, not just Monte Carlo precision – this is a useful sanity check whenever you point empiricaloptions at a new underlying.

parity_error <- with(results, (model_call - model_put) - (S0 - K * exp(-r_flat * T)))
cat("Max abs parity error (price units):", max(abs(parity_error)), "\n")
#> Max abs parity error (price units): 1.136868e-13
cat("As % of S0:", 100 * max(abs(parity_error)) / S0, "\n")
#> As % of S0: 1.484918e-14

6. A path-dependent example: arithmetic Asian calls

The same pool prices Asian options with no extra calibration step – useful precisely when no option market (vanilla or exotic) exists to calibrate a path-dependent model against.

asian_prices <- eq_price_asian(pool, S0 = S0, T = 0.5, K = K_grid, n_fix = 6)
asian_prices
#>         651         670         689         708         727         746 
#> 122.0958688 103.7061590  85.5298652  67.8150142  51.0122775  35.7802068 
#>         766         785         804         823         842         861 
#>  22.4071803  13.0752898   7.0727351   3.5940242   1.7529912   0.8143500 
#>         880 
#>   0.3600664

Summary

Given nothing but SPY’s own trailing price history and a discount rate, empiricaloptions produced a full European surface across four maturities, an arithmetic Asian quote, and an internal put-call parity check – all without touching a single quoted option price. Swap "SPY" for any quantmod-gettable ticker, or replace the getSymbols() call with your own CSV of prices for underlyings quantmod cannot reach (a private index, a fund NAV series, an insurance-linked index), and the rest of this vignette runs unchanged.