Skip to content

Instantly share code, notes, and snippets.

@taumuon
Created April 25, 2014 20:52
Show Gist options
  • Save taumuon/11302896 to your computer and use it in GitHub Desktop.
Save taumuon/11302896 to your computer and use it in GitHub Desktop.
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