//+------------------------------------------------------------------+
//|                    Robust_Volatility_Core.mqh                    |
//|                                                                  |
//|  Volatility-core module for the robust STARC / Z-Score, taken    |
//|  from StarcBands_Robust_MTF_Unified.mq4 (which compiles & loads).|
//|  It outputs, per bar, the volatility HALF-WIDTH in price units   |
//|  -> it REPLACES the ATR block of the Forex Station STARC.        |
//|                                                                  |
//|  ================ WHERE THIS PLUGS IN ========================   |
//|  Your STARC today:  atr   = YourChosenMA( TrueRange )            |
//|                     band  = midline +/- mult * atr               |
//|  With this module:  vol   = ComputeVolatilityAt(shift, seed)     |
//|                     band  = midline +/- mult * vol               |
//|  i.e. "ATR averaging type" + "ATR period" are REPLACED by        |
//|  'VolatilityMethod' + 'VolatilityPeriod'. The midline MA         |
//|  (your "Starc average type") is NOT changed here; it is only     |
//|  READ from g_center[].                                           |
//|                                                                  |
//|  The VolatilityMethod families have different dependencies:      |
//|   - True-Range based (median/Tukey/Huber of TR, log-TR):         |
//|       need only g_trueRange / g_logTrueRange                     |
//|   - RESIDUAL based (Qn/Sn/MAD/IQR/Tukey-scale of price-midline): |
//|       need the MIDLINE already in g_center                       |
//|   - RANGE based (Parkinson/Garman-Klass/Rogers-Satchell/Yang-    |
//|       Zhang): need OHLC in g_sourceOpen/High/Low/Close           |
//|   - RETURN based (close-to-close/EWMA/GARCH): use log returns    |
//|                                                                  |
//|  ==================== HOW TO USE (host) ======================   |
//|  1) In OnInit:  ResolvePairwiseMode();                           |
//|  2) When the bar count changes: RSE_ResizeVolatility(N);         |
//|  3) Size and FILL (index 0 = newest bar; series indexing):       |
//|       g_sourceOpen/High/Low/Close[], g_sourceRates               |
//|       g_price[]        = applied price (Close, HL2, ...)          |
//|       g_trueRange[]    = bar True Range                          |
//|       g_logTrueRange[] = ln( max(H,Cprev) / min(L,Cprev) )       |
//|       g_center[]       = YOUR midline (Starc average) at the bar |
//|  4) RESIDUAL-based methods (Qn/Sn/MAD/IQR/...): once the midline |
//|     is in g_center, call PrepareResidualData(oldest).            |
//|  5) vol = ComputeVolatilityAt(shift, seed);                      |
//|       - window/residual/range methods: seed can be false.       |
//|       - RECURSIVE methods (Wilder/EMA-TR/EWMA/GARCH): pass       |
//|         seed=true on the 1st (oldest) bar, false afterwards.     |
//|  6) band_k = center +/- mult_k * vol            (points)         |
//|     percent: center*(1 +/- mult_k*vol/center)                   |
//|     log/geometric: center*exp(+/- mult_k*vol/center)            |
//|     (GetVolatilityScaleSpace() tells you points/relative/log)   |
//|                                                                  |
//|  Minimum history per bar: use VolatilityLookbackRequirement()    |
//|  and VolatilityWarmupRequirement() to size the calculation.      |
//|  If VolatilityNeedsCenterWindow()==true, the midline must exist  |
//|  for the same bars (residual-based methods).                    |
//|                                                                  |
//|  NOTE: taken verbatim from a file that compiles; still COMPILE   |
//|  it in your project and rename identifiers if any collide.       |
//|  Nothing here was validated on your specific host.               |
//+------------------------------------------------------------------+
#ifndef __ROBUST_VOLATILITY_CORE_MQH__
#define __ROBUST_VOLATILITY_CORE_MQH__

#define EPSILON_VALUE 1.0e-12
#define PI_VALUE      3.14159265358979323846

//==================== ENUMS ====================
enum ENUM_VOLATILITY_METHOD
  {
   VOL_WILDER_ATR = 0,                // Wilder ATR (RMA of True Range)
   VOL_SMA_TR,                        // SMA of True Range
   VOL_EMA_TR,                        // EMA of True Range
   VOL_MEDIAN_TR,                     // Median True Range
   VOL_UPPER_TRIMMED_TR,              // Upper-tail trimmed True Range
   VOL_SYMMETRIC_TRIMMED_TR,          // Symmetric trimmed True Range
   VOL_UPPER_WINSORIZED_TR,           // Upper-tail winsorized True Range
   VOL_SYMMETRIC_WINSORIZED_TR,       // Symmetric winsorized True Range
   VOL_HUBER_TR,                      // Huber M-location of True Range
   VOL_TUKEY_TR,                      // Tukey biweight location of True Range
   VOL_NATR_WILDER,                   // Wilder-smoothed normalized True Range
   VOL_LOG_TR_WILDER,                 // Wilder-smoothed logarithmic True Range
   VOL_MEDIAN_LOG_TR,                 // Median logarithmic True Range
   VOL_UPPER_WINSORIZED_LOG_TR,       // Upper-tail winsorized logarithmic True Range
   VOL_HUBER_LOG_TR,                  // Huber M-location of logarithmic True Range
   VOL_TUKEY_LOG_TR,                  // Tukey biweight logarithmic True Range
   VOL_FAST_SLOW_MAX_LOG_TR,          // Max(fast EMA log-TR, slow robust log-TR)
   VOL_RESIDUAL_RMS,                  // Residual RMS around selected center
   VOL_EWMA_RESIDUAL_STD,             // EWMA standard deviation of residuals
   VOL_HUBER_EWMA_RESIDUAL_STD,       // Huber-clipped EWMA residual deviation
   VOL_MAD_RESIDUAL,                  // MAD scale of residuals
   VOL_HUBER_SCALE_RESIDUAL,          // Huber M-scale of residuals
   VOL_TUKEY_SCALE_RESIDUAL,          // Tukey biweight scale of residuals
   VOL_QN_RESIDUAL,                   // Rousseeuw-Croux Qn residual scale
   VOL_SN_RESIDUAL,                   // Rousseeuw-Croux Sn residual scale
   VOL_ABS_RESIDUAL_QUANTILE,         // Quantile of absolute residuals
   VOL_CLOSE_TO_CLOSE,                // Rolling close-to-close log volatility
   VOL_EWMA_LOG_RETURN,               // EWMA log-return volatility
   VOL_GARCH_11,                      // GARCH(1,1) log-return volatility
   VOL_PARKINSON,                     // Parkinson high-low volatility
   VOL_GARMAN_KLASS,                  // Garman-Klass OHLC volatility
   VOL_ROGERS_SATCHELL,               // Rogers-Satchell OHLC volatility
   VOL_YANG_ZHANG,                    // Yang-Zhang OHLC volatility
   VOL_IQR_RESIDUAL                   // IQR of residuals / 1.349 (~sigma) [new v4.3]
  };

enum ENUM_APPLIED_PRICE_ADVANCED
  {
   APPLIED_CLOSE = 0,                 // Close
   APPLIED_OPEN,                      // Open
   APPLIED_HIGH,                      // High
   APPLIED_LOW,                       // Low
   APPLIED_MEDIAN,                    // Median price (H+L)/2
   APPLIED_TYPICAL,                   // Typical price (H+L+C)/3
   APPLIED_WEIGHTED                   // Weighted price (H+L+2C)/4
  };

enum ENUM_RESIDUAL_SPACE
  {
   RESIDUAL_IN_POINTS = 0,            // Applied price - center
   RESIDUAL_IN_LOG_RATIO              // log(applied price / center)
  };

enum ENUM_SCALE_SPACE
  {
   SCALE_POINTS = 0,
   SCALE_RELATIVE,
   SCALE_LOG
  };

//==================== MODULE PARAMETERS ====================
// (the host exposes these on #include; rename if they collide)
input ENUM_APPLIED_PRICE_ADVANCED AppliedPriceMode = APPLIED_CLOSE;  // Preco aplicado (usado por AppliedPriceAt)
input ENUM_VOLATILITY_METHOD VolatilityMethod = VOL_QN_RESIDUAL;     // VOLATILITY CORE (the main selector)
input int    VolatilityPeriod        = 100;   // Core lookback (= "ATR period")
input ENUM_RESIDUAL_SPACE ResidualSpace = RESIDUAL_IN_LOG_RATIO;     // Residual space (linear vs log)
input int    FastVolatilityPeriod    = 24;    // Used by the Max(fast/slow log-TR) method
input double ResidualQuantile        = 0.95;  // "quantile of absolute residuals" method
input double EwmaLambda              = 0.94;  // EWMA (RiskMetrics ~0.94)
input double GarchAlpha              = 0.05;  // GARCH(1,1) alpha
input double GarchBeta               = 0.90;  // GARCH(1,1) beta (alpha+beta<1)
input double SymmetricTrimFraction   = 0.10;  // Trim fraction (trimmed)
input double UpperTailFraction       = 0.05;  // Upper-tail fraction
input double TukeyConstant           = 6.0;   // Tukey biweight tuning constant
input double HuberConstant           = 1.345; // Huber tuning constant
input int    PairwiseScaleMaximumWindow = 200;// Window cap for Qn/Sn (O(n^2) cost)
input int    RecursiveWarmupBars     = 300;   // Warmup for recursive methods

// Speed/accuracy mode for the pairwise estimators (Qn/Sn).
enum ENUM_PAIRWISE_MODE
  {
   PAIRWISE_EXACT = 0,      // Exact: full window every bar (slowest)
   PAIRWISE_BALANCED,       // Balanced: exact on recent bars, sampled history
   PAIRWISE_FAST,           // Fast: heavily sampled history (lightest)
   PAIRWISE_MANUAL          // Manual: uses the fields below
  };
input ENUM_PAIRWISE_MODE PairwiseMode = PAIRWISE_BALANCED;
input int    PairwiseHistoricalSampleSize = 0; // (MANUAL mode) 0 = full window
input int    PairwiseHistoricalStride     = 1; // (MANUAL mode) 1 = every bar

//==================== MODULE INTERNAL STATE ====================
double g_selectScratch[];
double g_deviationScratch[];
double g_pairwiseSampleScratch[];
double g_qnDistanceScratch[];
double g_windowScratch[];       // window scratch (auto-resizes)
double g_residual[];            // (price - midline) cache; filled by PrepareResidualData
double g_volState1[];           // recursive volatility state (Wilder/EMA-TR/EWMA/GARCH)
double g_volState2[];
double g_volState3[];
int    g_effPairwiseSample = 0;
int    g_effPairwiseStride = 1;

// Resize the bar-count-dependent internal state. Call whenever the size
// changes, BEFORE filling the data and calling ComputeVolatilityAt.
void RSE_ResizeVolatility(const int size)
  {
   ArrayResize(g_residual,size);
   ArrayResize(g_volState1,size);
   ArrayResize(g_volState2,size);
   ArrayResize(g_volState3,size);
   // g_windowScratch and the Qn/Sn scratch auto-resize via EnsureScratchSize.
  }

//========= HOST-PROVIDED DATA (fill before calling) =========
// Index 0 = newest bar (series indexing). The host sizes/fills these.
double g_sourceOpen[];
double g_sourceHigh[];
double g_sourceLow[];
double g_sourceClose[];
int    g_sourceRates = 0;
double g_price[];          // applied price per bar
double g_trueRange[];      // True Range per bar
double g_logTrueRange[];   // ln(max(H,Cprev)/min(L,Cprev)) per bar
double g_center[];         // the host MIDLINE (Starc average) per bar

// Resolve the effective Qn/Sn sampling/stride values (call in OnInit).
void ResolvePairwiseMode()
  {
   switch(PairwiseMode)
     {
      case PAIRWISE_EXACT:    g_effPairwiseSample=0;  g_effPairwiseStride=1; break;
      case PAIRWISE_BALANCED: g_effPairwiseSample=80; g_effPairwiseStride=2; break;
      case PAIRWISE_FAST:     g_effPairwiseSample=50; g_effPairwiseStride=3; break;
      default:                g_effPairwiseSample=MathMax(0,PairwiseHistoricalSampleSize);
                              g_effPairwiseStride=MathMax(1,PairwiseHistoricalStride); break;
     }
  }

//==================== HELPERS (basic/statistical/window) ====================

// BASIC HELPERS
//====================================================================

double ClampValue(const double value,const double minimum,const double maximum)
  {
   if(value<minimum) return minimum;
   if(value>maximum) return maximum;
   return value;
  }

bool IsFiniteValue(const double value)
  {
   return MathIsValidNumber(value);
  }

int SafePeriod(const int value,const int minimum)
  {
   if(value<minimum) return minimum;
   return value;
  }

void ResizeInternalArray(double &array[],const int size)
  {
   ArrayResize(array,size);
   // Internal arrays deliberately use normal indexing. We explicitly treat index 0 as newest.
   ArraySetAsSeries(array,false);
  }

void InitializeStateArray(double &array[])
  {
   ArrayInitialize(array,EMPTY_VALUE);
  }

void ShiftSeriesArrayOneBar(double &array[],const int lastIndex)
  {
   int size=ArraySize(array);
   if(size<=1) return;
   int last=MathMin(lastIndex,size-1);
   for(int i=last;i>=1;i--) array[i]=array[i-1];
   array[0]=EMPTY_VALUE;
  }

void SortAscending(double &array[],const int count)
  {
   if(count>1) ArraySort(array,count,0,MODE_ASCEND);
  }

void EnsureScratchSize(double &array[],const int requiredSize)
  {
   if(ArraySize(array)<requiredSize) ArrayResize(array,requiredSize);
  }

bool CopyWindow(double &source[],const int shift,const int length,double &window[])
  {
   if(length<=0 || shift<0 || shift+length>ArraySize(source)) return false;
   ArrayResize(window,length);
   for(int i=0;i<length;i++) window[i]=source[shift+i];
   return true;
  }

//====================================================================
// STATISTICAL HELPERS
//====================================================================

double MeanArray(double &values[],const int count)
  {
   if(count<=0) return 0.0;
   double sum=0.0;
   for(int i=0;i<count;i++) sum+=values[i];
   return sum/count;
  }

double MeanAt(double &source[],const int shift,const int length)
  {
   if(length<=0 || shift+length>ArraySize(source)) return EMPTY_VALUE;
   double sum=0.0;
   for(int i=0;i<length;i++) sum+=source[shift+i];
   return sum/length;
  }

double LinearWeightedMeanAt(double &source[],const int shift,const int length)
  {
   if(length<=0 || shift+length>ArraySize(source)) return EMPTY_VALUE;
   double weightedSum=0.0;
   double weightSum=0.0;
   for(int i=0;i<length;i++)
     {
      double weight=(double)(length-i);
      weightedSum+=weight*source[shift+i];
      weightSum+=weight;
     }
   if(weightSum<=0.0) return source[shift];
   return weightedSum/weightSum;
  }

double MedianArray(double &values[],const int count)
  {
   if(count<=0) return 0.0;
   EnsureScratchSize(g_selectScratch,count);
   for(int i=0;i<count;i++) g_selectScratch[i]=values[i];
   ArraySort(g_selectScratch,count,0,MODE_ASCEND);
   if((count%2)==1) return g_selectScratch[count/2];
   return 0.5*(g_selectScratch[count/2-1]+g_selectScratch[count/2]);
  }


double WeightedMedianAt(double &source[],const int shift,const int length)
  {
   if(length<=0 || shift<0 || shift+length>ArraySize(source)) return EMPTY_VALUE;
   double values[];
   double weights[];
   ArrayResize(values,length);
   ArrayResize(weights,length);
   for(int i=0;i<length;i++)
     {
      values[i]=source[shift+i];
      weights[i]=(double)(length-i); // newest observation receives the largest weight
     }

   for(int i=1;i<length;i++)
     {
      double valueKey=values[i];
      double weightKey=weights[i];
      int j=i-1;
      while(j>=0 && values[j]>valueKey)
        {
         values[j+1]=values[j];
         weights[j+1]=weights[j];
         j--;
        }
      values[j+1]=valueKey;
      weights[j+1]=weightKey;
     }

   double totalWeight=0.0;
   for(int i=0;i<length;i++) totalWeight+=weights[i];
   double threshold=0.5*totalWeight;
   double cumulative=0.0;
   for(int i=0;i<length;i++)
     {
      cumulative+=weights[i];
      if(cumulative>=threshold) return values[i];
     }
   return values[length-1];
  }

double QuantileArray(double &values[],const int count,double probability)
  {
   if(count<=0) return 0.0;
   probability=ClampValue(probability,0.0,1.0);
   double sorted[];
   ArrayResize(sorted,count);
   for(int i=0;i<count;i++) sorted[i]=values[i];
   SortAscending(sorted,count);
   if(count==1) return sorted[0];
   double position=probability*(count-1);
   int lower=(int)MathFloor(position);
   int upper=(int)MathCeil(position);
   if(lower==upper) return sorted[lower];
   double fraction=position-lower;
   return sorted[lower]+fraction*(sorted[upper]-sorted[lower]);
  }

double SymmetricTrimmedMeanArray(double &values[],const int count,double fraction)
  {
   if(count<=0) return 0.0;
   fraction=ClampValue(fraction,0.0,0.49);
   double sorted[];
   ArrayResize(sorted,count);
   for(int i=0;i<count;i++) sorted[i]=values[i];
   SortAscending(sorted,count);
   int cut=(int)MathFloor(count*fraction);
   if(2*cut>=count) return MedianArray(values,count);
   double sum=0.0;
   int used=0;
   for(int i=cut;i<count-cut;i++)
     {
      sum+=sorted[i];
      used++;
     }
   if(used<=0) return MedianArray(values,count);
   return sum/used;
  }

double UpperTrimmedMeanArray(double &values[],const int count,double fraction)
  {
   if(count<=0) return 0.0;
   fraction=ClampValue(fraction,0.0,0.49);
   double sorted[];
   ArrayResize(sorted,count);
   for(int i=0;i<count;i++) sorted[i]=values[i];
   SortAscending(sorted,count);
   int cut=(int)MathFloor(count*fraction);
   int used=count-cut;
   if(used<=0) return MedianArray(values,count);
   double sum=0.0;
   for(int i=0;i<used;i++) sum+=sorted[i];
   return sum/used;
  }

double SymmetricWinsorizedMeanArray(double &values[],const int count,double fraction)
  {
   if(count<=0) return 0.0;
   fraction=ClampValue(fraction,0.0,0.49);
   double sorted[];
   ArrayResize(sorted,count);
   for(int i=0;i<count;i++) sorted[i]=values[i];
   SortAscending(sorted,count);
   int cut=(int)MathFloor(count*fraction);
   if(2*cut>=count) return MedianArray(values,count);
   double lowLimit=sorted[cut];
   double highLimit=sorted[count-1-cut];
   double sum=0.0;
   for(int i=0;i<count;i++) sum+=ClampValue(sorted[i],lowLimit,highLimit);
   return sum/count;
  }

double UpperWinsorizedMeanArray(double &values[],const int count,double fraction)
  {
   if(count<=0) return 0.0;
   fraction=ClampValue(fraction,0.0,0.49);
   double sorted[];
   ArrayResize(sorted,count);
   for(int i=0;i<count;i++) sorted[i]=values[i];
   SortAscending(sorted,count);
   int cut=(int)MathFloor(count*fraction);
   int capIndex=count-1-cut;
   if(capIndex<0) capIndex=0;
   double highLimit=sorted[capIndex];
   double sum=0.0;
   for(int i=0;i<count;i++) sum+=MathMin(sorted[i],highLimit);
   return sum/count;
  }

double MadScaleArray(double &values[],const int count)
  {
   if(count<=1) return 0.0;
   double location=MedianArray(values,count);
   EnsureScratchSize(g_deviationScratch,count);
   for(int i=0;i<count;i++) g_deviationScratch[i]=MathAbs(values[i]-location);
   return 1.482602218505602*MedianArray(g_deviationScratch,count);
  }

double HuberLocationArray(double &values[],const int count,double tuningConstant)
  {
   if(count<=0) return 0.0;
   tuningConstant=MathMax(tuningConstant,0.10);
   double location=MedianArray(values,count);
   double scale=MadScaleArray(values,count);
   if(scale<=EPSILON_VALUE) return location;

   int iterations=MathMax(1,RobustMaximumIterations);
   double tolerance=MathMax(RobustConvergenceTolerance,1.0e-14);
   for(int iteration=0;iteration<iterations;iteration++)
     {
      double numerator=0.0;
      double denominator=0.0;
      for(int i=0;i<count;i++)
        {
         double residual=(values[i]-location)/scale;
         double absoluteResidual=MathAbs(residual);
         double weight=1.0;
         if(absoluteResidual>tuningConstant) weight=tuningConstant/absoluteResidual;
         numerator+=weight*values[i];
         denominator+=weight;
        }
      if(denominator<=EPSILON_VALUE) break;
      double nextLocation=numerator/denominator;
      if(MathAbs(nextLocation-location)<=tolerance*(MathAbs(location)+1.0))
        {
         location=nextLocation;
         break;
        }
      location=nextLocation;
     }
   return location;
  }

double TukeyLocationArray(double &values[],const int count,double tuningConstant)
  {
   if(count<=0) return 0.0;
   tuningConstant=MathMax(tuningConstant,0.10);
   double location=MedianArray(values,count);
   double scale=MadScaleArray(values,count);
   if(scale<=EPSILON_VALUE) return location;

   int iterations=MathMax(1,RobustMaximumIterations);
   double tolerance=MathMax(RobustConvergenceTolerance,1.0e-14);
   for(int iteration=0;iteration<iterations;iteration++)
     {
      double numerator=0.0;
      double denominator=0.0;
      for(int i=0;i<count;i++)
        {
         double u=(values[i]-location)/(tuningConstant*scale);
         if(MathAbs(u)<1.0)
           {
            double oneMinus=1.0-u*u;
            double weight=oneMinus*oneMinus;
            numerator+=weight*values[i];
            denominator+=weight;
           }
        }
      if(denominator<=EPSILON_VALUE) break;
      double nextLocation=numerator/denominator;
      if(MathAbs(nextLocation-location)<=tolerance*(MathAbs(location)+1.0))
        {
         location=nextLocation;
         break;
        }
      location=nextLocation;
     }
   return location;
  }


// Standard normal CDF approximation used only for Huber scale normalization.
double NormalCdf(const double value)
  {
   double x=value;
   double sign=1.0;
   if(x<0.0) { sign=-1.0; x=-x; }
   double t=1.0/(1.0+0.2316419*x);
   double density=0.3989422804014327*MathExp(-0.5*x*x);
   double polynomial=t*(0.319381530+t*(-0.356563782+t*(1.781477937+t*(-1.821255978+t*1.330274429))));
   double cdf=1.0-density*polynomial;
   if(sign<0.0) cdf=1.0-cdf;
   return ClampValue(cdf,0.0,1.0);
  }

double HuberScaleArray(double &values[],const int count,double tuningConstant)
  {
   if(count<2) return 0.0;
   tuningConstant=MathMax(tuningConstant,0.20);
   double location=HuberLocationArray(values,count,tuningConstant);
   double scale=MadScaleArray(values,count);
   if(scale<=EPSILON_VALUE) return 0.0;

   double phi=0.3989422804014327*MathExp(-0.5*tuningConstant*tuningConstant);
   double tail=1.0-NormalCdf(tuningConstant);
   double expectedClippedSquare=2.0*(NormalCdf(tuningConstant)-0.5-tuningConstant*phi)+
                                2.0*tuningConstant*tuningConstant*tail;
   expectedClippedSquare=MathMax(expectedClippedSquare,EPSILON_VALUE);

   int iterations=MathMax(1,RobustMaximumIterations);
   double tolerance=MathMax(RobustConvergenceTolerance,1.0e-14);
   for(int iteration=0;iteration<iterations;iteration++)
     {
      double cap=tuningConstant*scale;
      double sum=0.0;
      for(int i=0;i<count;i++)
        {
         double residual=values[i]-location;
         double clipped=ClampValue(residual,-cap,cap);
         sum+=clipped*clipped;
        }
      double nextScale=MathSqrt(MathMax(0.0,(sum/count)/expectedClippedSquare));
      if(MathAbs(nextScale-scale)<=tolerance*(scale+1.0))
        {
         scale=nextScale;
         break;
        }
      scale=nextScale;
      if(scale<=EPSILON_VALUE) break;
     }
   return MathMax(0.0,scale);
  }

double TukeyBiweightScaleArray(double &values[],const int count,double tuningConstant)
  {
   if(count<3) return 0.0;
   tuningConstant=MathMax(tuningConstant,1.0);
   double location=MedianArray(values,count);
   double mad=MadScaleArray(values,count);
   if(mad<=EPSILON_VALUE) return 0.0;

   double numerator=0.0;
   double denominator=0.0;
   int used=0;
   for(int i=0;i<count;i++)
     {
      double residual=values[i]-location;
      double u=residual/(tuningConstant*mad);
      if(MathAbs(u)<1.0)
        {
         double oneMinus=1.0-u*u;
         double oneMinus2=oneMinus*oneMinus;
         numerator+=residual*residual*oneMinus2*oneMinus2;
         denominator+=oneMinus*(1.0-5.0*u*u);
         used++;
        }
     }
   if(used<2 || MathAbs(denominator)<=EPSILON_VALUE) return mad;
   double variance=count*numerator/(denominator*denominator);
   return MathSqrt(MathMax(0.0,variance));
  }

// Finite-sample correction used by the asymptotically normalized Qn scale.
double QnFiniteSampleCorrection(const int count)
  {
   double correction=1.0;
   if(count<=9)
     {
      if(count==2) correction=0.399;
      else if(count==3) correction=0.994;
      else if(count==4) correction=0.512;
      else if(count==5) correction=0.844;
      else if(count==6) correction=0.611;
      else if(count==7) correction=0.857;
      else if(count==8) correction=0.669;
      else if(count==9) correction=0.872;
     }
   else
     {
      if((count%2)==1) correction=(double)count/(count+1.4);
      else             correction=(double)count/(count+3.8);
     }
   return correction;
  }

// Exact Rousseeuw-Croux Qn for a supplied sample. The live/recent path uses
// the full residual window; optional history sampling is applied elsewhere.
double QnScaleArray(double &values[],const int count)
  {
   if(count<2) return 0.0;
   int pairCount=count*(count-1)/2;
   EnsureScratchSize(g_qnDistanceScratch,pairCount);
   int position=0;
   for(int i=0;i<count-1;i++)
     {
      double left=values[i];
      for(int j=i+1;j<count;j++)
         g_qnDistanceScratch[position++]=MathAbs(left-values[j]);
     }

   ArraySort(g_qnDistanceScratch,pairCount,0,MODE_ASCEND);
   int h=count/2+1;
   int order=h*(h-1)/2; // one-based order statistic
   int index=MathMax(0,MathMin(pairCount-1,order-1));
   return 2.2219*QnFiniteSampleCorrection(count)*g_qnDistanceScratch[index];
  }

// Asymptotically normalized Rousseeuw-Croux Sn.
double SnScaleArray(double &values[],const int count)
  {
   if(count<2) return 0.0;
   double rowMedians[];
   ArrayResize(rowMedians,count);
   double distances[];
   ArrayResize(distances,count);
   for(int i=0;i<count;i++)
     {
      for(int j=0;j<count;j++) distances[j]=MathAbs(values[i]-values[j]);
      rowMedians[i]=MedianArray(distances,count);
     }
   return 1.1926*MedianArray(rowMedians,count);
  }

//====================================================================
// REGRESSION HELPERS
//====================================================================

bool OrdinaryLinearRegression(double &values[],const int count,double &intercept,double &slope)
  {
   if(count<2) return false;
   double sumX=0.0,sumY=0.0,sumXX=0.0,sumXY=0.0;
   for(int i=0;i<count;i++)
     {
      double x=-(double)i; // x=0 is the current edge
      double y=values[i];
      sumX+=x;
      sumY+=y;
      sumXX+=x*x;
      sumXY+=x*y;
     }
   double denominator=count*sumXX-sumX*sumX;
   if(MathAbs(denominator)<=EPSILON_VALUE)
     {
      intercept=MeanArray(values,count);
      slope=0.0;
      return false;
     }
   slope=(count*sumXY-sumX*sumY)/denominator;
   intercept=(sumY-slope*sumX)/count;
   return true;
  }

double RobustLinearRegressionEdge(double &values[],const int count,const bool useTukey)
  {
   if(count<3) return MeanArray(values,count);
   double intercept=0.0,slope=0.0;
   OrdinaryLinearRegression(values,count,intercept,slope);

   double residuals[];
   ArrayResize(residuals,count);
   int iterations=MathMax(1,RobustMaximumIterations);
   double tolerance=MathMax(RobustConvergenceTolerance,1.0e-14);
   double tuning=useTukey?MathMax(TukeyConstant,0.10):MathMax(HuberConstant,0.10);

   for(int iteration=0;iteration<iterations;iteration++)
     {
      for(int i=0;i<count;i++)
        {
         double x=-(double)i;
         residuals[i]=values[i]-(intercept+slope*x);
        }
      double scale=MadScaleArray(residuals,count);
      if(scale<=EPSILON_VALUE) break;

      double sw=0.0,swx=0.0,swy=0.0,swxx=0.0,swxy=0.0;
      for(int i=0;i<count;i++)
        {
         double x=-(double)i;
         double standardized=residuals[i]/scale;
         double absoluteValue=MathAbs(standardized);
         double weight=1.0;
         if(useTukey)
           {
            double u=standardized/tuning;
            if(MathAbs(u)>=1.0) weight=0.0;
            else
              {
               double oneMinus=1.0-u*u;
               weight=oneMinus*oneMinus;
              }
           }
         else if(absoluteValue>tuning)
            weight=tuning/absoluteValue;

         sw+=weight;
         swx+=weight*x;
         swy+=weight*values[i];
         swxx+=weight*x*x;
         swxy+=weight*x*values[i];
        }

      double denominator=sw*swxx-swx*swx;
      if(sw<=EPSILON_VALUE || MathAbs(denominator)<=EPSILON_VALUE) break;
      double nextSlope=(sw*swxy-swx*swy)/denominator;
      double nextIntercept=(swy-nextSlope*swx)/sw;
      double change=MathMax(MathAbs(nextIntercept-intercept),MathAbs(nextSlope-slope));
      intercept=nextIntercept;
      slope=nextSlope;
      if(change<=tolerance*(MathAbs(intercept)+MathAbs(slope)+1.0)) break;
     }
   return intercept;
  }

bool SolveAugmentedSystem4(double &matrix[][5],const int size,double &solution[])
  {
   ArrayResize(solution,size);
   for(int column=0;column<size;column++)
     {
      int pivot=column;
      double largest=MathAbs(matrix[column][column]);
      for(int row=column+1;row<size;row++)
        {
         double candidate=MathAbs(matrix[row][column]);
         if(candidate>largest)
           {
            largest=candidate;
            pivot=row;
           }
        }
      if(largest<=EPSILON_VALUE) return false;
      if(pivot!=column)
        {
         for(int j=column;j<=size;j++)
           {
            double temporary=matrix[column][j];
            matrix[column][j]=matrix[pivot][j];
            matrix[pivot][j]=temporary;
           }
        }

      double divisor=matrix[column][column];
      for(int j=column;j<=size;j++) matrix[column][j]/=divisor;

      for(int row=0;row<size;row++)
        {
         if(row==column) continue;
         double factor=matrix[row][column];
         for(int j=column;j<=size;j++) matrix[row][j]-=factor*matrix[column][j];
        }
     }
   for(int i=0;i<size;i++) solution[i]=matrix[i][size];
   return true;
  }

double CausalSavitzkyGolayEdge(double &values[],const int count,int order)
  {
   order=MathMax(1,MathMin(3,order));
   if(count<order+2) return MeanArray(values,count);
   int size=order+1;
   double matrix[4][5];
   for(int r=0;r<4;r++) for(int c=0;c<5;c++) matrix[r][c]=0.0;

   for(int row=0;row<size;row++)
     {
      for(int column=0;column<size;column++)
        {
         double sum=0.0;
         int power=row+column;
         for(int i=0;i<count;i++) sum+=MathPow(-(double)i,power);
         matrix[row][column]=sum;
        }
      double rightSide=0.0;
      for(int i=0;i<count;i++) rightSide+=values[i]*MathPow(-(double)i,row);
      matrix[row][size]=rightSide;
     }

   double coefficients[];
   if(!SolveAugmentedSystem4(matrix,size,coefficients)) return MeanArray(values,count);
   return coefficients[0]; // fitted value at x=0
  }

//====================================================================
// PRICE AND WINDOW HELPERS
//====================================================================

double AppliedPriceAt(const int shift)
  {
   switch(AppliedPriceMode)
     {
      case APPLIED_CLOSE:    return g_sourceClose[shift];
      case APPLIED_OPEN:     return g_sourceOpen[shift];
      case APPLIED_HIGH:     return g_sourceHigh[shift];
      case APPLIED_LOW:      return g_sourceLow[shift];
      case APPLIED_MEDIAN:   return 0.5*(g_sourceHigh[shift]+g_sourceLow[shift]);
      case APPLIED_TYPICAL:  return (g_sourceHigh[shift]+g_sourceLow[shift]+g_sourceClose[shift])/3.0;
      case APPLIED_WEIGHTED: return (g_sourceHigh[shift]+g_sourceLow[shift]+2.0*g_sourceClose[shift])/4.0;
     }
   return g_sourceClose[shift];
  }

double WindowMedianAt(double &source[],const int shift,const int length)
  {
   if(length<=0 || shift<0 || shift+length>ArraySize(source)) return EMPTY_VALUE;
   EnsureScratchSize(g_windowScratch,length);
   for(int i=0;i<length;i++) g_windowScratch[i]=source[shift+i];
   return MedianArray(g_windowScratch,length);
  }

double WindowSymmetricTrimmedAt(double &source[],const int shift,const int length,const double fraction)
  {
   double window[];
   if(!CopyWindow(source,shift,length,window)) return EMPTY_VALUE;
   return SymmetricTrimmedMeanArray(window,length,fraction);
  }

double WindowUpperTrimmedAt(double &source[],const int shift,const int length,const double fraction)
  {
   double window[];
   if(!CopyWindow(source,shift,length,window)) return EMPTY_VALUE;
   return UpperTrimmedMeanArray(window,length,fraction);
  }

double WindowSymmetricWinsorizedAt(double &source[],const int shift,const int length,const double fraction)
  {
   double window[];
   if(!CopyWindow(source,shift,length,window)) return EMPTY_VALUE;
   return SymmetricWinsorizedMeanArray(window,length,fraction);
  }

double WindowUpperWinsorizedAt(double &source[],const int shift,const int length,const double fraction)
  {
   double window[];
   if(!CopyWindow(source,shift,length,window)) return EMPTY_VALUE;
   return UpperWinsorizedMeanArray(window,length,fraction);
  }

double WindowHuberAt(double &source[],const int shift,const int length)
  {
   if(length<=0 || shift<0 || shift+length>ArraySize(source)) return EMPTY_VALUE;
   EnsureScratchSize(g_windowScratch,length);
   for(int i=0;i<length;i++) g_windowScratch[i]=source[shift+i];
   return HuberLocationArray(g_windowScratch,length,HuberConstant);
  }

double WindowTukeyAt(double &source[],const int shift,const int length)
  {
   if(length<=0 || shift<0 || shift+length>ArraySize(source)) return EMPTY_VALUE;
   EnsureScratchSize(g_windowScratch,length);
   for(int i=0;i<length;i++) g_windowScratch[i]=source[shift+i];
   return TukeyLocationArray(g_windowScratch,length,TukeyConstant);
  }

double WindowQuantileAt(double &source[],const int shift,const int length,const double probability)
  {
   double window[];
   if(!CopyWindow(source,shift,length,window)) return EMPTY_VALUE;
   return QuantileArray(window,length,probability);
  }

double WindowMinimum(double &source[],const int shift,const int length)
  {
   if(length<=0 || shift+length>ArraySize(source)) return EMPTY_VALUE;
   double result=source[shift];
   for(int i=1;i<length;i++) result=MathMin(result,source[shift+i]);
   return result;
  }

double WindowMaximum(double &source[],const int shift,const int length)
  {
   if(length<=0 || shift+length>ArraySize(source)) return EMPTY_VALUE;
   double result=source[shift];
   for(int i=1;i<length;i++) result=MathMax(result,source[shift+i]);
   return result;
  }


//==================== VOLATILITY CORE ====================
// VOLATILITY HELPERS
//====================================================================

void PrepareResidualData(const int maximumIndex)
  {
   int available=MathMin(ArraySize(g_residual),ArraySize(g_center));
   int last=MathMin(maximumIndex,available-1);
   for(int shift=0;shift<=last;shift++)
     {
      double center=g_center[shift];
      if(!IsFiniteValue(center) || center==EMPTY_VALUE)
        {
         g_residual[shift]=EMPTY_VALUE;
         continue;
        }
      if(ResidualSpace==RESIDUAL_IN_LOG_RATIO)
        {
         if(g_price[shift]>0.0 && center>0.0) g_residual[shift]=MathLog(g_price[shift]/center);
         else g_residual[shift]=0.0;
        }
      else g_residual[shift]=g_price[shift]-center;
     }
  }

double ResidualAt(const int shift)
  {
   if(shift>=0 && shift<ArraySize(g_residual))
     {
      double cached=g_residual[shift];
      if(IsFiniteValue(cached) && cached!=EMPTY_VALUE) return cached;
     }
   if(shift<0 || shift>=ArraySize(g_center)) return 0.0;
   double center=g_center[shift];
   if(!IsFiniteValue(center) || center==EMPTY_VALUE) return 0.0;
   if(ResidualSpace==RESIDUAL_IN_LOG_RATIO)
     {
      if(g_price[shift]>0.0 && center>0.0) return MathLog(g_price[shift]/center);
      return 0.0;
     }
   return g_price[shift]-center;
  }

double LogReturnAt(const int shift)
  {
   if(shift+1>=g_sourceRates) return 0.0;
   if(g_sourceClose[shift]>0.0 && g_sourceClose[shift+1]>0.0) return MathLog(g_sourceClose[shift]/g_sourceClose[shift+1]);
   return 0.0;
  }

bool BuildResidualWindow(const int shift,const int length,double &window[])
  {
   if(length<=0 || shift<0 || shift+length>ArraySize(g_residual)) return false;
   ArrayResize(window,length);
   for(int i=0;i<length;i++)
     {
      double value=g_residual[shift+i];
      if(!IsFiniteValue(value) || value==EMPTY_VALUE) return false;
      window[i]=value;
     }
   return true;
  }

bool BuildNormalizedTrWindow(const int shift,const int length,double &window[])
  {
   if(length<=0 || shift+length>ArraySize(g_center)) return false;
   ArrayResize(window,length);
   for(int i=0;i<length;i++)
     {
      double center=g_center[shift+i];
      if(!IsFiniteValue(center) || center==EMPTY_VALUE) return false;
      double denominator=MathAbs(center);
      window[i]=(denominator>EPSILON_VALUE)?g_trueRange[shift+i]/denominator:0.0;
     }
   return true;
  }

bool BuildLogTrWindow(const int shift,const int length,double &window[])
  {
   return CopyWindow(g_logTrueRange,shift,length,window);
  }

double ResidualRmsAt(const int shift,const int length)
  {
   double residuals[];
   if(!BuildResidualWindow(shift,length,residuals)) return 0.0;
   double sumSquares=0.0;
   for(int i=0;i<length;i++) sumSquares+=residuals[i]*residuals[i];
   int denominator=MathMax(1,length-1);
   return MathSqrt(MathMax(0.0,sumSquares/denominator));
  }

double ResidualMadAt(const int shift,const int length)
  {
   double residuals[];
   if(!BuildResidualWindow(shift,length,residuals)) return 0.0;
   return MadScaleArray(residuals,length);
  }

double ResidualHuberScaleAt(const int shift,const int length)
  {
   double residuals[];
   if(!BuildResidualWindow(shift,length,residuals)) return 0.0;
   return HuberScaleArray(residuals,length,HuberConstant);
  }

double ResidualTukeyScaleAt(const int shift,const int length)
  {
   double residuals[];
   if(!BuildResidualWindow(shift,length,residuals)) return 0.0;
   return TukeyBiweightScaleArray(residuals,length,TukeyConstant);
  }

// IQR of residuals, scaled to ~sigma under normality.
double ResidualIqrAt(const int shift,const int length)
  {
   double residuals[];
   if(!BuildResidualWindow(shift,length,residuals)) return 0.0;
   double q1=QuantileArray(residuals,length,0.25);
   double q3=QuantileArray(residuals,length,0.75);
   double iqr=q3-q1;
   // 1.349 = 2 * 0.6745 (IQR ~= 1.349*sigma for the normal). Known constant.
   return MathMax(0.0,iqr/1.349);
  }

double ResidualQnAt(const int shift,const int requestedLength)
  {
   int maximum=SafePeriod(PairwiseScaleMaximumWindow,10);
   int fullLength=MathMin(requestedLength,maximum);
   if(fullLength<2 || shift<0 || shift+fullLength>ArraySize(g_residual)) return 0.0;

   int sampleCount=fullLength;
   int requestedSample=g_effPairwiseSample;
   if(shift>MathMax(0,PairwiseExactRecentBars) && requestedSample>=10 && requestedSample<fullLength)
      sampleCount=requestedSample;

   EnsureScratchSize(g_pairwiseSampleScratch,sampleCount);
   if(sampleCount==fullLength)
     {
      for(int i=0;i<sampleCount;i++)
        {
         double value=g_residual[shift+i];
         if(!IsFiniteValue(value) || value==EMPTY_VALUE) return 0.0;
         g_pairwiseSampleScratch[i]=value;
        }
     }
   else
     {
      double step=(double)(fullLength-1)/(double)(sampleCount-1);
      for(int i=0;i<sampleCount;i++)
        {
         int offset=(int)MathRound(i*step);
         offset=MathMax(0,MathMin(fullLength-1,offset));
         double value=g_residual[shift+offset];
         if(!IsFiniteValue(value) || value==EMPTY_VALUE) return 0.0;
         g_pairwiseSampleScratch[i]=value;
        }
     }

   int pairCount=sampleCount*(sampleCount-1)/2;
   EnsureScratchSize(g_qnDistanceScratch,pairCount);
   int position=0;
   for(int i=0;i<sampleCount-1;i++)
     {
      double left=g_pairwiseSampleScratch[i];
      for(int j=i+1;j<sampleCount;j++)
         g_qnDistanceScratch[position++]=MathAbs(left-g_pairwiseSampleScratch[j]);
     }

   ArraySort(g_qnDistanceScratch,pairCount,0,MODE_ASCEND);
   int h=sampleCount/2+1;
   int order=h*(h-1)/2;
   int index=MathMax(0,MathMin(pairCount-1,order-1));
   return 2.2219*QnFiniteSampleCorrection(sampleCount)*g_qnDistanceScratch[index];
  }

double ResidualSnAt(const int shift,const int requestedLength)
  {
   int maximum=SafePeriod(PairwiseScaleMaximumWindow,10);
   int length=MathMin(requestedLength,maximum);
   double residuals[];
   if(!BuildResidualWindow(shift,length,residuals)) return 0.0;
   return SnScaleArray(residuals,length);
  }

double AbsoluteResidualQuantileAt(const int shift,const int length)
  {
   double residuals[];
   if(!BuildResidualWindow(shift,length,residuals)) return 0.0;
   for(int i=0;i<length;i++) residuals[i]=MathAbs(residuals[i]);
   return QuantileArray(residuals,length,ClampValue(ResidualQuantile,0.50,0.9999));
  }

double CloseToCloseVolatilityAt(const int shift,const int length)
  {
   if(length<2 || shift+length>=g_sourceRates) return 0.0;
   double sum=0.0,sumSquares=0.0;
   for(int i=0;i<length;i++)
     {
      double value=LogReturnAt(shift+i);
      sum+=value;
      sumSquares+=value*value;
     }
   double mean=sum/length;
   double variance=(sumSquares-length*mean*mean)/MathMax(1,length-1);
   return MathSqrt(MathMax(0.0,variance));
  }

double ParkinsonVolatilityAt(const int shift,const int length)
  {
   if(length<=0 || shift+length>g_sourceRates) return 0.0;
   double sum=0.0;
   double divisor=4.0*MathLog(2.0);
   int used=0;
   for(int i=0;i<length;i++)
     {
      int bar=shift+i;
      if(g_sourceHigh[bar]<=0.0 || g_sourceLow[bar]<=0.0) continue;
      double range=MathLog(g_sourceHigh[bar]/g_sourceLow[bar]);
      sum+=range*range/divisor;
      used++;
     }
   if(used<=0) return 0.0;
   return MathSqrt(MathMax(0.0,sum/used));
  }

double GarmanKlassVolatilityAt(const int shift,const int length)
  {
   if(length<=0 || shift+length>g_sourceRates) return 0.0;
   double sum=0.0;
   double coefficient=2.0*MathLog(2.0)-1.0;
   int used=0;
   for(int i=0;i<length;i++)
     {
      int bar=shift+i;
      if(g_sourceHigh[bar]<=0.0 || g_sourceLow[bar]<=0.0 || g_sourceOpen[bar]<=0.0 || g_sourceClose[bar]<=0.0) continue;
      double highLow=MathLog(g_sourceHigh[bar]/g_sourceLow[bar]);
      double closeOpen=MathLog(g_sourceClose[bar]/g_sourceOpen[bar]);
      sum+=0.5*highLow*highLow-coefficient*closeOpen*closeOpen;
      used++;
     }
   if(used<=0) return 0.0;
   return MathSqrt(MathMax(0.0,sum/used));
  }

double RogersSatchellVarianceAtBar(const int bar)
  {
   if(g_sourceHigh[bar]<=0.0 || g_sourceLow[bar]<=0.0 || g_sourceOpen[bar]<=0.0 || g_sourceClose[bar]<=0.0) return 0.0;
   return MathLog(g_sourceHigh[bar]/g_sourceClose[bar])*MathLog(g_sourceHigh[bar]/g_sourceOpen[bar])+
          MathLog(g_sourceLow[bar]/g_sourceClose[bar])*MathLog(g_sourceLow[bar]/g_sourceOpen[bar]);
  }

double RogersSatchellVolatilityAt(const int shift,const int length)
  {
   if(length<=0 || shift+length>g_sourceRates) return 0.0;
   double sum=0.0;
   int used=0;
   for(int i=0;i<length;i++)
     {
      double value=RogersSatchellVarianceAtBar(shift+i);
      if(IsFiniteValue(value))
        {
         sum+=value;
         used++;
        }
     }
   if(used<=0) return 0.0;
   return MathSqrt(MathMax(0.0,sum/used));
  }

double YangZhangVolatilityAt(const int shift,const int length)
  {
   if(length<2 || shift+length>=g_sourceRates) return 0.0;
   double sumOvernight=0.0,sumOvernightSquares=0.0;
   double sumOpenClose=0.0,sumOpenCloseSquares=0.0;
   double sumRs=0.0;
   int used=0;

   for(int i=0;i<length;i++)
     {
      int bar=shift+i;
      if(g_sourceOpen[bar]<=0.0 || g_sourceHigh[bar]<=0.0 || g_sourceLow[bar]<=0.0 || g_sourceClose[bar]<=0.0 || g_sourceClose[bar+1]<=0.0) continue;
      double overnight=MathLog(g_sourceOpen[bar]/g_sourceClose[bar+1]);
      double openClose=MathLog(g_sourceClose[bar]/g_sourceOpen[bar]);
      double rs=RogersSatchellVarianceAtBar(bar);
      sumOvernight+=overnight;
      sumOvernightSquares+=overnight*overnight;
      sumOpenClose+=openClose;
      sumOpenCloseSquares+=openClose*openClose;
      sumRs+=rs;
      used++;
     }

   if(used<2) return 0.0;
   double meanOvernight=sumOvernight/used;
   double meanOpenClose=sumOpenClose/used;
   double overnightVariance=(sumOvernightSquares-used*meanOvernight*meanOvernight)/(used-1);
   double openCloseVariance=(sumOpenCloseSquares-used*meanOpenClose*meanOpenClose)/(used-1);
   double rsVariance=sumRs/used;
   double k=0.34/(1.34+(double)(used+1)/(double)(used-1));
   double variance=overnightVariance+k*openCloseVariance+(1.0-k)*rsVariance;
   return MathSqrt(MathMax(0.0,variance));
  }

double RollingLogReturnVarianceAt(const int shift,const int length)
  {
   double sigma=CloseToCloseVolatilityAt(shift,length);
   return sigma*sigma;
  }

ENUM_SCALE_SPACE GetVolatilityScaleSpace()
  {
   switch(VolatilityMethod)
     {
      case VOL_NATR_WILDER:
         return SCALE_RELATIVE;

      case VOL_LOG_TR_WILDER:
      case VOL_MEDIAN_LOG_TR:
      case VOL_UPPER_WINSORIZED_LOG_TR:
      case VOL_HUBER_LOG_TR:
      case VOL_TUKEY_LOG_TR:
      case VOL_FAST_SLOW_MAX_LOG_TR:
      case VOL_CLOSE_TO_CLOSE:
      case VOL_EWMA_LOG_RETURN:
      case VOL_GARCH_11:
      case VOL_PARKINSON:
      case VOL_GARMAN_KLASS:
      case VOL_ROGERS_SATCHELL:
      case VOL_YANG_ZHANG:
         return SCALE_LOG;

      case VOL_RESIDUAL_RMS:
      case VOL_EWMA_RESIDUAL_STD:
      case VOL_HUBER_EWMA_RESIDUAL_STD:
      case VOL_MAD_RESIDUAL:
      case VOL_HUBER_SCALE_RESIDUAL:
      case VOL_TUKEY_SCALE_RESIDUAL:
      case VOL_QN_RESIDUAL:
      case VOL_SN_RESIDUAL:
      case VOL_ABS_RESIDUAL_QUANTILE:
         return (ResidualSpace==RESIDUAL_IN_LOG_RATIO)?SCALE_LOG:SCALE_POINTS;
     }
   return SCALE_POINTS;
  }

double ComputeVolatilityAt(const int shift,const bool seed)
  {
   int length=SafePeriod(VolatilityPeriod,2);
   double result=0.0;

   switch(VolatilityMethod)
     {
      case VOL_WILDER_ATR:
        {
         if(seed) g_volState1[shift]=MeanAt(g_trueRange,shift,length);
         else
           {
            double previous=g_volState1[shift+1];
            if(!IsFiniteValue(previous) || previous==EMPTY_VALUE) previous=MeanAt(g_trueRange,shift+1,length);
            g_volState1[shift]=previous+(g_trueRange[shift]-previous)/length;
           }
         result=g_volState1[shift];
         break;
        }

      case VOL_SMA_TR:
         result=MeanAt(g_trueRange,shift,length);
         break;

      case VOL_EMA_TR:
        {
         double alpha=2.0/(length+1.0);
         if(seed) g_volState1[shift]=MeanAt(g_trueRange,shift,length);
         else
           {
            double previous=g_volState1[shift+1];
            if(!IsFiniteValue(previous) || previous==EMPTY_VALUE) previous=MeanAt(g_trueRange,shift+1,length);
            g_volState1[shift]=alpha*g_trueRange[shift]+(1.0-alpha)*previous;
           }
         result=g_volState1[shift];
         break;
        }

      case VOL_MEDIAN_TR:
         result=WindowMedianAt(g_trueRange,shift,length);
         break;

      case VOL_UPPER_TRIMMED_TR:
         result=WindowUpperTrimmedAt(g_trueRange,shift,length,UpperTailFraction);
         break;

      case VOL_SYMMETRIC_TRIMMED_TR:
         result=WindowSymmetricTrimmedAt(g_trueRange,shift,length,SymmetricTrimFraction);
         break;

      case VOL_UPPER_WINSORIZED_TR:
         result=WindowUpperWinsorizedAt(g_trueRange,shift,length,UpperTailFraction);
         break;

      case VOL_SYMMETRIC_WINSORIZED_TR:
         result=WindowSymmetricWinsorizedAt(g_trueRange,shift,length,SymmetricTrimFraction);
         break;

      case VOL_HUBER_TR:
         result=WindowHuberAt(g_trueRange,shift,length);
         break;

      case VOL_TUKEY_TR:
         result=WindowTukeyAt(g_trueRange,shift,length);
         break;

      case VOL_NATR_WILDER:
        {
         double normalizedWindow[];
         if(!BuildNormalizedTrWindow(shift,length,normalizedWindow)) break;
         double current=normalizedWindow[0];
         if(seed) g_volState1[shift]=MeanArray(normalizedWindow,length);
         else
           {
            double previous=g_volState1[shift+1];
            if(!IsFiniteValue(previous) || previous==EMPTY_VALUE) previous=MeanArray(normalizedWindow,length);
            g_volState1[shift]=previous+(current-previous)/length;
           }
         result=g_volState1[shift];
         break;
        }

      case VOL_LOG_TR_WILDER:
        {
         if(seed) g_volState1[shift]=MeanAt(g_logTrueRange,shift,length);
         else
           {
            double previous=g_volState1[shift+1];
            if(!IsFiniteValue(previous) || previous==EMPTY_VALUE) previous=MeanAt(g_logTrueRange,shift+1,length);
            g_volState1[shift]=previous+(g_logTrueRange[shift]-previous)/length;
           }
         result=g_volState1[shift];
         break;
        }

      case VOL_MEDIAN_LOG_TR:
         result=WindowMedianAt(g_logTrueRange,shift,length);
         break;

      case VOL_UPPER_WINSORIZED_LOG_TR:
         result=WindowUpperWinsorizedAt(g_logTrueRange,shift,length,UpperTailFraction);
         break;

      case VOL_HUBER_LOG_TR:
         result=WindowHuberAt(g_logTrueRange,shift,length);
         break;

      case VOL_TUKEY_LOG_TR:
         result=WindowTukeyAt(g_logTrueRange,shift,length);
         break;

      case VOL_FAST_SLOW_MAX_LOG_TR:
        {
         int fastLength=SafePeriod(FastVolatilityPeriod,2);
         double fastAlpha=2.0/(fastLength+1.0);
         if(seed) g_volState1[shift]=MeanAt(g_logTrueRange,shift,fastLength);
         else
           {
            double previousFast=g_volState1[shift+1];
            if(!IsFiniteValue(previousFast) || previousFast==EMPTY_VALUE) previousFast=MeanAt(g_logTrueRange,shift+1,fastLength);
            g_volState1[shift]=fastAlpha*g_logTrueRange[shift]+(1.0-fastAlpha)*previousFast;
           }
         double slow=WindowUpperWinsorizedAt(g_logTrueRange,shift,length,UpperTailFraction);
         result=MathMax(g_volState1[shift],slow);
         break;
        }

      case VOL_RESIDUAL_RMS:
         result=ResidualRmsAt(shift,length);
         break;

      case VOL_EWMA_RESIDUAL_STD:
        {
         double lambda=ClampValue(EwmaLambda,0.001,0.9999);
         double residual=ResidualAt(shift);
         if(seed)
           {
            double initial=ResidualRmsAt(shift,length);
            g_volState1[shift]=initial*initial;
           }
         else
           {
            double previousVariance=g_volState1[shift+1];
            if(!IsFiniteValue(previousVariance) || previousVariance==EMPTY_VALUE || previousVariance<0.0)
              {
               double initial=ResidualRmsAt(shift+1,length);
               previousVariance=initial*initial;
              }
            g_volState1[shift]=lambda*previousVariance+(1.0-lambda)*residual*residual;
           }
         result=MathSqrt(MathMax(0.0,g_volState1[shift]));
         break;
        }

      case VOL_HUBER_EWMA_RESIDUAL_STD:
        {
         double lambda=ClampValue(EwmaLambda,0.001,0.9999);
         double residual=ResidualAt(shift);
         if(seed)
           {
            double initial=ResidualMadAt(shift,length);
            if(initial<=EPSILON_VALUE) initial=ResidualRmsAt(shift,length);
            g_volState1[shift]=initial*initial;
           }
         else
           {
            double previousVariance=g_volState1[shift+1];
            if(!IsFiniteValue(previousVariance) || previousVariance==EMPTY_VALUE || previousVariance<0.0)
              {
               double initial=ResidualRmsAt(shift+1,length);
               previousVariance=initial*initial;
              }
            double previousScale=MathSqrt(MathMax(previousVariance,0.0));
            double cap=MathMax(HuberConstant,0.10)*MathMax(previousScale,EPSILON_VALUE);
            double clipped=ClampValue(residual,-cap,cap);
            g_volState1[shift]=lambda*previousVariance+(1.0-lambda)*clipped*clipped;
           }
         result=MathSqrt(MathMax(0.0,g_volState1[shift]));
         break;
        }

      case VOL_MAD_RESIDUAL:
         result=ResidualMadAt(shift,length);
         break;

      case VOL_HUBER_SCALE_RESIDUAL:
         result=ResidualHuberScaleAt(shift,length);
         break;

      case VOL_TUKEY_SCALE_RESIDUAL:
         result=ResidualTukeyScaleAt(shift,length);
         break;

      case VOL_QN_RESIDUAL:
         result=ResidualQnAt(shift,length);
         break;

      case VOL_SN_RESIDUAL:
         result=ResidualSnAt(shift,length);
         break;

      case VOL_ABS_RESIDUAL_QUANTILE:
         result=AbsoluteResidualQuantileAt(shift,length);
         break;

      case VOL_CLOSE_TO_CLOSE:
         result=CloseToCloseVolatilityAt(shift,length);
         break;

      case VOL_EWMA_LOG_RETURN:
        {
         double lambda=ClampValue(EwmaLambda,0.001,0.9999);
         double returnValue=LogReturnAt(shift);
         if(seed) g_volState1[shift]=RollingLogReturnVarianceAt(shift,length);
         else
           {
            double previousVariance=g_volState1[shift+1];
            if(!IsFiniteValue(previousVariance) || previousVariance==EMPTY_VALUE || previousVariance<0.0)
               previousVariance=RollingLogReturnVarianceAt(shift+1,length);
            g_volState1[shift]=lambda*previousVariance+(1.0-lambda)*returnValue*returnValue;
           }
         result=MathSqrt(MathMax(0.0,g_volState1[shift]));
         break;
        }

      case VOL_GARCH_11:
        {
         double alpha=ClampValue(GarchAlpha,0.0,0.999);
         double beta=ClampValue(GarchBeta,0.0,0.999);
         if(alpha+beta>=0.999) beta=MathMax(0.0,0.999-alpha);
         double longRunVariance=RollingLogReturnVarianceAt(shift,length);
         if(seed) g_volState1[shift]=longRunVariance;
         else
           {
            double previousVariance=g_volState1[shift+1];
            if(!IsFiniteValue(previousVariance) || previousVariance==EMPTY_VALUE || previousVariance<0.0)
               previousVariance=RollingLogReturnVarianceAt(shift+1,length);
            double previousReturn=LogReturnAt(shift+1);
            double omega=(1.0-alpha-beta)*longRunVariance;
            g_volState1[shift]=omega+alpha*previousReturn*previousReturn+beta*previousVariance;
           }
         result=MathSqrt(MathMax(0.0,g_volState1[shift]));
         break;
        }

      case VOL_PARKINSON:
         result=ParkinsonVolatilityAt(shift,length);
         break;

      case VOL_GARMAN_KLASS:
         result=GarmanKlassVolatilityAt(shift,length);
         break;

      case VOL_ROGERS_SATCHELL:
         result=RogersSatchellVolatilityAt(shift,length);
         break;

      case VOL_YANG_ZHANG:
         result=YangZhangVolatilityAt(shift,length);
         break;

      case VOL_IQR_RESIDUAL:
         result=ResidualIqrAt(shift,length);
         break;
     }

   if(!IsFiniteValue(result) || result==EMPTY_VALUE || result<0.0) result=0.0;
   return result;
  }


//==================== SIZING (history/warmup) ====================
int VolatilityLookbackRequirement()
  {
   int length=SafePeriod(VolatilityPeriod,2);
   int requirement=length+5;
   if(VolatilityMethod==VOL_FAST_SLOW_MAX_LOG_TR)
      requirement=MathMax(length,SafePeriod(FastVolatilityPeriod,2))+5;
   if(VolatilityMethod==VOL_QN_RESIDUAL || VolatilityMethod==VOL_SN_RESIDUAL)
      requirement=MathMin(length,SafePeriod(PairwiseScaleMaximumWindow,10))+5;
   if(VolatilityMethod==VOL_CLOSE_TO_CLOSE || VolatilityMethod==VOL_EWMA_LOG_RETURN ||
      VolatilityMethod==VOL_GARCH_11 || VolatilityMethod==VOL_YANG_ZHANG)
      requirement=length+6;
   return requirement;
  }

bool VolatilityNeedsCenterWindow()
  {
   if(VolatilityMethod==VOL_NATR_WILDER) return true;
   if(VolatilityMethod>=VOL_RESIDUAL_RMS && VolatilityMethod<=VOL_ABS_RESIDUAL_QUANTILE) return true;
   if(VolatilityMethod==VOL_IQR_RESIDUAL) return true;
   return false;
  }

bool VolatilityMethodIsRecursive()
  {
   return (VolatilityMethod==VOL_WILDER_ATR || VolatilityMethod==VOL_EMA_TR ||
           VolatilityMethod==VOL_NATR_WILDER || VolatilityMethod==VOL_LOG_TR_WILDER ||
           VolatilityMethod==VOL_FAST_SLOW_MAX_LOG_TR ||
           VolatilityMethod==VOL_EWMA_RESIDUAL_STD ||
           VolatilityMethod==VOL_HUBER_EWMA_RESIDUAL_STD ||
           VolatilityMethod==VOL_EWMA_LOG_RETURN || VolatilityMethod==VOL_GARCH_11);
  }

int VolatilityWarmupRequirement()
  {
   if(!VolatilityMethodIsRecursive()) return 0;
   return MathMax(50,RecursiveWarmupBars);
  }

#endif // __ROBUST_VOLATILITY_CORE_MQH__