In quantitative finance, pricing derivatives often requires working under a risk-neutral measure (Q) rather than the real-world physical measure (P). While the Girsanov theorem provides a theoretical framework for this change of measure, practical implementation can be challenging—especially for complex models like stochastic volatility with jumps.

This blog post demonstrates a semi-parametric approach that bridges classical time series modeling with risk-neutral pricing. We’ll show how to:

1. Simulate stock price paths under three different models (GBM, SVJD, Heston) using the physical measure
2. Extract the risk premium using ARIMA modeling
3. Transform physical paths to risk-neutral paths through residual resampling
4. Price European options under the risk-neutral measure

Models Compared

• Geometric Brownian Motion (GBM): The classic Black-Scholes model with constant volatility
• Stochastic Volatility Jump Diffusion (SVJD): Incorporates both stochastic volatility and price jumps
• Heston Model: Stochastic volatility without jumps (a special case of SVJD)

Methodological Approach

Our semi-parametric method involves:

1. Physical Simulation: Generate paths under the real-world measure with expected return μ
2. Risk Premium Extraction: Fit ARIMA models to discounted price increments to capture serial dependence
3. Residual Resampling: Use Gaussian density estimation to resample centered ARIMA residuals
4. Risk-Neutral Path Generation: Combine fitted ARIMA models with resampled residuals to create martingale paths
5. Option Pricing: Compute option prices as discounted expected payoffs under Q

R Packages Used

• esgtoolkit: For financial simulations and risk-neutral transformations
• ahead: For time series forecasting and residual resampling
• forecast (auto.arima): For automatic ARIMA model selection

The complete reproducible code is presented below, organized in logical sections from simulation to option pricing.

## ----setup, include=FALSE------------------------------------------------------
knitr::opts_chunk$set(echo = TRUE)


## ----1-simulate-SVJD, cache=TRUE-----------------------------------------------
# ARIMA-Black-Scholes: Semi-Parametric Risk-Neutral Pricing
# T. Moudiki
# 2025-12-04

library(esgtoolkit)
library(forecast)
library(ahead)

# =============================================================================
# 1 - SVJD SIMULATION (Physical Measure)
# =============================================================================

set.seed(123)
n <- 250L
h <- 5
freq <- "daily"
r <- 0.05
maturity <- 5
S0 <- 100
mu <- 0.08
sigma <- 0.04

# Simulate under physical measure with stochastic volatility and jumps
sim_GBM <- esgtoolkit::simdiff(n=n, horizon = h, frequency = freq, x0=S0, theta1 = mu, theta2 = sigma)
sim_SVJD <- esgtoolkit::rsvjd(n=n, r0=mu)
sim_Heston <- esgtoolkit::rsvjd(n=n, r0=mu, 
                                lambda = 0,
                                mu_J = 0,
                                sigma_J = 0)

cat("Simulation dimensions:\n")
cat("Start:", start(sim_SVJD), "\n")
cat("End:", end(sim_SVJD), "\n")
cat("Paths:", ncol(sim_SVJD), "\n")
cat("Time steps:", nrow(sim_SVJD), "\n\n")

par(mfrow=c(1, 3))
# Plot historical (physical measure) paths
esgtoolkit::esgplotbands(sim_GBM, main="GBM Paths under the Physical Measure", 
                         xlab="Time",
                         ylab="Stock prices")
esgtoolkit::esgplotbands(sim_SVJD, main="SVJD Paths under the Physical Measure", 
                         xlab="Time",
                         ylab="Stock prices")
esgtoolkit::esgplotbands(sim_Heston, main="Heston Paths under the Physical Measure", 
                         xlab="Time",
                         ylab="Stock prices")                         

# Summary statistics
cat("Physical measure statistics (GBM):\n")
terminal_prices_physical_GBM <- sim_GBM[nrow(sim_GBM), ]
cat("Mean terminal price:", mean(terminal_prices_physical_GBM), "\n")
cat("Std terminal price:", sd(terminal_prices_physical_GBM), "\n")
cat("Expected under P:", S0 * exp(mu * maturity), "\n\n")
 
cat("Physical measure statistics (SVJD):\n")
terminal_prices_physical_SVJD <- sim_SVJD[nrow(sim_SVJD), ]
cat("Mean terminal price:", mean(terminal_prices_physical_SVJD), "\n")
cat("Std terminal price:", sd(terminal_prices_physical_SVJD), "\n")
cat("Expected under P:", S0 * exp(mu * maturity), "\n\n")

cat("Physical measure statistics (Heston):\n")
terminal_prices_physical_Heston <- sim_Heston[nrow(sim_Heston), ]
cat("Mean terminal price:", mean(terminal_prices_physical_Heston), "\n")
cat("Std terminal price:", sd(terminal_prices_physical_Heston), "\n")
cat("Expected under P:", S0 * exp(mu * maturity), "\n\n")


## ----2-compute-discounted, cache=TRUE------------------------------------------
# =============================================================================
# 2 - COMPUTE DISCOUNTED PRICES (Transform to Martingale Domain)
# =============================================================================

discounted_prices_GBM <- esgtoolkit::esgdiscountfactor(r=r, X=sim_GBM)
discounted_prices_SVJD <- esgtoolkit::esgdiscountfactor(r=r, X=sim_SVJD)
discounted_prices_Heston <- esgtoolkit::esgdiscountfactor(r=r, X=sim_Heston)
martingale_diff_GBM <- discounted_prices_GBM - S0 
martingale_diff_SVJD <- discounted_prices_SVJD - S0 
martingale_diff_Heston <- discounted_prices_Heston - S0 

cat("Martingale differences dimensions (GBM):", dim(martingale_diff_GBM), "\n")
cat("Mean martingale diff (should be ≠ 0 under P):\n")
print(t.test(rowMeans(martingale_diff_GBM)))

cat("\nMartingale differences dimensions (SVJD):", dim(martingale_diff_SVJD), "\n")
cat("Mean martingale diff (should be ≠ 0 under P):\n")
print(t.test(rowMeans(martingale_diff_SVJD)))

cat("\nMartingale differences dimensions (Heston):", dim(martingale_diff_Heston), "\n")
cat("Mean martingale diff (should be ≠ 0 under P):\n")
print(t.test(rowMeans(martingale_diff_Heston)))

# =============================================================================
# 3 - VISUALIZE RISK PREMIUM
# =============================================================================

par(mfrow=c(2,2))

matplot(discounted_prices_GBM, type='l', col=rgb(0,0,1,0.3),
        main="Discounted Stock Prices (Physical Measure - GBM)",
        ylab="exp(-rt) * S_t", xlab="Time step")
abline(h=S0, col='red', lwd=2, lty=2)

matplot(discounted_prices_SVJD, type='l', col=rgb(0,0,1,0.3),
        main="Discounted Stock Prices (Physical Measure - SVJD)",
        ylab="exp(-rt) * S_t", xlab="Time step")
abline(h=S0, col='red', lwd=2, lty=2)

matplot(discounted_prices_Heston, type='l', col=rgb(0,0,1,0.3),
        main="Discounted Stock Prices (Physical Measure - Heston)",
        ylab="exp(-rt) * S_t", xlab="Time step")
abline(h=S0, col='red', lwd=2, lty=2)

par(mfrow=c(1,1))

mean_disc_path_GBM <- rowMeans(discounted_prices_GBM)
times_plot <- as.numeric(time(discounted_prices_GBM)) 
plot(times_plot, mean_disc_path_GBM, type='l', lwd=2, col='blue',
     main="Risk Premium in Discounted Prices (GBM)",
     xlab="Time (years)", ylab="E[exp(-rt)*S_t]")
abline(h=S0, col='red', lwd=2, lty=2)
lines(times_plot, S0 * exp((mu-r)*times_plot), col='green', lwd=2, lty=3)
legend("topleft",
       legend=c("Empirical mean", "S0", "Theoretical (μ-r drift)"),
       col=c('blue','red','green'), lty=c(1,2,3), lwd=2)

mean_disc_path_SVJD <- rowMeans(discounted_prices_SVJD)
times_plot <- as.numeric(time(discounted_prices_SVJD)) 
plot(times_plot, mean_disc_path_SVJD, type='l', lwd=2, col='blue',
     main="Risk Premium in Discounted Prices (SVJD)",
     xlab="Time (years)", ylab="E[exp(-rt)*S_t]")
abline(h=S0, col='red', lwd=2, lty=2)
lines(times_plot, S0 * exp((mu-r)*times_plot), col='green', lwd=2, lty=3)
legend("topleft",
       legend=c("Empirical mean", "S0", "Theoretical (μ-r drift)"),
       col=c('blue','red','green'), lty=c(1,2,3), lwd=2)

mean_disc_path_Heston <- rowMeans(discounted_prices_Heston)
times_plot <- as.numeric(time(discounted_prices_Heston)) 
plot(times_plot, mean_disc_path_Heston, type='l', lwd=2, col='blue',
     main="Risk Premium in Discounted Prices (Heston)",
     xlab="Time (years)", ylab="E[exp(-rt)*S_t]")
abline(h=S0, col='red', lwd=2, lty=2)
lines(times_plot, S0 * exp((mu-r)*times_plot), col='green', lwd=2, lty=3)
legend("topleft",
       legend=c("Empirical mean", "S0", "Theoretical (μ-r drift)"),
       col=c('blue','red','green'), lty=c(1,2,3), lwd=2)     


## ----3-fit-ARIMA, cache=TRUE, eval=TRUE----------------------------------------
# =============================================================================
# 4 - FIT ARIMA MODELS TO EXTRACT RISK PREMIUM
# =============================================================================

n_periods <- nrow(martingale_diff_GBM)
n_paths <- ncol(martingale_diff_GBM)

martingale_increments_GBM <- diff(martingale_diff_GBM)
martingale_increments_SVJD <- diff(martingale_diff_SVJD)
martingale_increments_Heston <- diff(martingale_diff_Heston)

# Initialize storage arrays
arima_residuals_GBM <- array(NA, dim = c(nrow(martingale_increments_GBM), n_paths))
centered_arima_residuals_GBM <- array(NA, dim = c(nrow(martingale_increments_GBM), n_paths))
means_arima_residuals_GBM <- rep(NA, n_paths)
arima_models_GBM <- list()

arima_residuals_SVJD <- array(NA, dim = c(nrow(martingale_increments_SVJD), n_paths))
centered_arima_residuals_SVJD <- array(NA, dim = c(nrow(martingale_increments_SVJD), n_paths))
means_arima_residuals_SVJD <- rep(NA, n_paths)
arima_models_SVJD <- list()

arima_residuals_Heston <- array(NA, dim = c(nrow(martingale_increments_Heston), n_paths))
centered_arima_residuals_Heston <- array(NA, dim = c(nrow(martingale_increments_Heston), n_paths))
means_arima_residuals_Heston <- rep(NA, n_paths)
arima_models_Heston <- list()

# Fit ARIMA models to GBM
cat("Fitting ARIMA models to", n_paths, "GBM paths...\n")
for (i in 1:n_paths) {
  y <- as.numeric(martingale_increments_GBM[, i])
  fit <- forecast::auto.arima(y, allowmean = FALSE)
  arima_models_GBM[[i]] <- fit
  
  res <- as.numeric(residuals(fit))
  arima_residuals_GBM[, i] <- res
  
  centre_arima_residuals <- scale(res, center = TRUE, scale = FALSE)
  means_arima_residuals_GBM[i] <- attr(centre_arima_residuals, "scaled:center")
  centered_arima_residuals_GBM[, i] <- centre_arima_residuals[,1]
}

# Fit ARIMA models to SVJD
cat("Fitting ARIMA models to", n_paths, "SVJD paths...\n")
for (i in 1:n_paths) {
  y <- as.numeric(martingale_increments_SVJD[, i])
  fit <- forecast::auto.arima(y, allowmean = FALSE)
  arima_models_SVJD[[i]] <- fit
  
  res <- as.numeric(residuals(fit))
  arima_residuals_SVJD[, i] <- res
  
  centre_arima_residuals <- scale(res, center = TRUE, scale = FALSE)
  means_arima_residuals_SVJD[i] <- attr(centre_arima_residuals, "scaled:center")
  centered_arima_residuals_SVJD[, i] <- centre_arima_residuals[,1]
}

# Fit ARIMA models to Heston
cat("Fitting ARIMA models to", n_paths, "Heston paths...\n")
for (i in 1:n_paths) {
  y <- as.numeric(martingale_increments_Heston[, i])
  fit <- forecast::auto.arima(y, allowmean = FALSE)
  arima_models_Heston[[i]] <- fit
  
  res <- as.numeric(residuals(fit))
  arima_residuals_Heston[, i] <- res
  
  centre_arima_residuals <- scale(res, center = TRUE, scale = FALSE)
  means_arima_residuals_Heston[i] <- attr(centre_arima_residuals, "scaled:center")
  centered_arima_residuals_Heston[, i] <- centre_arima_residuals[,1]
}

cat("\nARIMA model summary (first 5 GBM paths):\n")
for (i in 1:min(5, n_paths)) {
  cat("Path", i, ":", as.character(arima_models_GBM[[i]]), "\n")
}

# Box-Ljung tests
pvalues_GBM <- sapply(1:ncol(centered_arima_residuals_GBM), 
                  function(i) Box.test(centered_arima_residuals_GBM[,i])$p.value)
cat("\nBox-Ljung test p-values (GBM):\n")
cat("Mean p-value:", mean(pvalues_GBM), "\n")
cat("Proportion > 0.05:", mean(pvalues_GBM > 0.05), "\n")

pvalues_SVJD <- sapply(1:ncol(centered_arima_residuals_SVJD), 
                  function(i) Box.test(centered_arima_residuals_SVJD[,i])$p.value)
cat("\nBox-Ljung test p-values (SVJD):\n")
cat("Mean p-value:", mean(pvalues_SVJD), "\n")
cat("Proportion > 0.05:", mean(pvalues_SVJD > 0.05), "\n")

pvalues_Heston <- sapply(1:ncol(centered_arima_residuals_Heston), 
                  function(i) Box.test(centered_arima_residuals_Heston[,i])$p.value)
cat("\nBox-Ljung test p-values (Heston):\n")
cat("Mean p-value:", mean(pvalues_Heston), "\n")
cat("Proportion > 0.05:", mean(pvalues_Heston > 0.05), "\n")

par(mfrow=c(1,3))
hist(pvalues_GBM, breaks=20, col='lightgreen',
     main="Box-Ljung P-values (GBM)",
     xlab="P-value")
abline(v=0.05, col='red', lwd=2, lty=2)

hist(pvalues_SVJD, breaks=20, col='lightblue',
     main="Box-Ljung P-values (SVJD)",
     xlab="P-value")
abline(v=0.05, col='red', lwd=2, lty=2)

hist(pvalues_Heston, breaks=20, col='lightcoral',
     main="Box-Ljung P-values (Heston)",
     xlab="P-value")
abline(v=0.05, col='red', lwd=2, lty=2)
par(mfrow=c(1,1))


## ----4-generate-rn-paths, cache=TRUE, eval=TRUE--------------------------------
# =============================================================================
# 5 - GENERATE RISK-NEUTRAL PATHS
# =============================================================================

cat("\n\nGenerating risk-neutral paths from ALL historical paths...\n")
n_sim_per_path <- 20  # Generate 10 paths per historical path
times <- seq(0, maturity, length.out = n_periods)
discount_factor <- exp(r * times)

# Storage for all risk-neutral paths
all_S_tilde_GBM <- list()
all_S_tilde_SVJD <- list()
all_S_tilde_Heston <- list()

# Generate GBM risk-neutral paths
for (i in 1:n_paths) {
  resampled_residuals <- ahead::rgaussiandens(centered_arima_residuals_GBM[, i], 
                                              p = n_sim_per_path)
  fit <- arima_models_GBM[[i]]
  fitted_increments <- as.numeric(fitted(fit))
  
  discounted_path <- matrix(0, nrow = n_periods, ncol = n_sim_per_path)
  discounted_path[1, ] <- S0
  
  increments <- matrix(scale(fitted_increments, center = TRUE, 
                             scale = FALSE)[,1], nrow = n_periods - 1, ncol = n_sim_per_path) + 
                resampled_residuals[1:(n_periods - 1), ]
  
  discounted_path[-1, ] <- S0 + apply(increments, 2, cumsum)
  S_tilde_price <- discounted_path * discount_factor
  all_S_tilde_GBM[[i]] <- S_tilde_price
}

# Generate SVJD risk-neutral paths
for (i in 1:n_paths) {
  resampled_residuals <- ahead::rgaussiandens(centered_arima_residuals_SVJD[, i], 
                                              p = n_sim_per_path)
  fit <- arima_models_SVJD[[i]]
  fitted_increments <- as.numeric(fitted(fit))
  
  discounted_path <- matrix(0, nrow = n_periods, ncol = n_sim_per_path)
  discounted_path[1, ] <- S0
  
  increments <- matrix(scale(fitted_increments, center = TRUE, 
                             scale = FALSE)[,1], nrow = n_periods - 1, ncol = n_sim_per_path) + 
                resampled_residuals[1:(n_periods - 1), ]
  
  discounted_path[-1, ] <- S0 + apply(increments, 2, cumsum)
  S_tilde_price <- discounted_path * discount_factor
  all_S_tilde_SVJD[[i]] <- S_tilde_price
}

# Generate Heston risk-neutral paths
for (i in 1:n_paths) {
  resampled_residuals <- ahead::rgaussiandens(centered_arima_residuals_Heston[, i], 
                                              p = n_sim_per_path)
  fit <- arima_models_Heston[[i]]
  fitted_increments <- as.numeric(fitted(fit))
  
  discounted_path <- matrix(0, nrow = n_periods, ncol = n_sim_per_path)
  discounted_path[1, ] <- S0
  
  increments <- matrix(scale(fitted_increments, center = TRUE, 
                             scale = FALSE)[,1], nrow = n_periods - 1, ncol = n_sim_per_path) + 
                resampled_residuals[1:(n_periods - 1), ]
  
  discounted_path[-1, ] <- S0 + apply(increments, 2, cumsum)
  S_tilde_price <- discounted_path * discount_factor
  all_S_tilde_Heston[[i]] <- S_tilde_price
}

# Combine all paths
S_tilde_combined_GBM <- do.call(cbind, all_S_tilde_GBM)
cat("Total GBM risk-neutral paths generated:", ncol(S_tilde_combined_GBM), "\n")

S_tilde_combined_SVJD <- do.call(cbind, all_S_tilde_SVJD)
cat("Total SVJD risk-neutral paths generated:", ncol(S_tilde_combined_SVJD), "\n")

S_tilde_combined_Heston <- do.call(cbind, all_S_tilde_Heston)
cat("Total Heston risk-neutral paths generated:", ncol(S_tilde_combined_Heston), "\n\n")

# Convert to time series
S_tilde_ts_GBM <- ts(S_tilde_combined_GBM, start = start(sim_GBM), 
                     frequency = frequency(sim_GBM))
S_tilde_ts_SVJD <- ts(S_tilde_combined_SVJD, start = start(sim_SVJD), 
                      frequency = frequency(sim_SVJD))
S_tilde_ts_Heston <- ts(S_tilde_combined_Heston, start = start(sim_Heston), 
                        frequency = frequency(sim_Heston))

# Visualize risk-neutral paths
par(mfrow=c(1,3))
esgtoolkit::esgplotbands(S_tilde_ts_GBM, 
                        main="Risk-Neutral Paths - GBM")
esgtoolkit::esgplotbands(S_tilde_ts_SVJD, 
                        main="Risk-Neutral Paths - SVJD")
esgtoolkit::esgplotbands(S_tilde_ts_Heston, 
                        main="Risk-Neutral Paths - Heston")

# Sample plots
par(mfrow=c(1,3))
matplot(S_tilde_combined_GBM[, sample(ncol(S_tilde_combined_GBM), 200)], 
        type='l', col=rgb(0,0,1,0.1),
        main="Risk-Neutral Paths (GBM)",
        xlab="Time step", ylab="Stock Price")
lines(rowMeans(S_tilde_combined_GBM), col='red', lwd=3)
abline(h=S0, col='green', lwd=2, lty=2)

matplot(S_tilde_combined_SVJD[, sample(ncol(S_tilde_combined_SVJD), 200)], 
        type='l', col=rgb(0,0,1,0.1),
        main="Risk-Neutral Paths (SVJD)",
        xlab="Time step", ylab="Stock Price")
lines(rowMeans(S_tilde_combined_SVJD), col='red', lwd=3)
abline(h=S0, col='green', lwd=2, lty=2)

matplot(S_tilde_combined_Heston[, sample(ncol(S_tilde_combined_Heston), 200)], 
        type='l', col=rgb(0,0,1,0.1),
        main="Risk-Neutral Paths (Heston)",
        xlab="Time step", ylab="Stock Price")
lines(rowMeans(S_tilde_combined_Heston), col='red', lwd=3)
abline(h=S0, col='green', lwd=2, lty=2)
par(mfrow=c(1,1))


## ----5-rn-verif, cache=TRUE, eval=TRUE-----------------------------------------
# =============================================================================
# 6 - VERIFY RISK-NEUTRAL PROPERTY
# =============================================================================

cat("\n=== RISK-NEUTRAL VERIFICATION ===\n\n")

terminal_prices_rn_GBM <- S_tilde_combined_GBM[n_periods, ]
terminal_prices_rn_SVJD <- S_tilde_combined_SVJD[n_periods, ]
terminal_prices_rn_Heston <- S_tilde_combined_Heston[n_periods, ]
capitalized_stock_price <- S0 * exp(r * maturity)

cat("GBM Risk-Neutral Verification:\n")
cat("Expected terminal price (Q):", capitalized_stock_price, "\n")
cat("Empirical mean:", mean(terminal_prices_rn_GBM), "\n")
cat("Difference:", mean(terminal_prices_rn_GBM) - capitalized_stock_price, "\n")
print(t.test(terminal_prices_rn_GBM - capitalized_stock_price))

cat("\nSVJD Risk-Neutral Verification:\n")
cat("Expected terminal price (Q):", capitalized_stock_price, "\n")
cat("Empirical mean:", mean(terminal_prices_rn_SVJD), "\n")
cat("Difference:", mean(terminal_prices_rn_SVJD) - capitalized_stock_price, "\n")
print(t.test(terminal_prices_rn_SVJD - capitalized_stock_price))

cat("\nHeston Risk-Neutral Verification:\n")
cat("Expected terminal price (Q):", capitalized_stock_price, "\n")
cat("Empirical mean:", mean(terminal_prices_rn_Heston), "\n")
cat("Difference:", mean(terminal_prices_rn_Heston) - capitalized_stock_price, "\n")
print(t.test(terminal_prices_rn_Heston - capitalized_stock_price))

# Visualization comparison
par(mfrow=c(3, 2))
hist(terminal_prices_physical_GBM, breaks=30, col=rgb(1,0,0,0.5),
     main="Terminal Prices: Physical (GBM)",
     xlab="Price", xlim=c(50, 300))
abline(v=mean(terminal_prices_physical_GBM), col='red', lwd=2)
abline(v=S0*exp(mu*maturity), col='blue', lwd=2, lty=2)

hist(terminal_prices_rn_GBM, breaks=30, col=rgb(0,0,1,0.5),
     main="Terminal Prices: Risk-Neutral (GBM)",
     xlab="Price", xlim=c(50, 300))
abline(v=mean(terminal_prices_rn_GBM), col='blue', lwd=2)
abline(v=S0*exp(r*maturity), col='red', lwd=2, lty=2)

hist(terminal_prices_physical_SVJD, breaks=30, col=rgb(1,0,0,0.5),
     main="Terminal Prices: Physical (SVJD)",
     xlab="Price", xlim=c(50, 300))
abline(v=mean(terminal_prices_physical_SVJD), col='red', lwd=2)
abline(v=S0*exp(mu*maturity), col='blue', lwd=2, lty=2)

hist(terminal_prices_rn_SVJD, breaks=30, col=rgb(0,0,1,0.5),
     main="Terminal Prices: Risk-Neutral (SVJD)",
     xlab="Price", xlim=c(50, 300))
abline(v=mean(terminal_prices_rn_SVJD), col='blue', lwd=2)
abline(v=S0*exp(r*maturity), col='red', lwd=2, lty=2)

hist(terminal_prices_physical_Heston, breaks=30, col=rgb(1,0,0,0.5),
     main="Terminal Prices: Physical (Heston)",
     xlab="Price", xlim=c(50, 300))
abline(v=mean(terminal_prices_physical_Heston), col='red', lwd=2)
abline(v=S0*exp(mu*maturity), col='blue', lwd=2, lty=2)

hist(terminal_prices_rn_Heston, breaks=30, col=rgb(0,0,1,0.5),
     main="Terminal Prices: Risk-Neutral (Heston)",
     xlab="Price", xlim=c(50, 300))
abline(v=mean(terminal_prices_rn_Heston), col='blue', lwd=2)
abline(v=S0*exp(r*maturity), col='red', lwd=2, lty=2)
par(mfrow=c(1,1))


## ----6-option-pricing, cache=TRUE, eval=TRUE-----------------------------------
# =============================================================================
# 7 - OPTION PRICING
# =============================================================================

cat("\n=== OPTION PRICING ===\n\n")

bs_price <- function(S, K, r, sigma, T, q = 0) {

  d1 <- (log(S / K) + (r - q + 0.5 * sigma^2) * T) / (sigma * sqrt(T))
  d2 <- d1 - sigma * sqrt(T)

  call <- S * exp(-q * T) * pnorm(d1) - K * exp(-r * T) * pnorm(d2)
  put  <- K * exp(-r * T) * pnorm(-d2) - S * exp(-q * T) * pnorm(-d1)

  list(call = call, put = put)
}


strikes <- seq(80, 160, by=10)
d_f <- exp(-r * maturity)

# Function to price options
price_options <- function(terminal_prices, strikes, discount_factor) {
  n_strikes <- length(strikes)
  call_prices <- numeric(n_strikes)
  bs_call_prices <- numeric(n_strikes)
  put_prices <- numeric(n_strikes)
  bs_put_prices <- numeric(n_strikes)
  call_se <- numeric(n_strikes)
  put_se <- numeric(n_strikes)
  
  for (k in 1:n_strikes) {
    K <- strikes[k]
    call_payoffs <- pmax(terminal_prices - K, 0)
    call_prices[k] <- mean(call_payoffs) * discount_factor
    call_se[k] <- sd(call_payoffs) / sqrt(length(call_payoffs)) * discount_factor
    
    put_payoffs <- pmax(K - terminal_prices, 0)
    put_prices[k] <- mean(put_payoffs) * discount_factor
    put_se[k] <- sd(put_payoffs) / sqrt(length(put_payoffs)) * discount_factor
    
    bs_prices <- bs_price(S0, K, r, sigma=sigma, T=5, q = 0)
    bs_call_prices[k] <- bs_prices$call
    bs_put_prices[k] <- bs_prices$put
  }
  
  list(call_prices = call_prices, put_prices = put_prices,
       bs_call_prices = bs_call_prices, bs_put_prices = bs_put_prices,
       call_se = call_se, put_se = put_se)
}

# Price options for all models
options_GBM <- price_options(terminal_prices_rn_GBM, strikes, d_f)
options_SVJD <- price_options(terminal_prices_rn_SVJD, strikes, d_f)
options_Heston <- price_options(terminal_prices_rn_Heston, strikes, d_f)


## ------------------------------------------------------------------------------
kableExtra::kable(as.data.frame(options_GBM))
kableExtra::kable(as.data.frame(options_Heston))
kableExtra::kable(as.data.frame(options_SVJD))

Full code and additional details are available on GitHub:

https://github.com/thierrymoudiki/2025-12-07-risk-neutralization-with-ARIMA

The preprint of the associated research paper can be found here:

https://www.researchgate.net/publication/398427354_An_ARIMA-Based_Semi-Parametric_Approach_to_Market_Price_of_Risk_estimation_and_Risk-Neutral_Pricing