Monte Carlo pricing of Exotic Options in F#
open MathNet.Numerics.Distributions | |
open MathNet.Numerics.Statistics | |
let callPayoff strike price = max (price - strike) 0.0 | |
let europeanPayoff payoff assetPath = assetPath |> Seq.last |> payoff | |
let europeanCallPayoff strike assetPath = assetPath |> europeanPayoff (callPayoff strike) | |
let asianArithmeticMeanCallPayoff strike (assetPath:seq<float>) = assetPath.Mean() |> callPayoff strike | |
let upAndOutPayoff payoff upperBarrier assetPath = | |
if assetPath |> Seq.exists (fun x -> x > upperBarrier) then 0.0 else payoff assetPath | |
let upAndOutEuropeanCallPayoff strike = upAndOutPayoff (europeanCallPayoff strike) | |
let doubleBarrierPayoff payoff upperBarrier lowerBarrier assetPath = | |
if assetPath |> Seq.exists (fun x -> x > upperBarrier || x < lowerBarrier) then 0.0 else payoff assetPath | |
let doubleBarrierEuropeanCallPayoff strike = doubleBarrierPayoff (europeanCallPayoff strike) | |
let getAssetPath S0 r deltaT sigma (normal:Normal) numSamples = | |
Seq.unfold (fun S -> let sNew = (S * exp(((r - (0.5 * sigma * sigma)) * deltaT) + (sigma * sqrt(deltaT) * normal.Sample()))) | |
Some(sNew, sNew)) S0 | |
|> Seq.take numSamples | |
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 = new Normal() | |
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 = new Normal(0.0, 1.0) | |
let deltaT = T / float numSamples | |
let payoffs = seq { for n in 1 .. numTrajectories do | |
let assetPath = getAssetPath S0 r deltaT sigma normal numSamples |> Seq.toList | |
yield assetPath |> 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 = new Normal(0.0, 1.0) | |
let deltaT = T / float numSamples | |
let payoffs = seq { for n in 1 .. numTrajectories do | |
let assetPath = getAssetPath S0 r deltaT sigma normal numSamples |> Seq.toList | |
let price = assetPath |> doubleBarrierEuropeanCallPayoff strike upperBarrier lowerBarrier | |
yield price | |
} | |
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 : unit -> seq<float>) controlExact (payoff : seq<float> -> float) (controlPayoff) = | |
let samplePayoffs = seq { for n in 1 .. numTrajectories do | |
let assetPath = samplesFactory() |> Seq.toList | |
let payoff = (payoff assetPath, controlPayoff assetPath) | |
yield payoff | |
} | |
|> Seq.toList | |
let payoffs = samplePayoffs |> Seq.map (fun i -> fst i) |> Seq.toArray | |
let payoffsControlVariate = samplePayoffs |> Seq.map (fun i -> snd i) |> Seq.toArray | |
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 = new Normal(0.0, 1.0) | |
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 = new Normal(0.0, 1.0) | |
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 = | |
asian | |
barrier | |
0 |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
This comment has been minimized.
taumuon commentedApr 25, 2014
http://taumuon-jabuka.blogspot.co.uk/2014/04/monte-carlo-pricing-of-exotic-options.html