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

#ifndef GK_GARMANKOHLHAGEN_MQH
#define GK_GARMANKOHLHAGEN_MQH

//+------------------------------------------------------------------+
//| Option right. Kept numerically identical to MQL5's own           |
//| ENUM_SYMBOL_OPTION_RIGHT so the native provider can pass the     |
//| server's value straight through without a translation table.     |
//+------------------------------------------------------------------+
enum ENUM_OPT_RIGHT
  {
   OPT_CALL = 0,   // right to buy
   OPT_PUT  = 1    // right to sell
  };

//+------------------------------------------------------------------+
//| FX delta convention. FX vols are quoted against a delta, not a   |
//| strike, and the delta itself is convention-dependent: it may be  |
//| measured on the spot or the forward, and it may or may not be    |
//| premium-adjusted. The four combinations below place the same     |
//| 25-delta at four different strikes, so picking the wrong one     |
//| silently misprices the whole smile. Which one a pair uses is a   |
//| market convention, tabulated per pair, not a free choice.        |
//+------------------------------------------------------------------+
enum ENUM_GK_DELTA
  {
   GK_DELTA_SPOT,       // spot delta, unadjusted (pips)
   GK_DELTA_FORWARD,    // forward delta, unadjusted (pips)
   GK_DELTA_SPOT_PA,    // spot delta, premium-adjusted
   GK_DELTA_FORWARD_PA  // forward delta, premium-adjusted
  };

//+------------------------------------------------------------------+
//| ATM convention. "At-the-money" is not one strike. The delta-     |
//| neutral straddle is the strike where a call and a put share the  |
//| same absolute delta, and it depends on the delta convention in   |
//| use; the at-the-money-forward is simply the forward. A quoted    |
//| ATM vol means nothing until the strike it attaches to is fixed.  |
//+------------------------------------------------------------------+
enum ENUM_GK_ATM
  {
   GK_ATM_DNS,      // delta-neutral straddle
   GK_ATM_FORWARD   // at-the-money forward (K = F)
  };

//+------------------------------------------------------------------+
//| Standard normal CDF, via the Abramowitz & Stegun 7.1.26          |
//| approximation of erf. Accurate to about 1e-7, far tighter than   |
//| any FX vol-quote noise, and branch-free enough to stay cheap in  |
//| the strike-solver's inner loop.                                  |
//+------------------------------------------------------------------+
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));
  }

//+------------------------------------------------------------------+
//| The forward FX rate. A currency position earns the foreign rate  |
//| and funds at the domestic rate, so the drift is (rd - rf) and    |
//| the forward is the spot carried at that differential. Every      |
//| Garman-Kohlhagen quantity below is cleanest expressed through    |
//| this forward, so we compute it once here.                        |
//+------------------------------------------------------------------+
double GKForward(const double S, const double rd, const double rf, const double T)
  {
   return(S * MathExp((rd - rf) * T));
  }

//+------------------------------------------------------------------+
//| Garman-Kohlhagen price of a European FX option. This is Black-   |
//| Scholes with the foreign interest rate playing the role of a     |
//| continuous dividend yield on the base currency: the base leg is  |
//| discounted at rf, the strike leg at the domestic rate rd. At or  |
//| past expiry, or for a non-positive vol, value is pure intrinsic. |
//+------------------------------------------------------------------+
double GKPrice(const ENUM_OPT_RIGHT right, const double S, const double K,
               const double rd, const double rf, 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) + (rd - rf + 0.5 * sigma * sigma) * T) / (sigma * sqrtT);
   double d2    = d1 - sigma * sqrtT;
   double dfDom = MathExp(-rd * T);
   double dfFor = MathExp(-rf * T);
   if(right == OPT_CALL)
      return(S * dfFor * NormCDF(d1) - K * dfDom * NormCDF(d2));
   else
      return(K * dfDom * NormCDF(-d2) - S * dfFor * NormCDF(-d1));
  }

//+------------------------------------------------------------------+
//| Vega: sensitivity to a 1.00 (100 vol-point) change in sigma.     |
//| Both a call and a put on the same strike share it, and it is the |
//| Newton denominator in the inversion below. It also flags a       |
//| worthless quote: where vega vanishes, the implied vol is         |
//| numerically meaningless and the solver must not trust a step.    |
//+------------------------------------------------------------------+
double GKVega(const double S, const double K, const double rd, const double rf,
              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) + (rd - rf + 0.5 * sigma * sigma) * T) / (sigma * sqrtT);
   return(S * MathExp(-rf * T) * NormPDF(d1) * sqrtT);
  }

//+------------------------------------------------------------------+
//| Implied volatility by inverting the Garman-Kohlhagen price. This |
//| is a hybrid Newton-Raphson / bisection: Newton is fast where     |
//| vega is healthy, but deep in/out-of-the-money quotes have tiny   |
//| vega and a raw Newton step there can fly off to nonsense. We     |
//| keep a bracket [loVol, hiVol] that must contain the root         |
//| (price rises monotonically with sigma) and fall back to          |
//| bisection whenever Newton would leave it.                        |
//|                                                                  |
//| Returns the implied vol, or a negative sentinel when no positive |
//| sigma can reproduce the quote (e.g. a price below intrinsic).    |
//+------------------------------------------------------------------+
double GKImpliedVol(const ENUM_OPT_RIGHT right, const double price, const double S,
                    const double K, const double rd, const double rf, 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(-rf * T) - K * MathExp(-rd * T), 0.0)
                      : MathMax(K * MathExp(-rd * T) - S * MathExp(-rf * 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.1;                 // 10% is a sensible FX starting guess

   for(int i = 0; i < 100; i++)
     {
      double model = GKPrice(right, S, K, rd, rf, sigma, T);
      double diff  = model - price;
      if(MathAbs(diff) < 1e-8)
         return(sigma);

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

      double vega = GKVega(S, K, rd, rf, 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-10)
         return(next);
      sigma = next;
     }
   return(sigma);
  }

//+------------------------------------------------------------------+
//| Delta under a chosen FX convention. The unadjusted spot delta,   |
//| e^(-rf T) N(d1), is the raw hedge ratio in base-currency units.  |
//| Dropping the e^(-rf T) discount gives the forward delta. The     |
//| premium-adjusted variants subtract the option's own base-        |
//| currency premium from the hedge, which turns the N(d1) term into |
//| a (K/F) N(d2) term. Puts carry the mirror sign throughout.       |
//+------------------------------------------------------------------+
double GKDelta(const ENUM_OPT_RIGHT right, const ENUM_GK_DELTA conv, const double S,
               const double K, const double rd, const double rf, const double sigma, const double T)
  {
   if(T <= 0.0 || sigma <= 0.0 || S <= 0.0 || K <= 0.0)
      return(0.0);
   double sqrtT = MathSqrt(T);
   double F     = GKForward(S, rd, rf, T);
   double d1    = (MathLog(S / K) + (rd - rf + 0.5 * sigma * sigma) * T) / (sigma * sqrtT);
   double d2    = d1 - sigma * sqrtT;
   double phi   = (right == OPT_CALL) ? 1.0 : -1.0;
   double dfFor = MathExp(-rf * T);

   switch(conv)
     {
      case GK_DELTA_SPOT:
         return(phi * dfFor * NormCDF(phi * d1));
      case GK_DELTA_FORWARD:
         return(phi * NormCDF(phi * d1));
      case GK_DELTA_SPOT_PA:
         return(phi * dfFor * (K / F) * NormCDF(phi * d2));
      case GK_DELTA_FORWARD_PA:
         return(phi * (K / F) * NormCDF(phi * d2));
     }
   return(0.0);
  }

//+------------------------------------------------------------------+
//| Gamma: the rate of change of the spot delta per 1.0 move in the  |
//| spot. Discounted at the foreign rate, as the base-currency hedge |
//| is. A call and a put on the same strike/expiry share it.         |
//+------------------------------------------------------------------+
double GKGamma(const double S, const double K, const double rd, const double rf,
               const double sigma, const double T)
  {
   if(T <= 0.0 || sigma <= 0.0 || S <= 0.0)
      return(0.0);
   double sqrtT = MathSqrt(T);
   double d1    = (MathLog(S / K) + (rd - rf + 0.5 * sigma * sigma) * T) / (sigma * sqrtT);
   return(MathExp(-rf * T) * NormPDF(d1) / (S * sigma * sqrtT));
  }

//+------------------------------------------------------------------+
//| Theta: value lost per year as time passes (the sign is dV/dt, so |
//| a decaying long option returns a negative number). It carries    |
//| both carry terms an equity option lacks: the base leg bleeds at  |
//| rf and the strike leg accrues at rd, on top of the shared time-  |
//| decay of the volatility term.                                    |
//+------------------------------------------------------------------+
double GKTheta(const ENUM_OPT_RIGHT right, const double S, const double K, const double rd,
               const double rf, const double sigma, const double T)
  {
   if(T <= 0.0 || sigma <= 0.0 || S <= 0.0 || K <= 0.0)
      return(0.0);
   double sqrtT = MathSqrt(T);
   double d1    = (MathLog(S / K) + (rd - rf + 0.5 * sigma * sigma) * T) / (sigma * sqrtT);
   double d2    = d1 - sigma * sqrtT;
   double dfDom = MathExp(-rd * T);
   double dfFor = MathExp(-rf * T);
   double decay = -S * dfFor * NormPDF(d1) * sigma / (2.0 * sqrtT);
   if(right == OPT_CALL)
      return(decay + rf * S * dfFor * NormCDF(d1) - rd * K * dfDom * NormCDF(d2));
   else
      return(decay - rf * S * dfFor * NormCDF(-d1) + rd * K * dfDom * NormCDF(-d2));
  }

//+------------------------------------------------------------------+
//| Rho to the domestic rate: sensitivity of value to a 1.00 change  |
//| in rd, acting through the discounted strike leg. FX has two of   |
//| these, one per rate, which is the structural break from the      |
//| single-rho equity Greeks.                                        |
//+------------------------------------------------------------------+
double GKRhoDom(const ENUM_OPT_RIGHT right, const double S, const double K, const double rd,
                const double rf, const double sigma, const double T)
  {
   if(T <= 0.0 || sigma <= 0.0 || S <= 0.0 || K <= 0.0)
      return(0.0);
   double sqrtT = MathSqrt(T);
   double d1    = (MathLog(S / K) + (rd - rf + 0.5 * sigma * sigma) * T) / (sigma * sqrtT);
   double d2    = d1 - sigma * sqrtT;
   double dfDom = MathExp(-rd * T);
   if(right == OPT_CALL)
      return(K * T * dfDom * NormCDF(d2));
   else
      return(-K * T * dfDom * NormCDF(-d2));
  }

//+------------------------------------------------------------------+
//| Rho to the foreign rate: sensitivity of value to a 1.00 change   |
//| in rf, acting through the discounted base leg. It is the second  |
//| interest-rate Greek, and it moves opposite to the domestic one.  |
//+------------------------------------------------------------------+
double GKRhoFor(const ENUM_OPT_RIGHT right, const double S, const double K, const double rd,
                const double rf, const double sigma, const double T)
  {
   if(T <= 0.0 || sigma <= 0.0 || S <= 0.0 || K <= 0.0)
      return(0.0);
   double sqrtT = MathSqrt(T);
   double d1    = (MathLog(S / K) + (rd - rf + 0.5 * sigma * sigma) * T) / (sigma * sqrtT);
   double dfFor = MathExp(-rf * T);
   if(right == OPT_CALL)
      return(-S * T * dfFor * NormCDF(d1));
   else
      return(S * T * dfFor * NormCDF(-d1));
  }

//+------------------------------------------------------------------+
//| The ATM strike for a given ATM and delta convention. The delta-  |
//| neutral straddle sits half a variance above the forward under an |
//| unadjusted delta, and half a variance below it under a premium-  |
//| adjusted delta, because the premium term shifts the delta that   |
//| has to net to zero. The forward convention is just F. Getting    |
//| this strike wrong offsets the entire smile by its ATM anchor.    |
//+------------------------------------------------------------------+
double GKATMStrike(const ENUM_GK_ATM atmConv, const ENUM_GK_DELTA deltaConv,
                   const double F, const double sigma, const double T)
  {
   if(atmConv == GK_ATM_FORWARD)
      return(F);
//--- delta-neutral straddle: sign of the half-variance shift depends
//--- on whether the delta is premium-adjusted
   bool pa = (deltaConv == GK_DELTA_SPOT_PA || deltaConv == GK_DELTA_FORWARD_PA);
   double half = 0.5 * sigma * sigma * T;
   return(pa ? F * MathExp(-half) : F * MathExp(half));
  }

//+------------------------------------------------------------------+
//| Invert a target delta to its strike, holding sigma fixed. Used   |
//| to turn each quoted delta pillar (10d, 25d) into the strike the  |
//| smile actually lives on. The unadjusted delta is monotonic in    |
//| the strike, but the premium-adjusted call delta is not: it rises |
//| then falls, so a 25-delta call has two strikes and the market    |
//| convention takes the larger (out-of-the-money) one. Scanning the |
//| strike axis from high to low and bracketing the first sign       |
//| change returns exactly that root, and the unique root otherwise. |
//|                                                                  |
//| targetAbs is the delta magnitude (e.g. 0.25); the sign is taken  |
//| from the right. Returns a negative sentinel if no strike matches.|
//+------------------------------------------------------------------+
double GKStrikeFromDelta(const ENUM_OPT_RIGHT right, const ENUM_GK_DELTA conv, const double targetAbs,
                         const double S, const double rd, const double rf, const double sigma, const double T)
  {
   if(T <= 0.0 || sigma <= 0.0 || S <= 0.0 || targetAbs <= 0.0)
      return(-1.0);
   double target = (right == OPT_CALL) ? targetAbs : -targetAbs;
   double F      = GKForward(S, rd, rf, T);

//--- scan log-strike space around the forward, high to low, for a
//--- sign change of (delta - target); K in [F*e^-3, F*e^+3] is ample
   int    steps = 400;
   double kHi   = F * MathExp(3.0);
   double kLo   = F * MathExp(-3.0);
   double logHi = MathLog(kHi), logLo = MathLog(kLo);
   double prevK = kHi;
   double prevF = GKDelta(right, conv, S, prevK, rd, rf, sigma, T) - target;

   for(int i = 1; i <= steps; i++)
     {
      double frac = (double)i / steps;
      double curK = MathExp(logHi + (logLo - logHi) * frac);
      double curF = GKDelta(right, conv, S, curK, rd, rf, sigma, T) - target;
      if(prevF == 0.0)
         return(prevK);
      if(prevF * curF < 0.0)
        {
         //--- bracket found: bisect [curK, prevK] to a tight strike
         double a = curK, b = prevK, fa = curF;
         for(int j = 0; j < 100; j++)
           {
            double m  = 0.5 * (a + b);
            double fm = GKDelta(right, conv, S, m, rd, rf, sigma, T) - target;
            if(MathAbs(fm) < 1e-12 || (b - a) < 1e-10)
               return(m);
            if(fa * fm < 0.0)
               b = m;
            else
              {
               a  = m;
               fa = fm;
              }
           }
         return(0.5 * (a + b));
        }
      prevK = curK;
      prevF = curF;
     }
   return(-1.0);
  }

#endif // GK_GARMANKOHLHAGEN_MQH
//+------------------------------------------------------------------+