//+------------------------------------------------------------------+
//|                                                 BlackScholes.mqh |
//|                                      MMQ — Muhammad Minhas Qamar |
//|                                   www.mql5.com/en/articles/23385 |
//+------------------------------------------------------------------+
#property copyright "MMQ — Muhammad Minhas Qamar"
#property link      "https://www.mql5.com/en/articles/23385"
#property version   "1.00"

#ifndef IVSURFACE_BLACKSCHOLES_MQH
#define IVSURFACE_BLACKSCHOLES_MQH

//+------------------------------------------------------------------+
//| Option right, matching MQL5's ENUM_SYMBOL_OPTION_RIGHT so both   |
//| data providers can pass the value straight through.              |
//+------------------------------------------------------------------+
enum ENUM_OPT_RIGHT
  {
   OPT_CALL = 0,   // right to buy
   OPT_PUT  = 1    // right to sell
  };

//+------------------------------------------------------------------+
//| Standard normal CDF, via the Abramowitz & Stegun 7.1.26          |
//| approximation of erf. Accurate to about 1e-7, far tighter than   |
//| option-quote noise.                                              |
//+------------------------------------------------------------------+
double NormCDF(const double x)
  {
   double t = 1.0 / (1.0 + 0.2316419 * MathAbs(x));
   double d = 0.3989422804014327 * MathExp(-x * x / 2.0);   // 1/sqrt(2pi) * e^(-x^2/2)
   double p = d * t * (0.319381530 + t * (-0.356563782 + t * (1.781477937 + t * (-1.821255978 + t * 1.330274429))));
   return((x > 0.0) ? 1.0 - p : p);
  }

//+------------------------------------------------------------------+
//| Standard normal PDF.                                             |
//+------------------------------------------------------------------+
double NormPDF(const double x)
  {
   return(0.3989422804014327 * MathExp(-x * x / 2.0));
  }

//+------------------------------------------------------------------+
//| Black-Scholes price of a European option.                        |
//|   S     spot price of the underlying                             |
//|   K     strike price                                             |
//|   r     risk-free rate (annual, continuous)                      |
//|   q     dividend / carry yield (annual, continuous)              |
//|   sigma volatility (annual)                                      |
//|   T     time to expiry in years                                  |
//+------------------------------------------------------------------+
double BSPrice(const ENUM_OPT_RIGHT right, const double S, const double K,
               const double r, const double q, const double sigma, const double T)
  {
   if(T <= 0.0 || sigma <= 0.0)
     {
      //--- at/after expiry, value is pure intrinsic
      double intrinsic = (right == OPT_CALL) ? MathMax(S - K, 0.0) : MathMax(K - S, 0.0);
      return(intrinsic);
     }
   double sqrtT  = MathSqrt(T);
   double d1     = (MathLog(S / K) + (r - q + 0.5 * sigma * sigma) * T) / (sigma * sqrtT);
   double d2     = d1 - sigma * sqrtT;
   double disc_r = MathExp(-r * T);
   double disc_q = MathExp(-q * T);
   if(right == OPT_CALL)
      return(S * disc_q * NormCDF(d1) - K * disc_r * NormCDF(d2));
   else
      return(K * disc_r * NormCDF(-d2) - S * disc_q * NormCDF(-d1));
  }

//+------------------------------------------------------------------+
//| Vega (price sensitivity to a 1.00 change in sigma). Used both as |
//| the Newton step denominator and to reject near-worthless options |
//| whose IV is numerically meaningless.                             |
//+------------------------------------------------------------------+
double BSVega(const double S, const double K, const double r, const double q,
              const double sigma, const double T)
  {
   if(T <= 0.0 || sigma <= 0.0)
      return(0.0);
   double sqrtT = MathSqrt(T);
   double d1    = (MathLog(S / K) + (r - q + 0.5 * sigma * sigma) * T) / (sigma * sqrtT);
   return(S * MathExp(-q * T) * NormPDF(d1) * sqrtT);
  }

//+------------------------------------------------------------------+
//| Implied volatility by inverting the Black-Scholes price. Uses    |
//| Newton-Raphson (fast) and falls back to bisection whenever a     |
//| Newton step leaves the bracket, which keeps the solver robust for|
//| deep in/out-of-the-money quotes where vega is tiny.              |
//|                                                                  |
//| Returns the implied volatility, or a negative value if no valid  |
//| IV exists for the given inputs (e.g. price below intrinsic).     |
//+------------------------------------------------------------------+
double ImpliedVol(const ENUM_OPT_RIGHT right, const double price, const double S,
                  const double K, const double r, const double q, const double T)
  {
   if(price <= 0.0 || S <= 0.0 || K <= 0.0 || T <= 0.0)
      return(-1.0);

//--- price must sit above intrinsic value or no positive sigma solves it
   double intrinsic = (right == OPT_CALL) ? MathMax(S * MathExp(-q * T) - K * MathExp(-r * T), 0.0)
                      : MathMax(K * MathExp(-r * T) - S * MathExp(-q * T), 0.0);
   if(price < intrinsic - 1e-8)
      return(-1.0);

   double loVol = 1e-4, hiVol = 5.0;   // 0.01% .. 500% brackets the search
   double sigma = 0.2;                 // a sensible starting guess (20%)

   for(int i = 0; i < 100; i++)
     {
      double model = BSPrice(right, S, K, r, q, sigma, T);
      double diff  = model - price;
      if(MathAbs(diff) < 1e-6)
         return(sigma);

      //--- shrink the bracket using the sign of the error (price rises monotonically with sigma)
      if(diff > 0.0)
         hiVol = sigma;
      else
         loVol = sigma;

      double vega = BSVega(S, K, r, q, sigma, T);
      double next;
      if(vega > 1e-8)
         next = sigma - diff / vega;       // Newton step
      else
         next = 0.5 * (loVol + hiVol);     // vega too small: bisect

      //--- if Newton jumped outside the bracket, fall back to bisection
      if(next <= loVol || next >= hiVol)
         next = 0.5 * (loVol + hiVol);

      if(MathAbs(next - sigma) < 1e-8)
         return(next);
      sigma = next;
     }
   return(sigma);
  }

#endif // IVSURFACE_BLACKSCHOLES_MQH
//+------------------------------------------------------------------+