//+------------------------------------------------------------------+
//|                                                   AlgLibTest.mq5 |
//|                                  Copyright 2022, MetaQuotes Ltd. |
//|                                             https://www.mql5.com |
//+------------------------------------------------------------------+
#property copyright "Arthur Albano"
#property link      "https://www.mql5.com"
#property version   "1.00"


//+------------------------------------------------------------------+
//| Script program start function                                    |
//+------------------------------------------------------------------+

#include <Math\Alglib\alglib.mqh>

//+------------------------------------------------------------------+
//| Derived class from CNDimensional_PFunc                           |
//+------------------------------------------------------------------+
class CNDimensional_CX_2_Func : public CNDimensional_PFunc
  {
public:
                     CNDimensional_CX_2_Func(void);
                    ~CNDimensional_CX_2_Func(void);
   virtual void      PFunc(double &c[], double &K[], double &func, CObject &obj);
  };
//+------------------------------------------------------------------+
//| Constructor                                                      |
//+------------------------------------------------------------------+
CNDimensional_CX_2_Func::CNDimensional_CX_2_Func(void) {}
//+------------------------------------------------------------------+
//| Destructor                                                       |
//+------------------------------------------------------------------+
CNDimensional_CX_2_Func::~CNDimensional_CX_2_Func(void) {}
//+------------------------------------------------------------------+
//| This callback calculates option premiums not distinguishing american/european  |
//+------------------------------------------------------------------+
void CNDimensional_CX_2_Func::PFunc(double &m[], double &K[], double &func, CObject &obj)
  {
   double d1 = (m[0]*K[0] - m[1]*_S + m[3]);
   double c = 4.0*(_T)*pow(m[2]*_S,2.0);
   double d2 = (m[0]*K[0] - m[1]*_S + m[3]);
   func = (1.0/2.0)*sqrt(pow(d1,2.0)+c) - _Z * (d2) / 2.0;
  }
//+------------------------------------------------------------------+
//| Nonlinear fitting using function value only                      |
//+------------------------------------------------------------------+
double LSFit_OPTIONS()
  {
//--- create variables
   CMatrixDouble           K;    //strikes (K)
   double                  y[];  //premiums (X0)
   double                  c[];  //overfitting coeffs
   double                  w[];  //weights (deals, for day close)
   CObject                 obj;
   CNDimensional_CX_2_Func fcx2func;
   CNDimensional_Rep       frep;
//--- load options symbols
   MqlTick _MqlTicks[];
   int lTicksReceived=0;
   int iPutCount=0;
   string strISIN,strSymbol;
   string strS=StringSubstr(SymbolInfoString(_SYMBOL,SYMBOL_ISIN),2,5);
   if(StringSubstr(strS,4,1)=="9")
      strS=StringSubstr(SymbolInfoString(_SYMBOL,SYMBOL_ISIN),2,4)+"11";
   SymbolSelect(strS,true);
   double dbPrice=0.0;
   for(int i=0; i<SymbolsTotal(false); i++)
     {
      strSymbol=SymbolName(i,false);
      strISIN=SymbolInfoString(strSymbol,SYMBOL_ISIN);
      if(
         (StringSubstr(strISIN,0,8)==StringSubstr(SymbolInfoString(_SYMBOL,SYMBOL_ISIN),0,8))&&
         (SymbolInfoDouble(strSymbol,SYMBOL_OPTION_STRIKE)>0.0)&&
         (SymbolInfoInteger(strSymbol,SYMBOL_EXPIRATION_TIME)==SymbolInfoInteger(_SYMBOL,SYMBOL_EXPIRATION_TIME))&&
         (SymbolInfoInteger(strSymbol,SYMBOL_OPTION_MODE)==SymbolInfoInteger(_SYMBOL,SYMBOL_OPTION_MODE))
      )
        {
         SymbolSelect(strSymbol,true);
         lTicksReceived=CopyTicks(strSymbol,_MqlTicks,COPY_TICKS_ALL,Today());
         dbPrice=SymbolPrice(strSymbol,inpCalcMode);
         if(
            (dbPrice>0) &&
            (Day(SymbolInfoInteger(strSymbol,SYMBOL_TIME))==Day(SymbolInfoInteger(strS,SYMBOL_TIME))) &&
            ((InpFilterBE==true)?(_Z*(SymbolInfoDouble(strS,SYMBOL_LAST)-SymbolInfoDouble(strSymbol,SYMBOL_OPTION_STRIKE))<dbPrice):true)
          )
           {
            iPutCount++;
            K.Resize(iPutCount,1);                                                                 //-- allocation x
            K[iPutCount-1].Set(0,SymbolInfoDouble(strSymbol,SYMBOL_OPTION_STRIKE));                //-- strikes
            ArrayResize(y,iPutCount);                                                              //-- allocation y
            //dbPrice = MathMax(_Z*(SymbolInfoDouble(strS,SYMBOL_LAST)-SymbolInfoDouble(strSymbol,SYMBOL_OPTION_STRIKE)),dbPrice);
            y[iPutCount-1]=dbPrice;                                                                //-- market premiums
            ArrayResize(w,iPutCount);                                                              //-- allocation w (weights)
            w[iPutCount-1]=(SymbolInfoDouble(strSymbol,SYMBOL_ASK)+SymbolInfoDouble(strSymbol,SYMBOL_BID))/fabs(SymbolInfoDouble(strSymbol,SYMBOL_ASK)-SymbolInfoDouble(strSymbol,SYMBOL_BID));         //-- session deals is the right thing to do?
            if(InpVerbose==true)
               Print(strSymbol,";",
                     K[iPutCount-1][0],";",
                     DoubleToString(y[iPutCount-1],(int)SymbolInfoInteger(strSymbol,SYMBOL_DIGITS)),";",
                     DoubleToString(SymbolInfoDouble(strSymbol,SYMBOL_BID),(int)SymbolInfoInteger(strSymbol,SYMBOL_DIGITS)),";",
                     DoubleToString(SymbolInfoDouble(strSymbol,SYMBOL_ASK),(int)SymbolInfoInteger(strSymbol,SYMBOL_DIGITS)),";",
                     DoubleToString(SymbolInfoDouble(strSymbol,SYMBOL_LAST),(int)SymbolInfoInteger(strSymbol,SYMBOL_DIGITS)),";",
                     (datetime)SymbolInfoInteger(strSymbol,SYMBOL_TIME));
           }
        }
     }
//--- allocation
   ArrayResize(c,5);
//--- initialization
   c[0] = 1.0; //
   c[1] = 1.0; //
   c[2] = 0.2; //
   c[3] = 0.0; //
   c[4] = 1.0; //
   _S   = SymbolPrice(strS,inpCalcMode);
   _T   = MathMax(0.0,(SymbolInfoInteger(_SYMBOL,SYMBOL_EXPIRATION_TIME)-SymbolInfoInteger(strS,SYMBOL_TIME))/(365.0*86400.0));
   _Z   = ((SymbolInfoInteger(_SYMBOL,SYMBOL_OPTION_RIGHT)==SYMBOL_OPTION_RIGHT_CALL)?1.0:-1.0);
   _K   = SymbolInfoDouble(_SYMBOL,SYMBOL_OPTION_STRIKE);
   _X   = 0.0;
//--- create variables
   int maxits=0;
   int info;
   CLSFitStateShell state;
   double diffstep=0.0001;
//--- Fitting withweights
   CAlglib::LSFitCreateF(K,y,c,diffstep,state);
//--- function call
   double epsf = 0.000001;
   double epsx = 0.000001;
   CAlglib::LSFitSetCond(state,epsf,epsx,maxits);
//--- function call
   CAlglib::LSFitFit(state,fcx2func,frep,0,obj); //puts
//--- function call
   CLSFitReportShell rep;
   CAlglib::LSFitResults(state,info,c,rep);
//--- Clear weights
   ArrayFill(w,0,ArraySize(w),1.0);
//--- function call
   CAlglib::LSFitCreateWF(K,y,w,c,diffstep,state);
// allocation of lower boundaries
   double bndl[5];
// all lower boundaries must be positive
   bndl[0] = 0.5;
   bndl[1] = 0.5;
   bndl[2] = 0.0;
   bndl[3] = -100;
   bndl[4] = 0.5;
// allocation of upper boundaries
   double bndu[5];
   bndu[0] = 1.5;
   bndu[1] = 1.5;
   bndu[2] = 1.0;
   bndu[3] = +100;
   bndu[4] = 1.5;
//--- function call
   CAlglib::LSFitSetBC(state,bndl,bndu);
//--- function call
   CAlglib::LSFitSetCond(state,epsf,epsx,maxits);
//--- function call
   CAlglib::LSFitFit(state,fcx2func,frep,0,obj);
//--- function call
   CAlglib::LSFitResults(state,info,c,rep);
   _X = (1.0/2.0)*sqrt(pow(c[0]*_K-c[1]*_S+c[3],2.0)+4*_T*pow(c[2]*_S,2.0))-_Z*(c[0]*_K-c[1]*_S+c[3])/2.0; //Premium
   _r = (((1.0/2.0)*sqrt(pow(c[0]*_S-c[1]*_S+c[3],2.0)+4*_T*pow(c[2]*_S,2.0))-_Z*(c[0]*_S-c[1]*_S+c[3])/2.0)/_S)/sqrt(_T); //Return
   if(InpVerbose==true)
     {
      printf("Option; %s", _SYMBOL);
      printf("K; %.2f", SymbolInfoDouble(_SYMBOL,SYMBOL_OPTION_STRIKE));
      printf("S; %.2f", _S);
      printf("T; %.5f", _T);
      printf("Z; %.5f", _Z);
      printf("r; %.5f ", _r);
      printf("X; %.2f",_X);
     };

   if(InpGraph)
     {
      double Graphic[][3];
      ArrayResize(Graphic,K.Size(),0);

      for(int i=0; i<K.Size(); i++)
        {
         Graphic[i][0]=K[i][0];
         Graphic[i][1]=y[i];
         Graphic[i][2]=(1.0/2.0)*sqrt(pow(c[0]*K[i][0]-c[1]*_S+c[3],2.0)+4*_T*pow(c[2]*_S,2.0))-_Z*(c[0]*K[i][0]-c[1]*_S+c[3])/2.0;
        };

      ArraySort(Graphic);

      const int size = ArraySize(Graphic)/3;

      double GraphicK[];
      ArrayResize(GraphicK,size,0);
      double GraphicX0[];
      ArrayResize(GraphicX0,size,0);
      double GraphicX[];
      ArrayResize(GraphicX,size,0);

      for(int i=0; i<size; i++)
        {
         GraphicK[i]  = Graphic[i][0];
         GraphicX0[i] = Graphic[i][1];
         GraphicX[i]  = Graphic[i][2];
        }
      
      if(ObjectFind(0,"OptionGraphic")==0) ObjectDelete(0,"OptionGraphic");
            
      CGraphic _CGraphic;
      _CGraphic.Create(0,"OptionGraphic",0,0,100,1600/4+0,900/4+100);
      _CGraphic.BackgroundMain(strS+" ("+DoubleToString(_S,2)+")"+" @ "+TimeToString(TimeCurrent(),TIME_DATE));
      _CGraphic.BackgroundMainSize(12);
      _CGraphic.XAxis().NameSize(12);
      _CGraphic.XAxis().Name("Strike (K)");
      _CGraphic.XAxis().ValuesFormat("%.2f");
      _CGraphic.YAxis().NameSize(12);
      _CGraphic.YAxis().Name("Premium (X||X0)");
      _CGraphic.YAxis().ValuesFormat("%.2f");

      CCurve *CCUrveX0  = _CGraphic.CurveAdd(GraphicK,GraphicX0,0x009900,CURVE_POINTS,"Premium");
      CCurve *CCUrveX   = _CGraphic.CurveAdd(GraphicK,GraphicX ,0x000099,CURVE_LINES,"OPM ("+DoubleToString(_r*100,1)+" %)");
      CCUrveX.LinesSmooth(true);
      
      double Sx[2]={_S,_S};
      double Sxmax = MathMax(MathMax(GraphicX0[0],GraphicX[0]),MathMax(GraphicX0[ArraySize(GraphicX0)-1],GraphicX[ArraySize(GraphicX)-1]));
      double Sy[2]={0.0,Sxmax};
      CCurve *CCurveS   = _CGraphic.CurveAdd(Sx,Sy,ColorToARGB(clrChartreuse,255),CURVE_STEPS,strS);
      
      double Kx[2]={_K,_K};
      double Kxmax = MathMax(MathMax(GraphicX0[0],GraphicX[0]),MathMax(GraphicX0[ArraySize(GraphicX0)-1],GraphicX[ArraySize(GraphicX)-1]));
      double Ky[2]={0.0,Kxmax};
      CCurve *CCurveKVert = _CGraphic.CurveAdd(Kx,Ky,ColorToARGB(clrRed,255),CURVE_STEPS,_SYMBOL);
      Kx[0]=GraphicK[0];
      Kx[1]=GraphicK[ArraySize(GraphicK)-1];
      Ky[0]=_X;
      Ky[1]=_X;
      CCurve *CCurveKHoriz = _CGraphic.CurveAdd(Kx,Ky,ColorToARGB(clrRed,255),CURVE_STEPS,_SYMBOL);
      
      _CGraphic.CurvePlotAll();
      _CGraphic.Update(true);
      _CGraphic.Redraw(true);

     };

   return(_X);
  }
//+------------------------------------------------------------------+