如何在R语言中估算Black-Scholes、GBM及布朗运动的参数（含漂移与波动率）

Hey there! Let's walk through exactly how to calculate the drift and volatility parameters for Geometric Brownian Motion (GBM)—which are the core inputs for the Black-Scholes model—using R. I'll cover both manual calculations (so you understand the math behind it) and convenient package-based methods to save time.

Quick Background #

GBM follows the stochastic differential equation:
$$dS_t = \mu S_t dt + \sigma S_t dW_t$$
When we discretize this for daily price data, the log returns ($r_t = \ln(S_t/S_{t-1})$) are normally distributed with mean $\mu - \sigma^2/2$ and variance $\sigma^2$. This is the key relationship we'll use to estimate our parameters.

Step 1: Prepare Your Data #

First, let's get some price data—either real historical stock prices or simulated GBM data for testing.

Real Stock Data (using quantmod) #

# Install/load required packages
if (!require("quantmod")) install.packages("quantmod")
library(quantmod)

# Pull adjusted closing prices for Apple (2020-2023)
getSymbols("AAPL", from = "2020-01-01", to = "2023-12-31")
prices <- Cl(AAPL)  # Extract closing prices

Simulated GBM Data (for testing) #

If you don't have real data, you can generate synthetic GBM prices to validate your estimates:

set.seed(123)  # For reproducibility
n_days <- 1000
true_mu <- 0.15    # True annual drift
true_sigma <- 0.2  # True annual volatility
S0 <- 100          # Initial price

# Generate GBM path
time <- seq(0, 1, length.out = n_days)
brownian <- c(0, cumsum(rnorm(n_days-1, 0, sqrt(1/(n_days-1)))))
sim_prices <- S0 * exp((true_mu - 0.5*true_sigma^2)*time + true_sigma*brownian)

# Convert to xts object (mimics real price data structure)
prices_sim <- xts(sim_prices, order.by = seq.Date(as.Date("2020-01-01"), by = "day", length.out = n_days))

Step 2: Calculate Log Returns #

Log returns are the foundation of our parameter estimates:

# For real data
log_returns <- diff(log(prices))[-1]  # Remove the first NA value

# For simulated data
log_returns_sim <- diff(log(prices_sim))[-1]

Step 3: Estimate Volatility ($\sigma$) #

Volatility is the standard deviation of log returns, scaled to an annualized value (assuming 252 trading days per year):

# Daily volatility
sigma_daily <- sd(log_returns)

# Annualized volatility
sigma_annual <- sigma_daily * sqrt(252)

cat("Annualized Volatility (Manual):", round(sigma_annual, 4), "\n")

Step 4: Estimate Drift ($\mu$) #

The drift term accounts for the expected return of the asset. Using the GBM relationship, we adjust the mean of log returns by half the variance of returns to get the true drift:

# Daily drift component
mu_daily <- mean(log_returns) + 0.5 * sigma_daily^2

# Annualized drift
mu_annual <- mu_daily * 252

cat("Annualized Drift (Manual):", round(mu_annual, 4), "\n")

Using Packages for Faster Estimation #

If you want to skip the manual math, the fOptions package has a built-in function to fit GBM parameters directly:

if (!require("fOptions")) install.packages("fOptions")
library(fOptions)

# Fit GBM to price data (dt = 1/252 for daily data)
gbm_fit <- GBMFit(x = as.vector(prices), dt = 1/252)

cat("Annualized Drift (fOptions):", round(gbm_fit@mu, 4), "\n")
cat("Annualized Volatility (fOptions):", round(gbm_fit@sigma, 4), "\n")

Key Notes #

• Adjusted Prices: Always use split/dividend-adjusted prices for real stock data—quantmod pulls adjusted prices by default, but double-check if you're using other data sources.
• Drift vs Risk-Free Rate: For Black-Scholes pricing, many practitioners use the risk-free rate instead of historical drift (since the model assumes no arbitrage). Use historical drift if you're simulating future price paths, not pricing options.
• Scaling: If you're using intraday data (e.g., hourly), adjust the scaling factor: use sqrt(252*24) for hourly data instead of sqrt(252).

内容的提问来源于stack exchange，提问作者Mahdi Moghimi