Trying to use Black Scholes on market data estimating the volatility parameter

asian.jl

# Monte Carlo pricing for Asian options


function montecarlo_price_asian(payoff, S0, T, r, sigma, N)
    # Generate N random samples from a standard normal distribution
    mc_average = 0.0
    for i in 1:N
        Z = randn()
    
       
        # Simulate the stock price at maturity T
       
        S = simulate_geometric_brownian_motion_path(S0, r, sigma, T, 100)[2] 
        integralT  = sum(S) * T / length(S) / T # Average price over the path   
             
        # Discount the average payoff back to present value
        mc_average += exp(-r * T) * payoff(S[end], integralT)
    end
    return mc_average / N
end

mc_call_price_asian = montecarlo_price_asian((ST, integralT) -> max(integralT - K, 0), S0, T, r, sigma, N)

blackscholes.jl

using Distributions
using Downloads
using DelimitedFiles
using MarketData
using Dates


"""
    business_days(start_date, end_date)

Counts business days (Mon–Fri) between start_date (exclusive)
and end_date (inclusive).
"""
function business_days(start_date::Date, end_date::Date)
    start_date < end_date || return 0
    count = 0
    d = start_date + Day(1)
    while d <= end_date
        if dayofweek(d) ≤ 5   # 1=Mon, ..., 5=Fri
            count += 1
        end
        d += Day(1)
    end
    return count
end

"""
    load_yahoo_symbol(symbol; startdate=nothing, enddate=today())

Downloads daily data from Yahoo via MarketData.jl
and returns:

    ts :: Vector{Float64}   # time in years since first observation
    Ss :: Vector{Float64}   # closing prices

"""
function load_yahoo_symbol(symbol::String; startdate=nothing, enddate=today())

    # Download
    ta = startdate === nothing ?
        yahoo(Symbol(symbol), YahooOpt(period2 = DateTime(enddate))) :
        yahoo(Symbol(symbol), YahooOpt(period1 = DateTime(startdate)), period2 = DateTime(enddate))

    # Extract dates and closing prices
    dates = Date.(timestamp(ta))
    closes = values(ta)[:, findfirst(==(Symbol("Close")), colnames(ta))]

    Ss = Float64.(closes)

    return dates, Ss
end


d1(S, K, τ , r, σ) = (log(S / K) + (r + σ^2 / 2) * τ) / (σ * sqrt(τ))

d2(S, K, τ, r, σ) = d1(S, K, τ, r, σ) - σ * sqrt(τ)

function blackscholesprice(S, K, t, T, r, σ) 
    S * cdf(Normal(), d1(S, K, T - t, r, σ)) - K * exp(-r * (T - t)) * cdf(Normal(), d2(S, K, T - t, r, σ))
end

quadvar_logprice(S, σ, T) = σ^2 * T

function estimate_quadraticvariation(ts, Ss)
    logreturns = diff(log.(Ss))
    return sum(logreturns.^2)
end
    
function estimate_sigma(ts, Ss)
    T = ts[end] - ts[1]
    return sqrt(estimate_quadraticvariation(ts, Ss) / T)
end

function simulate_geometric_brownian_motion_path(S0, μ, σ, T, N)
    dt = T / N
    W = cumsum(sqrt(dt) * randn(N))
    t = 0:dt:T
    S = S0 * exp.((μ - σ^2 / 2) * t[1:end-1] + σ * W)
    return t[1:end-1], S
end 


S0 = 100.0
μ = 0.07    
σ = 0.2
T = 1.0
N = 1_000_000
r = 0.05

ts, Ss = simulate_geometric_brownian_motion_path(S0, μ, σ, T, N)
estimated_sigma = estimate_sigma(ts, Ss)
println("Estimated σ: ", estimated_sigma)


K = 100.0
# blackscholesprice(S, K, t, T, r, σ)
price = blackscholesprice(S0, K, 0.0, T, r, σ)
println("Black–Scholes price: ", price) 









# Downliad SPY
ds, Ss = load_yahoo_symbol("SPY"; enddate=Date("2026-02-11"))
Ss = Ss[end-100+1:end] # last 100 observations

ts = (0:length(Ss) .- 1)/252 # time in years since first observation, assuming 252 trading days per year

estimated_sigma = estimate_sigma(ts, Ss)

println("Estimated σ for SPY (annualized): ", estimated_sigma)


#= 
Real call price, looked up on February 10:
Instrument: SPY February 24, 2026 700 Call
Option Symbol: SPY260224C00700000

Underlying: SPY
Expiration Date: February 24, 2026
Strike Price: 700
Option Type: Call
Implied Volatility	11.65%

Last traded price: 3.46
=#

# Compare with Black–Scholes price for the same parameters:
S0 = Ss[end]
K = 700.0     
T = business_days(Date("2026-02-10"), Date("2026-02-24"))/252 # time to expiration in years (13 trading days) 
r = 0.0175      # risk-free rate 


price2 = blackscholesprice(S0, K, 0.0, T, r, estimated_sigma)

mcgreeks.jl

function ∂montecarlo_call_price(S0, T, r, sigma, N)
    mc_average = 0.0
    for i in 1:N
        Z = randn()
        ST = S0 * exp((r - 0.5 * sigma^2) * T + sigma * sqrt(T) * Z)
        if ST > K
            mc_average += exp(-r * T) * exp((r - 0.5 * sigma^2) * T + sigma * sqrt(T) * Z)
        else 
            mc_average += 0.0
        end
    end
    return mc_average / N
end
S0, K, T, r, sigma, N = 10., 12.0, 3, 0.05, 0.13, 1_000_000



using ForwardDiff



mc_greek = ForwardDiff.derivative(x -> montecarlo_price(call_payoff, x, T, r, sigma, N), S0)
mc_greek2 = ∂montecarlo_call_price(S0, T, r, sigma, N)

bs_greek = ForwardDiff.derivative(x -> blackscholesprice(x, K, 0.0, T, r, sigma), S0)
Δ = cdf(Normal(), d1(S0, K, T, r, sigma))

bs_greek - Δ

mc_greek - Δ

monte_carlo.jl

# Monte Carlo estimator Black Scholes Call price 
using ForwardDiff
using Distributions

d1(S, K, τ , r, σ) = (log(S / K) + (r + σ^2 / 2) * τ) / (σ * sqrt(τ))

d2(S, K, τ, r, σ) = d1(S, K, τ, r, σ) - σ * sqrt(τ)

function blackscholesprice(S, K, t, T, r, σ) 
    S * cdf(Normal(), d1(S, K, T - t, r, σ)) - K * exp(-r * (T - t)) * cdf(Normal(), d2(S, K, T - t, r, σ))
end

function simulate_geometric_brownian_motion_path(S0, μ, σ, T, N)
    dt = T / N
    W = cumsum(sqrt(dt) * randn(N))
    t = 0:dt:T
    S = S0 * exp.((μ - σ^2 / 2) * t[1:end-1] + σ * W)
    return t[1:end-1], S
end 

function simulate_sde(S0, drift, diffusion, T, N)
    dt = T / N
    # Euler-Maruyama method
    S = zeros(N)
    S[1] = S0
    for i in 2:N
        dW = sqrt(dt) * randn()
        S[i] = S[i-1] + drift(S[i-1]) * dt + diffusion(S[i-1]) * dW
    end
    return S    
end

function montecarlo_price(payoff, S0, T, r, sigma, N)
    # Generate N random samples from a standard normal distribution
    mc_average = 0.0
    for i in 1:N
        Z = randn()
    
       
        # Simulate the stock price at maturity T
        #ST = S0 * exp((r - 0.5 * sigma^2) * T + sigma * sqrt(T) * Z)
       
        ST = simulate_geometric_brownian_motion_path(S0, r, sigma, T, 100)[2][end] # Simulate the stock price at maturity T

             
        # Discount the average payoff back to present value
        mc_average += exp(-r * T) * payoff(ST)
    end
    return mc_average / N
end

S0, K, T, r, sigma, N = 10., 12.0, 3, 0.05, 0.13, 1_000_000

call_payoff(ST) = max(ST - K, 0)
put_payoff(ST) = max(K - ST, 0)
    

mc_call_price = montecarlo_price(call_payoff, S0, T, r, sigma, N)

bs_call_price = blackscholesprice(S0, K, 0.0, T, r, sigma)

println("Monte Carlo Call Price: $mc_call_price")
println("Black-Scholes Call Price: $bs_call_price")
println("Difference: $(abs(mc_call_price - bs_call_price))")


mc_put_price = montecarlo_price(put_payoff, S0, T, r, sigma, N)

# Check the put-call parity
parity = mc_call_price - mc_put_price - (S0 - K * exp(-r * T))

pde.jl

# Black scholes PDE solver using finite difference method
using LinearAlgebra, Distributions

# Parameters
S0 = 95.0
T = 1.0        # Time to maturity
K = 100.0      # Strike price
r = 0.05       # Risk-free interest rate
σ = 0.2    # Volatility 
g(S) = max(S - K, 0)  # Payoff function for a call option

# Closed form solution for Black–Scholes price
using Distributions
d1(S, K, τ , r, σ) = (log(S / K) + (r + σ^2 / 2) * τ) / (σ * sqrt(τ))
d2(S, K, τ, r, σ) = d1(S, K, τ, r, σ) - σ * sqrt(τ)
function blackscholesprice(S, K, t, T, r, σ) 
    S * cdf(Normal(), d1(S, K, T - t, r, σ)) - K * exp(-r * (T - t)) * cdf(Normal(), d2(S, K, T - t, r, σ))
end





# Discretization parameters
S_max = 200.0  # Maximum stock price / roof
S_min = 0.0    # Minimum stock price
ts = range(0, T, length=40_001) # Time grid)
xs = range(S_min, S_max, length=1_001)
M = length(xs)
dx = (S_max - S_min) / (M - 1)
N = length(ts)
dt = T / (N - 1)

# Stability check for explicit Euler (diffusion-dominated constraint)
a_max = 0.5 * σ^2 * S_max^2
dt_max = dx^2 / (2 * a_max)   # ≈ dx^2 / (σ^2 * S_max^2)

if dt > dt_max
    println("Warning: dt = $dt > $dt_max (stability limit). The scheme is likely unstable.")
else
    println("dt = $dt <= $dt_max (diffusion stability limit satisfied).")
end




function ∂(v) 
    dv = [v[i+1] - v[i-1] for i in 2:(length(v)-1)]/(2*dx)
    dv = [dv[1]; dv; dv[end]] # padding boundary values
    return dv
end 

function ∂²(v) 
    d²v = [v[i+1] - 2*v[i] + v[i-1] for i in 2:(length(v)-1)]/(dx^2)
    d²v = [d²v[1]; d²v; d²v[end]] # padding boundary values
    return d²v
end


unitvector(M, i) = [j == i ? 1.0 : 0.0 for j in 1:M]
function function_as_matrix(f, M)
    fvectors = [f(unitvector(M, i)) for i in 1:M] # finite difference matrix for first derivative
    fmatrix = [fvectors[i][j] for j in 1:M, i in 1:M] # convert to matrix form
    return fmatrix
end
∂matrix = function_as_matrix(∂, M)
∂²matrix = function_as_matrix(∂², M)





# Time stepping backwards
# ∂v/∂t + r*S*∂v/∂S + 0.5*σ^2*S^2*∂²v/∂S² - r*v = 0
# V[i, j] ≈ v(t_i, S_j)

# (V[i+1, :] - V[i, :])/dt + map(S->r*S, xs).*∂(V[i+1, :]) + 0.5*σ^2*map(S->S^2, xs).*∂²(V[i+1, :]) - r*V[i+1, :] = 0

V = zeros(N, M)
V[end, :] = g.(xs) # terminal condition at maturity
using SparseArrays

L = sparse(r*Diagonal(xs)*∂matrix + 0.5*σ^2*Diagonal(xs.^2)*∂²matrix - r*I)
A = I + dt*L
MODE = :crank_nicolson # choose between :explicit, :matrix_explicit, :crank_nicolson

for i in (N-1):-1:1

    if MODE == :explicit
        # Solve for V[i, :]
        V[i, :] = V[i+1, :] + dt*(r*xs.*∂(V[i+1, :]) + 0.5*σ^2*map(S->S^2, xs).*∂²(V[i+1, :]) - r*V[i+1, :])
    elseif MODE == :matrix_explicit
        V[i, :] = A*V[i+1, :]
    elseif MODE == :crank_nicolson
        V[i, :] = ((I - 0.5*dt*L) \ ((I + 0.5*dt*L) * V[i+1, :]))
    end

    V[i, 1] = exp(-r*(T - ts[i])) * g(0.0) 
    V[i, end] = S_max - K*exp(-r*(T - ts[i])) # boundary condition at S=S_max from formula script
    # if g is not a call, you could use the following instead:
    gprime = (g(S_max) - g(S_max - dx)) / dx
    V[i, end] = gprime*S_max*exp(-r*(T - ts[i])) + (g(S_max) - gprime*S_max)*exp(-r*(T - ts[i]))
end

price1 = blackscholesprice(S0, K, 0.0, T, r, σ)
price2 = V[1, findfirst(x -> x >= S0, xs)]
println("Black–Scholes price (closed form): ", price1)
println("Black–Scholes price (PDE solver): ", price2)