Also I think it may be interesting that naive realization of Black-Sholes can be 32X faster. An example for the case with interest rate = 0, dividend = 0:

Also I think it may be interesting that naive realization of Black-Scholes can be 32X faster. An example for the case with interest rate = 0, dividend = 0:

Also I think it may be interesting that naive realization of Black-Sholes can be 16X faster. An example for the case with rate = 0:

call[strike_, spot_, time_, vola_] := Module[
  {
   d1 = (Log[spot/strike] + vola^2/2*time/365)/(vola*Sqrt[time/365]),
   d2 = (Log[spot/strike] - vola^2/2*time/365)/(vola*Sqrt[time/365]),
   cdf = .5 Erfc[-#/Sqrt[2]] &,
   dtcdf = Exp[-#^2/2]/Sqrt[2*Pi] &
   },
  {
   spot*cdf[d1] - strike*cdf[d2],
   (*delta*)cdf[d1],
   (*gamma*)dtcdf[d1]/(spot*vola*Sqrt[time/365]),
   (*vega*)spot*dtcdf[d1]*Sqrt[time/365]/100,
   (*theta*)-spot*dtcdf[d1]*vola/(2*Sqrt[time/365])/365
   }
  ]

put[strike_, spot_, time_, vola_] := Module[
  {
   d1 = (Log[spot/strike] + vola^2/2*time/365)/(vola*Sqrt[time/365]),
   d2 = (Log[spot/strike] - vola^2/2*time/365)/(vola*Sqrt[time/365]),
   cdf = .5 Erfc[-#/Sqrt[2]] &,
   dtcdf = Exp[-#^2/2]/Sqrt[2*Pi] &
   },
  {
   strike*cdf[-d2] - spot*cdf[-d1],
   (*delta*)cdf[d1] - 1,
   (*gamma*)dtcdf[d1]/(spot*vola*Sqrt[time/365]),
   (*vega*)spot*dtcdf[d1]*Sqrt[time/365]/100,
   (*theta*)-spot*dtcdf[d1]*vola/(2*Sqrt[time/365])/365
   }
  ]

Also I think it may be interesting that naive realization of Black-Sholes can be 32X faster. An example for the case with interest rate = 0, dividend = 0:

call[strike_, spot_, time_, vola_] := Module[
  {
   d1 = (Log[spot/strike] + vola^2/2*time/365)/(vola*Sqrt[time/365]),
   d2 = (Log[spot/strike] - vola^2/2*time/365)/(vola*Sqrt[time/365]),
   cdf = .5 Erfc[-#/Sqrt[2]] &,
   dtcdf = Exp[-#^2/2]/Sqrt[2*Pi] &
   },
  {
   spot*cdf[d1] - strike*cdf[d2],
   (*delta*)cdf[d1],
   (*gamma*)dtcdf[d1]/(spot*vola*Sqrt[time/365]),
   (*vega*)spot*dtcdf[d1]*Sqrt[time/365]/100,
   (*theta*)-spot*dtcdf[d1]*vola/(2*Sqrt[time/365])/365
   }
  ]

put[strike_, spot_, time_, vola_] := Module[
  {
   d1 = (Log[spot/strike] + vola^2/2*time/365)/(vola*Sqrt[time/365]),
   d2 = (Log[spot/strike] - vola^2/2*time/365)/(vola*Sqrt[time/365]),
   cdf = .5 Erfc[-#/Sqrt[2]] &,
   dtcdf = Exp[-#^2/2]/Sqrt[2*Pi] &
   },
  {
   strike*cdf[-d2] - spot*cdf[-d1],
   (*delta*)cdf[d1] - 1,
   (*gamma*)dtcdf[d1]/(spot*vola*Sqrt[time/365]),
   (*vega*)spot*dtcdf[d1]*Sqrt[time/365]/100,
   (*theta*)-spot*dtcdf[d1]*vola/(2*Sqrt[time/365])/365
   }
  ]

Checking the correctness:

In[9]:= put[87500, 89040, 13, .33966]

Out[9]= {1569.98, -0.380478, 0.0000667356, 64.0063, -83.6168}

In[10]:= FinancialDerivative[
  {"European", "Put"},
   {"StrikePrice" -> 87500, 
   "Expiration" -> 13/365},  {"InterestRate" -> 0., 
   "Volatility" -> 0.33966, "CurrentPrice" -> 89040, "Dividend" -> 0.},
  {"Value", "Delta", "Gamma", "Vega", "Theta"}]*{1, 1, 1, 1/100, 1/365}

Out[10]= {1569.98, -0.380478, 0.0000667356, 64.0062, -83.6164}

Benchmarking:

In[11]:= Do[put[87500, 89040, 13, .33966], {1000}] // AbsoluteTiming

Out[11]= {0.159938, Null}

In[12]:= Do[
  FinancialDerivative[
    {"European", "Put"},
     {"StrikePrice" -> 87500, 
     "Expiration" -> 13/365},  {"InterestRate" -> 0., 
     "Volatility" -> .33966, "CurrentPrice" -> 89040, 
     "Dividend" -> 0.},
    {"Value", "Delta", "Gamma", "Vega", "Theta"}]*{1, 1, 1, 1/100, 
    1/365}
  , {1000}] // AbsoluteTiming

Out[12]= {4.79136, Null}