Martingales by simulation (Julia)

martingales.jl

###############################################################################
# Martingales by simulation (Julia)
#
# This script illustrates several facts discussed in lecture:
#
# (1) Brownian motion W_t has independent increments and is a martingale.
# (2) Additive Brownian stock model S_t = S_0 + μ t + σ W_t.
# (3) The conditional expectation (filtered value)
#         V_t = E[(S_T - K)^+ | 𝔽_t]
#     is a martingale in t (tower property).
#     We illustrate this numerically by a binning regression:
#         E[V_t | S_s] ≈ V_s  (since V_s is a function of S_s).
# (4) Discrete stochastic integral / innovation sum:
#         M_n = Σ X_k (W_{t_{k+1}} - W_{t_k})
#     with predictable X_k is a martingale (mean 0, no drift).
#
# Requirements:
#   Julia ≥ 1.8 recommended
#   Packages: Random, Statistics, Distributions
###############################################################################

using Random
using Statistics
using Distributions

Random.seed!(1)  # reproducibility

###############################################################################
# Part A. Simulate Brownian motion and additive Brownian stock
###############################################################################

"""
    simulate_brownian(T, N; rng=Random.GLOBAL_RNG)

Simulate standard Brownian motion on [0,T] on an equidistant grid with N steps.
Returns (tgrid, W) where:
  - tgrid has length N+1 (includes 0 and T),
  - W has length N+1 with W[1]=0.
"""
function simulate_brownian(T::Float64, N::Int; rng=Random.GLOBAL_RNG)
    dt = T / N
    tgrid = collect(0.0:dt:T)
    W = zeros(N + 1)
    for n in 2:(N + 1)
        W[n] = W[n-1] + sqrt(dt) * randn(rng)
    end
    return tgrid, W
end

"""
    additive_stock_from_W(S0, μ, σ, tgrid, W)

Build S_t = S0 + μ t + σ W_t from an already simulated Brownian path.
"""
function additive_stock_from_W(S0, μ, σ, tgrid, W)
    @assert length(tgrid) == length(W)
    S = similar(W)
    for i in eachindex(W)
        S[i] = S0 + μ * tgrid[i] + σ * W[i]
    end
    return S
end

# Parameters
T  = 1.0
N  = 1000
S0 = 0.0
μ  = 0.5
σ  = 2.0

tgrid, W = simulate_brownian(T, N)
S = additive_stock_from_W(S0, μ, σ, tgrid, W)

println("One path: W(T) = $(W[end]), S(T) = $(S[end])")

###############################################################################
# Part B. Closed-form conditional expectation for a call in additive model
#
# Model: S_T = S_t + μ (T-t) + σ sqrt(T-t) Z,  Z ~ N(0,1) independent of 𝔽_t.
# Payoff: H = (S_T - K)^+.
#
# If Y ~ N(m, v^2), then E[(Y-K)^+] = (m-K) Φ(d) + v φ(d),  d=(m-K)/v.
###############################################################################

"""
    call_condexp_additive(S_t, t, T, K, μ, σ)

Closed-form value of V_t = E[(S_T - K)^+ | 𝔽_t] in the additive Brownian model,
given the current stock value S_t.

Returns V_t (a number).
"""
function call_condexp_additive(S_t::Float64, t::Float64, T::Float64,
                               K::Float64, μ::Float64, σ::Float64)
    τ = T - t
    if τ <= 0
        return max(S_t - K, 0.0)
    end
    m = S_t + μ * τ
    v = σ * sqrt(τ)
    d = (m - K) / v
    Φ = cdf(Normal(0,1), d)
    φ = pdf(Normal(0,1), d)
    return (m - K) * Φ + v * φ
end

###############################################################################
# Part C. Monte Carlo: conditional expectation by simulating the future
#
# Given a realized S_t, we can approximate V_t by simulating the conditional
# future increments from t to T:
#    S_T = S_t + μ (T-t) + σ (W_T - W_t)
# where (W_T - W_t) ~ N(0, T-t).
###############################################################################

"""
    mc_call_condexp_additive(S_t, t, T, K, μ, σ; M=50_000, rng=Random.GLOBAL_RNG)

Monte Carlo approximation of V_t = E[(S_T-K)^+ | S_t] in the additive model.
"""
function mc_call_condexp_additive(S_t::Float64, t::Float64, T::Float64,
                                  K::Float64, μ::Float64, σ::Float64;
                                  M::Int=50_000, rng=Random.GLOBAL_RNG)
    τ = T - t
    if τ <= 0
        return max(S_t - K, 0.0)
    end
    # Conditional: S_T = S_t + μ τ + σ sqrt(τ) Z
    s = 0.0
    for _ in 1:M
        Z = randn(rng)
        ST = S_t + μ * τ + σ * sqrt(τ) * Z
        s += max(ST - K, 0.0)
    end
    return s / M
end

# Demonstrate closed-form vs Monte Carlo at one time t = T/2 along our path
K = 0.0
t_mid = 0.5T
idx_mid = Int(round(t_mid / (T/N))) + 1  # grid index (1-based)
S_mid = S[idx_mid]

V_mid_closed = call_condexp_additive(S_mid, t_mid, T, K, μ, σ)
V_mid_mc     = mc_call_condexp_additive(S_mid, t_mid, T, K, μ, σ; M=100_000)

println("\nCall conditional expectation at t=T/2 along one path:")
println("  S_t = $(S_mid)")
println("  closed form V_t = $(V_mid_closed)")
println("  Monte Carlo V_t  = $(V_mid_mc)")