package main import ( "math" "math/rand" "time" ) // Option struct that holds all the details of an option type Option struct { OptType int64 Strike float64 ExpiryDate time.Time // Market data ValueDate time.Time Spot float64 Rfr float64 Vol float64 Sims int64 // Results FV float64 Delta float64 Vega float64 Rho float64 Gamma float64 Theta float64 Levels []float64 } // Single simulation of Geometric Brownian Motion func gbmSimulation(spot float64, rfr float64, vol float64, tte float64, randNum float64) float64 { var drift, stoch float64 drift = (rfr - math.Pow(vol, 2)/2) * tte stoch = vol * math.Pow(tte, 0.5) * randNum return spot * math.Exp(drift+stoch) } // PriceMonteCarlo is the exported method to run the simulations and value // a vanilla option using Monte Carlo methods func (opt *Option) PriceMonteCarlo() { var i int64 var level float64 opt.Levels = make([]float64, opt.Sims) tte := opt.ExpiryDate.Sub(opt.ValueDate).Hours() / (24 * 365) for i = 0; i < opt.Sims; i++ { randNum := rand.NormFloat64() // Base opt.Levels[i] = gbmSimulation(opt.Spot, opt.Rfr, opt.Vol, tte, randNum) opt.FV += math.Max((opt.Levels[i]-opt.Strike)*float64(opt.OptType), 0) // Delta level = gbmSimulation(opt.Spot+0.0001, opt.Rfr, opt.Vol, tte, randNum) opt.Delta += math.Max((level-opt.Strike)*float64(opt.OptType), 0) // Gamma -- TODO: Doesn't look right level = gbmSimulation(opt.Spot+0.0001, opt.Rfr, opt.Vol, tte, randNum) level += gbmSimulation(opt.Spot-0.0001, opt.Rfr, opt.Vol, tte, randNum) level -= 2 * gbmSimulation(opt.Spot, opt.Rfr, opt.Vol, tte, randNum) opt.Gamma += math.Max((level-opt.Strike)*float64(opt.OptType), 0) // Vega level = gbmSimulation(opt.Spot, opt.Rfr, opt.Vol+0.0001, tte, randNum) opt.Vega += math.Max((level-opt.Strike)*float64(opt.OptType), 0) // Theta level = gbmSimulation(opt.Spot, opt.Rfr, opt.Vol, tte-1./365, randNum) opt.Theta += math.Max((level-opt.Strike)*float64(opt.OptType), 0) // Rho -- TODO: Doesn't look right level = gbmSimulation(opt.Spot, opt.Rfr+0.0001, opt.Vol, tte, randNum) opt.Rho += math.Max((level-opt.Strike)*float64(opt.OptType), 0) } df := math.Exp(-opt.Rfr * tte) opt.FV = opt.FV / float64(opt.Sims) * df opt.Delta = (opt.Delta/float64(opt.Sims)*df - opt.FV) / 0.0001 opt.Gamma = (opt.Gamma/float64(opt.Sims)*df - opt.FV) / 10000 opt.Vega = (opt.Vega/float64(opt.Sims)*df - opt.FV) / 0.01 opt.Theta = (opt.Theta/float64(opt.Sims)*math.Exp(-opt.Rfr*(tte-1./365)) - opt.FV) / -(1. / 365) opt.Rho = (opt.Rho/float64(opt.Sims)*math.Exp(-(opt.Rfr+0.01)*tte) - opt.FV) }