Skip to content

Instantly share code, notes, and snippets.

@jack-pappas
Forked from taumuon/MonteCarloExoticOptions
Last active August 29, 2015 14:04
Show Gist options
  • Save jack-pappas/d629a767823ca2f17967 to your computer and use it in GitHub Desktop.
Save jack-pappas/d629a767823ca2f17967 to your computer and use it in GitHub Desktop.
// Original code: https://gist.github.com/taumuon/11302896
// This implementation optimized by Jack Pappas (https://github.com/jack-pappas)
open MathNet.Numerics.Distributions
open MathNet.Numerics.Statistics
open MathNet.Numerics.Random
open LanguagePrimitives
let inline callPayoff strike price =
max (price - strike) GenericZero
let europeanPayoff payoff (assetPath : _[]) =
// Or, using ExtCore: assetPath |> Array.last |> payoff
payoff assetPath.[assetPath.Length - 1]
let europeanCallPayoff strike assetPath =
assetPath |> europeanPayoff (callPayoff strike)
let asianArithmeticMeanCallPayoff strike (assetPath : float[]) =
assetPath |> Array.average |> callPayoff strike
let upAndOutPayoff payoff (upperBarrier : float) assetPath =
if assetPath |> Array.exists (fun x -> x > upperBarrier) then GenericZero else payoff assetPath
let inline upAndOutEuropeanCallPayoff strike upperBarrier assetPath =
upAndOutPayoff (europeanCallPayoff strike) upperBarrier assetPath
let doubleBarrierPayoff payoff upperBarrier (lowerBarrier : float) assetPath =
if assetPath |> Array.exists (fun x -> x > upperBarrier || x < lowerBarrier) then GenericZero else payoff assetPath
let inline doubleBarrierEuropeanCallPayoff strike upperBarrier lowerBarrier assetPath =
doubleBarrierPayoff (europeanCallPayoff strike) upperBarrier lowerBarrier assetPath
let getAssetPath S0 r deltaT sigma (normal:Normal) numSamples =
let mutable S = S0
let normalSamples =
let normalSamples = Array.zeroCreate numSamples
normal.Samples normalSamples
normalSamples
(S0, normalSamples)
||> Array.scan (fun S normalSample ->
S * exp(((r - (0.5 * sigma * sigma)) * deltaT) + (sigma * sqrt(deltaT) * normalSample)))
let phi x =
let normal = new Normal()
normal.CumulativeDistribution(x)
// Exact option pricing:
let priceEuropeanCall S0 strike r T sigma =
let discountedK = strike * exp(-r * T)
let totalVolatility = sigma * sqrt(T)
let d_minus = log(S0 / discountedK) / totalVolatility
let d_plus = d_minus + (0.5 * totalVolatility)
let d_minus2 = d_minus - (0.5 * totalVolatility)
(S0 * phi(d_plus)) - (discountedK * phi(d_minus2))
let computeD S0 strike r sigma T =
let d1 = (log (S0 / strike) + ((r + (sigma * sigma / 2.0)) * T)) / (sigma * sqrt(T))
let d0 = d1 - (sigma * sqrt(T))
(d0, d1)
let priceBarrierUpAndOutCall S0 strike B r T sigma =
let normal = Normal(MersenneTwister ())
let mu = r - (0.5 * sigma * sigma)
let nu = 2.0 * (mu / (sigma * sigma))
let aux1 = (B / S0) ** (nu + 2.0)
let aux2 = (B / S0) ** nu
let d = computeD S0 strike r sigma T
let price = S0 * normal.CumulativeDistribution(snd d) - (strike * exp(-r * T) * normal.CumulativeDistribution(fst d))
let d2 = computeD S0 B r sigma T
let price2 = price - (S0 * normal.CumulativeDistribution(snd d2)) + (strike * exp(-r * T) * normal.CumulativeDistribution(fst d2))
let d3 = computeD (1.0/S0) (1.0/B) r sigma T
let price3 = price2 + (aux1 * S0 * normal.CumulativeDistribution(snd d3)) - (aux2 * strike * exp(-r * T) * normal.CumulativeDistribution(fst d3))
let d4 = computeD (1.0 / (strike * S0)) (1.0 / (B * B)) r sigma T
let price4 = price3 - (aux1 * S0 * normal.CumulativeDistribution(snd d4)) + (aux2 * strike * exp(-r * T) * normal.CumulativeDistribution(fst d4))
price4
// Monte Carlo pricing
let priceAsianArithmeticMeanMC S0 strike r T sigma numTrajectories numSamples =
let normal = Normal(MersenneTwister ())
let deltaT = T / float numSamples
let payoffs =
Array.init numTrajectories <| fun _ ->
getAssetPath S0 r deltaT sigma normal numSamples
|> asianArithmeticMeanCallPayoff strike
let discountFactor = exp(-r * T)
let priceMC = discountFactor * payoffs.Mean()
let stddevMC = discountFactor * payoffs.StandardDeviation() / sqrt(float numTrajectories)
(priceMC, stddevMC)
let priceDoubleBarrierCallMC S0 strike upperBarrier lowerBarrier r T sigma numTrajectories numSamples =
let normal = Normal(MersenneTwister ())
let deltaT = T / float numSamples
let payoffs =
Array.init numTrajectories <| fun _ ->
getAssetPath S0 r deltaT sigma normal numSamples
|> doubleBarrierEuropeanCallPayoff strike upperBarrier lowerBarrier
let discountFactor = exp(-r * T)
let priceMC = discountFactor * payoffs.Mean()
let stddevMC = discountFactor * payoffs.StandardDeviation() / sqrt(float numTrajectories)
(priceMC, stddevMC)
let controlVariate numTrajectories discountFactor samplesFactory controlExact payoff controlPayoff =
let samplePayoffs =
Array.init numTrajectories <| fun _ ->
let assetPath : _[] = samplesFactory()
payoff assetPath, controlPayoff assetPath
let payoffs, payoffsControlVariate =
samplePayoffs
|> Array.unzip
let variance = ArrayStatistics.PopulationCovariance(payoffs, payoffs)
let varianceControlVariate = ArrayStatistics.PopulationCovariance(payoffsControlVariate, payoffsControlVariate)
let covariance = ArrayStatistics.PopulationCovariance(payoffs, payoffsControlVariate)
let correlation = covariance / sqrt(varianceControlVariate * variance)
let priceMC = discountFactor * payoffs.Mean()
let stddevMC = discountFactor * payoffs.StandardDeviation() / sqrt(float numTrajectories)
let priceControlVariateMC = discountFactor * payoffsControlVariate.Mean()
let stddevControlVariateMC = discountFactor * payoffsControlVariate.StandardDeviation() / sqrt(float numTrajectories)
let price = priceMC - ((covariance / varianceControlVariate) * (priceControlVariateMC - controlExact))
let stddev = stddevMC * sqrt(1.0 - (correlation * correlation))
(price, stddev)
let priceDoubleBarrierCall_MC_CV_EU S0 strike upperBarrier lowerBarrier r T sigma numTrajectories numSamples =
let normal = Normal(MersenneTwister ())
let discountFactor = exp(-r * T)
let deltaT = T / float numSamples;
let priceEuropeanExact = priceEuropeanCall S0 strike r T sigma
let europeanCallfn = (fun assetPath -> europeanCallPayoff strike assetPath)
let barrierCallPayoff = (fun assetPath -> doubleBarrierEuropeanCallPayoff strike upperBarrier lowerBarrier assetPath)
controlVariate numTrajectories
discountFactor
(fun () -> getAssetPath S0 r deltaT sigma normal numSamples)
priceEuropeanExact
barrierCallPayoff
europeanCallfn
let priceDoubleBarrierCall_MC_CV_UpAndOut S0 strike upperBarrier lowerBarrier r T sigma numTrajectories numSamples =
let normal = Normal(MersenneTwister ())
let discountFactor = exp(-r * T)
let deltaT = T / float numSamples;
let priceUpAndOutExact = priceBarrierUpAndOutCall S0 strike upperBarrier r T sigma
let barrierCallPayoff = (fun assetPath -> doubleBarrierEuropeanCallPayoff strike upperBarrier lowerBarrier assetPath)
let upAndOut = (fun assetPath -> upAndOutEuropeanCallPayoff strike upperBarrier assetPath)
controlVariate numTrajectories
discountFactor
(fun () -> getAssetPath S0 r deltaT sigma normal numSamples)
priceUpAndOutExact
barrierCallPayoff
upAndOut
let asian () =
let S0 = 100.0
let strike = 90.0
let r = 0.05
let T = 1.0
let sigma = 0.2
let numTrajectories = 100000
let numSamples = 12
let result = priceAsianArithmeticMeanMC S0 strike r T sigma numTrajectories numSamples
printfn "Asian arithmetic mean:%f stddev:%f" (fst result) (snd result)
let barrier () =
let S0 = 100.0
let strike = 90.0
let upperBarrier = 160.0
let lowerBarrier = 75.0
let r = 0.05
let T = 0.5
let sigma = 0.4
let numTrajectories = 1000
let numSamples = 10000
let result = priceDoubleBarrierCallMC S0 strike upperBarrier lowerBarrier r T sigma numTrajectories numSamples
printfn "Monte Carlo double barrier call:%f stddev:%f" (fst result) (snd result)
let resultMCCVEU = priceDoubleBarrierCall_MC_CV_EU S0 strike upperBarrier lowerBarrier r T sigma numTrajectories numSamples
printfn "Control Variate with European Call price:%f stddev:%f" (fst resultMCCVEU) (snd resultMCCVEU)
let resultMCCVUpAndOut = priceDoubleBarrierCall_MC_CV_UpAndOut S0 strike upperBarrier lowerBarrier r T sigma numTrajectories numSamples
printfn "Control Variate with Up And Out price:%f stddev:%f" (fst resultMCCVUpAndOut) (snd resultMCCVUpAndOut)
[<EntryPoint>]
let main argv =
let watch = System.Diagnostics.Stopwatch.StartNew ()
asian ()
let asianTime = watch.Elapsed
watch.Restart ()
barrier ()
watch.Stop ()
let barrierTime = watch.Elapsed
printfn "Asian option compute time: %04fms" asianTime.TotalMilliseconds
printfn "Barrier option compute time: %04fms" barrierTime.TotalMilliseconds
0
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment