Risk neutral derivative pricing

fit.jl

using DelimitedFiles, GLMakie, Distributions

# Black-Scholes formula
normcdf(x) = cdf(Normal(), x)
function blcall(S0, K, r, T, σ)
    B = exp(r*T)
    F = S0 * B
    d1 = log(F/K) / (σ*sqrt(T)) + σ*sqrt(T)/2
    (F*normcdf(d1) - K*normcdf((d1 - σ*sqrt(T))))/B
end

# load data from 10 year csv downloaded at https://fred.stlouisfed.org/series/SP500
data_, header = readdlm("SP500.csv", ',', header=true)
Y = 5 # number of years
data = data_[end - Y*end÷10 + 1:end, :]

# data processing
T = Y
S = [a for a in data[:, 2] if a isa Number]
logS = log.(S)
N = length(S) # number of business days in Y years
ts = range(0, Y, length=N)
inc(x) = x[end] - x[begin]
dt = ts[2] - ts[1] # length of a business day as fraction of the year

# Visually check Gaussianity of the daily log returns
qqplot(Normal, diff(logS))

# fit geometric diffusion dS = μSdt + σSdW
σ2 = sum((log(S[i+1]) - log(S[i]))^2 for i in 1:length(S)-1)/ts[end]
σ = sqrt(σ2)
μ = inc(logS)/T
@show μ, σ

# Prediction
ΔT = 1.0 # in a year
D = LogNormal(log(S[end]) + μ*ΔT, σ*√(ΔT))
@show mean(D), std(D)

# Risk neutral process dS = rSdt + σSdW
r = log(1.043) # using rate 43% on US bonds

# Price of a call with strike K today in 5 business days
K = 6120 
blcall(S[end], K, r, 5/250, σ)

tma285fk.jl

using Distributions
n = 100^2
m = 50
xs = range(-5, 5, length=m)
ts = range(0, 0.5, length=n)
Δt = ts[2] - ts[1]
Δx = xs[2] - xs[1]

i0 = m÷2  
x0 = xs[i0] # starting value of Brownian motion (to generalize,...)

# ∂ₜ u + 0.5 ∂²ₓ u = 0
# u(T, x) = g(x) 

g(x) = pdf(Normal(), x)

# Financial derivative: Pay out of g(X) at time T, X = W 

# W(T) ~ N(0, T) E g(W(T)) = ∫ g(x) pdf(Normal(0, T), x) dx    


# to generalize, replace this by forward simulation of a stochastic differential equation starting in x0 
expect = sum(g(x)*pdf(Normal(0, sqrt(ts[end])), x)*Δx for x in xs)




u = zeros(n, m)
begin 
    u[end, :] = g.(xs)  # right boundary condition 

    for i in n:-1:2
        u[i-1, 1] = u[i-1, end] = 0.0 # to generalize, this needs care if g(x) > 0 at the boundary
        for j in 2:m-1
            # to generalize, replace by feynman kac formula for the discounted price
            u[i-1, j] =-(-u[i, j] - 0.5*Δt*(u[i, j+1] - 2u[i,j] + u[i, j-1])/(Δx^2))
        end
    end
end
heatmap(u)
u[1, i0] # u(0, x0)
expect

tma285mc.jl

using Distributions
using ForwardDiff
using ProgressMeter

# compute the value of the call using the Black Scholes formula
normcdf(x) = cdf(Normal(), x)
function blcall(S0, K, r, T, σ)
    B = exp(r*T)
    F = S0 * B
    d1 = log(F/K) / (σ*sqrt(T)) + σ*sqrt(T)/2
    (F*normcdf(d1) - K*normcdf((d1 - σ*sqrt(T))))/B
end

# Market parameter 
T = 1.0
N = 500
dt = 1/N
ts = range(0, 1, length=N+1)
iters = 2000

μ(t, x) = 0.2
r(t) = 0.05
σ(t, x) = 1.0

S0 = 2.

# European call
K = 3.0

# payout
Ycall(x) = max(x - K, 0.)
Yput(x) = max(K - x, 0.)

# value of the call according to black scholes (valid for constant r, σ)
V(t,x) = blcall(x, K, 0.05, T-t, 1.0)


# Riskless asset
B = zeros(N+1)
B0 = 1.0
B[1] = B0
for i in 1:N
    B[i+1] = B[i] + r(ts[i])*B[i]*dt
end
B[end]

# discount process
D(Δ) = exp(-0.05*Δ)
@show 1/D(T-0), B[end] 


## simulations using the physical measure

# simulate a stock and some calls
Ycalls = Float64[]
Yputs = Float64[]
S = zeros(N+1)
V(t,x) = blcall(x, K, 0.05, T-t, 1.0)
Vs = Float64[]

# replicating portfolio
a = ForwardDiff.derivative(x -> V(0.0, x), S0)
b = (V(0.0, S0) - a*S0)/B0
V0 = a*S0 + b*B0

@showprogress for iter in 1:iters
    global S
    S[1] = S0
    Viter = V0
    for i in 1:N
        dW = sqrt(dt)*rand(Normal())
        S[i+1] = S[i] + μ(ts[i], S[i])*S[i]*dt + σ(ts[i], S[i])*S[i]*dW  
        a = ForwardDiff.derivative(x -> V(ts[i]::Float64, x), S[i]::Float64)
        b = (V(ts[i], S[i]) - a*S[i])/B[i]
        Viter += a*(S[i+1]-S[i]) + b*(B[i+1]-B[i])
    end
    push!(Ycalls, Ycall(S[end]))
    push!(Yputs, Yput(S[end]))
    push!(Vs, Viter)
end

# average payout and standard deviation
mean(Ycalls), std(Ycalls)

# discounted call value at 0 not taking into account risk
mean(D(T-0)*Ycalls)

# replication: show the payout of the replicating portfolio equals the payout of the call
Vs[1:10]
Ycalls[1:10]


# discounted call value at 0 taking into account risk
utility(x) = x^0.9 # risk averse
inv_utility(x) = x^(1/0.9)
mean(utility.(D(T-0)*Ycalls))
# as much utility as the following amount cash a time 0
inv_utility(mean(utility.(D(T-0)*Ycalls)))


## simulations using the risk neutral measure

# simulate a stock under the risk neutral measure
Ỹcalls = Float64[]
Ỹputs = Float64[]
Ṽs = Float64[]

# again the replicating portfolio
a = ForwardDiff.derivative(x -> V(0.0, x), S0)
b = (V(0.0, S0) - a*S0)/B0
V0 = a*S0 + b*B0

@showprogress for iter in 1:iters
    global S
    S[1] = S0
    Viter = V0
    for i in 1:N
        dW̃ = sqrt(dt)*rand(Normal())
        S[i+1] = S[i] + r(ts[i])*S[i]*dt + σ(ts[i], S[i])*S[i]*dW̃

        a = ForwardDiff.derivative(x -> V(ts[i]::Float64, x), S[i]::Float64)
        b = (V(ts[i], S[i]) - a*S[i])/B[i]

        # self-financing by construction
        Viter += a*(S[i+1]-S[i]) + b*(B[i+1]-B[i])
    end
    push!(Ỹcalls, Ycall(S[end]))
    push!(Ỹputs, Yput(S[end]))
    push!(Ṽs, Viter) 
end

# discounted payout under the risk neutral measure
mean(D(T-0)*Ỹcalls)
blcall(S0, K, 0.05, T, 1.0)

# put call parity 
mean(D(T-0)*Ỹcalls) - mean(D(T-0)*Ỹputs), S0 - K*D(T-0) 

# replicating portfolio 
Ṽs[1:10]
Ỹcalls[1:10]
mean(D(T-0)*Vs)