//------------------------------------------------------------------
#property copyright "© mladen, 2016, MetaQuotes Software Corp."
#property link      "www.forex-tsd.com, www.mql5.com"
//------------------------------------------------------------------
#property indicator_chart_window
#property indicator_buffers 7
#property indicator_plots   7
#property indicator_label1  "ema"
#property indicator_type1   DRAW_LINE
#property indicator_color1  clrDarkGray
#property indicator_style1  STYLE_DOT
#property indicator_label2  "ema - fitted past"
#property indicator_type2   DRAW_LINE
#property indicator_color2  clrDarkGray
#property indicator_width2  2
#property indicator_label3  "ema - fitted future"
#property indicator_type3   DRAW_LINE
#property indicator_color3  clrRed
#property indicator_width3  2
#property indicator_label4  "ema - fitted past band up"
#property indicator_type4   DRAW_LINE
#property indicator_color4  clrLimeGreen
#property indicator_label5  "ema - fitted past band down"
#property indicator_type5   DRAW_LINE
#property indicator_color5  clrRed
#property indicator_label6  "ema - fitted future band up"
#property indicator_type6   DRAW_LINE
#property indicator_color6  clrLimeGreen
#property indicator_label7  "ema - fitted future band down"
#property indicator_type7   DRAW_LINE
#property indicator_color7  clrRed

//
//
//
//
//

input int      period   = 5;        // EMA period
input double   Factor   = 2;        // Bands deviations
input int      LastBar  = 0;        // Last bar in the past data
input int      PastBars = 200;      // Number of past bars
input int      FutBars  = 50;       // Number of bars to predict 
input int      HarmNo   = 8;        // Total number of harmonics in model
input double   FreqTOL  = 0.0001;   // Tolerance of frequency calculations
input bool     EveryTick = true;    // Calculate on every tick?

double in[],pv[],fv[],upp[],dnp[],upf[],dnf[];

//------------------------------------------------------------------
//
//------------------------------------------------------------------
//
//
//
//
//

int OnInit()
{
   SetIndexBuffer(0,in ,INDICATOR_DATA); 
   SetIndexBuffer(1,pv ,INDICATOR_DATA); 
   SetIndexBuffer(2,fv ,INDICATOR_DATA); 
   SetIndexBuffer(3,upp,INDICATOR_DATA); 
   SetIndexBuffer(4,dnp,INDICATOR_DATA); 
   SetIndexBuffer(5,upf,INDICATOR_DATA); 
   SetIndexBuffer(6,dnf,INDICATOR_DATA); 
   return(0);
}

//------------------------------------------------------------------
//
//------------------------------------------------------------------
//
//
//
//
//

double x[],xp[],xf[];
int OnCalculate (const int rates_total,
                 const int prev_calculated,
                 const int begin,
                 const double& price[] )
{
   if (Bars(_Symbol,_Period)<rates_total) return(0);
      
      if (rates_total<PastBars) return(0);
      datetime time[1]; if (CopyTime(_Symbol,_Period,0,1,time)!=1) return(0);
      static datetime lastTime=0; 
                  if (lastTime!=0 && (lastTime==time[0] && (!EveryTick || LastBar>0))) return(rates_total);
                                      lastTime =time[0];

      //
      //
      //
      //
      //
                  
         int lastBar  = MathMax(LastBar,0);
         int pastBars = MathMin(PastBars,rates_total);
         int futBars  = MathMin(pastBars-1,FutBars);
         if (ArraySize(x) != pastBars) { ArrayResize(x ,pastBars); ArrayResize(xp,pastBars); }
         if (ArraySize(xf)!= futBars+1)  ArrayResize(xf,futBars+1); 

            double av=0; int startAt=MathMax(rates_total-pastBars-lastBar+1,rates_total-1);
               for (int i=(int)MathMax(prev_calculated-1,0); i<rates_total; i++) in[i] = iEma(price[i],period,i,rates_total); 
               for (int i=startAt                          ; i<rates_total; i++) av+=in[i]; av/=pastBars;
                     ArrayInitialize(xp,av);
                     ArrayInitialize(xf,av);
               for (int i=lastBar; i<pastBars+lastBar; i++) x[i-lastBar] = in[rates_total-i-1];
   
   //
   //
   //
   //
   //
   
   double w,m,c,s,dev=0;
   for(int harm=1; harm<=HarmNo && !IsStopped();harm++)
   {
      iQuiFerFreq(pastBars,FreqTOL,x,xp,w,m,c,s);
      for(int i=0; i<pastBars; i++) 
      {
                           xp[i] += m+c*cos(w*i)+s*sin(w*i);
         if (i<=futBars)   xf[i] += m+c*cos(w*i)-s*sin(w*i);
         if (harm==HarmNo) dev   += MathPow(x[i]-xp[i],2);
      }       
   }
   dev = MathSqrt(dev/(pastBars+1))*Factor;
   for (int i=0; i<pastBars; i++) 
   {
      pv[rates_total-i-1]  = xp[i]; ArrayCopy(fv,xf,rates_total-futBars); 
      upp[rates_total-i-1] = pv[rates_total-i-1] + dev;
      dnp[rates_total-i-1] = pv[rates_total-i-1] - dev;
      upf[rates_total-i-1] = (i<=futBars) ? fv[rates_total-i-1] + dev : EMPTY_VALUE;
      dnf[rates_total-i-1] = (i<+futBars) ? fv[rates_total-i-1] - dev : EMPTY_VALUE;
   }      
   
   //
   //
   //
   //
   //
   
   PlotIndexSetInteger(1,PLOT_SHIFT,       -lastBar);   PlotIndexSetInteger(1,PLOT_DRAW_BEGIN,rates_total-pastBars);
   PlotIndexSetInteger(3,PLOT_SHIFT,       -lastBar);   PlotIndexSetInteger(3,PLOT_DRAW_BEGIN,rates_total-pastBars);
   PlotIndexSetInteger(4,PLOT_SHIFT,       -lastBar);   PlotIndexSetInteger(4,PLOT_DRAW_BEGIN,rates_total-pastBars);
   PlotIndexSetInteger(2,PLOT_SHIFT,futBars-lastBar-1); PlotIndexSetInteger(2,PLOT_DRAW_BEGIN,rates_total-futBars);
   PlotIndexSetInteger(5,PLOT_SHIFT,futBars-lastBar-1); PlotIndexSetInteger(5,PLOT_DRAW_BEGIN,rates_total-futBars);
   PlotIndexSetInteger(6,PLOT_SHIFT,futBars-lastBar-1); PlotIndexSetInteger(6,PLOT_DRAW_BEGIN,rates_total-futBars);
   return(rates_total);
}

//------------------------------------------------------------------
//
//------------------------------------------------------------------
//
//
//
//
//

double workEma[][1];
double iEma(double price, double tperiod, int r, int bars, int instanceNo=0)
{
   if (ArrayRange(workEma,0)!= bars) ArrayResize(workEma,bars); 

   //
   //
   //
   //
   //

   if (price== EMPTY_VALUE) price=0;
   workEma[r][instanceNo] = price;
   if (r>0 && tperiod>1)
          workEma[r][instanceNo] = workEma[r-1][instanceNo]+(2.0/(1.0+tperiod))*(price-workEma[r-1][instanceNo]);
   return(workEma[r][instanceNo]);
}

//---------------------------------------------------------------------------------------------
//
//    Quinn and Fernandes algorithm - originally for mt4 coded by qpwr, this version by mladen
//
//---------------------------------------------------------------------------------------------
//
//
//
//
//
//

double _z[];
void iQuiFerFreq(int np, double freqTol, double& _x[], double& values[], double& w, double& m, double& c, double& s)
{
   if (ArraySize(_z)!=np) ArrayResize(_z,np);
   
   //
   //
   //
   //
   //
   
   double a=0.0, b=2.0; _z[0]=_x[0]-values[0];
   while(MathAbs(a-b)>freqTol)
   {
      a=b;
            _z[1]=_x[1]-values[1]+a*_z[0];
               double num=_z[0]*_z[1];
               double den=_z[0]*_z[0];
               for(int i=2;i<np;i++)
               {
                  _z[i] = _x[i]-values[i]+a*_z[i-1]-_z[i-2];
                     num  += _z[i-1]*(_z[i]+_z[i-2]);
                     den  += _z[i-1]*_z[i-1];
               }
               b = (den!=0) ? num/den : 0;
   }
   w=MathArccos(MathMin(MathMax(b/2.0,-1),1));

   //
   //
   //
   //
   //
      
   double Sc=0,Ss=0,Scc=0,Sss=0,Scs=0,Sx=0,Sxc=0,Sxs=0;;
      for(int i=0;i<np;i++)
      {
         double wcos = cos(w*i);
         double wsin = sin(w*i);
                Sc  += wcos;
                Ss  += wsin;
                Scc += wcos*wcos;
                Sss += wsin*wsin;
                Scs += wcos*wsin;
                Sx  += (_x[i]-values[i]);
                Sxc += (_x[i]-values[i])*wcos;
                Sxs += (_x[i]-values[i])*wsin;
      }
      Sc/=np; Ss/=np; Scc/=np; Sss/=np; Scs/=np; Sx/=np; Sxc/=np; Sxs/=np; m=Sx; c=0.0;  s=0.0;
      if(w!=0.0)
      {
         double den = MathPow(Scs-Sc*Ss,2)-(Scc-Sc*Sc)*(Sss-Ss*Ss);
                c = ((Sxs-Sx*Ss)*(Scs-Sc*Ss)-(Sxc-Sx*Sc)*(Sss-Ss*Ss))/den;
                s = ((Sxc-Sx*Sc)*(Scs-Sc*Ss)-(Sxs-Sx*Ss)*(Scc-Sc*Sc))/den;
                m = Sx-c*Sc-s*Ss;
      }
   return;
}