//+------------------------------------------------------------------+
//|                                   AutocorrelationPeriodogram.mq4 |
//|                                         Copyright © 2014, Kirill |
//|                                          StockProgrammer@mail.ru |
//+------------------------------------------------------------------+
#property copyright "Copyright © 2014, Kirill"
#property link      "StockProgrammer@mail.ru"

//+------------------------------------------------------------------+
//| Description                                                      
//|
//| AutocorrelationPeriodogram
//| 
//| v3.0+ - my version: added else {ExtMapBuffer[pos] = ExtMapBuffer[pos+1];}
//|         this same step was taken by the author of CoronaCyclePeriod for MT4
//| v4.0  - added heatmap
//| v4.1  - coded heatmap as in Ehler's book. It works!
//| v4.5  - v4.4 + delete objects
//| v4.8  - corrected recalculation errors see below
//| 
//+------------------------------------------------------------------+


#property indicator_separate_window
#property indicator_buffers 3
#property indicator_color1 Yellow
#property indicator_minimum 0
//---- indicator parameters
extern int     LPPeriod          =        10;             
extern int     HPPeriod          =        48;             
extern int     AvgLength         =        3;
extern color   LineColor         =        Yellow;
extern bool    HeatMap           =        true;
extern int     HeatMapBars       =        1000;



//---- important: these two are used on the next bar!
double R[]; 
double MaxPwr = 0;
datetime TimeOfLastCalculatedBar = 0;
//---- indicator buffers
double ExtMapBuffer[];
double HP[];           //used for calculations only
double RF[];           //used for calculations only
double Corr[]; 
double Pwr[]; 
//----
double  MinPwrLimit = 0.5;
string hex;
int win=-1;   
string short_name;  
int ExtCountedBars=0;
int BarsUsed = 3 + HPPeriod + AvgLength;                              
//---- constants
double alpha,a,b,c1,c2,c3;
double pi = 3.1415926535;
//+------------------------------------------------------------------+
//| Custom indicator initialization function                         |
//+------------------------------------------------------------------+
int init(){
//---- drawing settings
   SetIndexStyle(0,DRAW_LINE,STYLE_SOLID,2,LineColor);  
   
   SetIndexStyle(1,DRAW_NONE);
   SetIndexStyle(2,DRAW_NONE);
         
   IndicatorDigits(MarketInfo(Symbol(),MODE_DIGITS));
   
   SetIndexDrawBegin(0,BarsUsed);                                      
//---- indicator short name
   short_name = "Autocorrelation Periodogram("+LPPeriod+","+HPPeriod+","+AvgLength+","+HeatMapBars+")";   
   IndicatorShortName(short_name);
//---- indicator buffers mapping
   SetIndexBuffer(0,ExtMapBuffer);
   SetIndexBuffer(1,HP);    
   SetIndexBuffer(2,RF);  
//---- prep arrays
ArrayResize(Corr, HPPeriod+1); ArrayInitialize(Corr, 0.0);
ArrayResize(R, HPPeriod+1); ArrayInitialize(R, 0.0);
ArrayResize(Pwr, HPPeriod+1); 
//---- prep constants
   alpha = (MathCos(0.707*2*pi/HPPeriod) + MathSin(0.707*2*pi/HPPeriod) - 1) / MathCos(0.707*2*pi/HPPeriod);
   a = MathExp(-1.414 * pi / LPPeriod);
   b = 2 * a * MathCos(1.414 * pi / LPPeriod);
   c2 = b;
   c3 = -1 * a * a;
   c1 = 1 - c2 - c3;
//---- initialization done
   return(0);
}
int deinit(){
   ObjectsDeleteAll(win);
   return(0);
}
//+------------------------------------------------------------------+
//|                                                                  |
//+------------------------------------------------------------------+
int start(){
   win = WindowFind(short_name); 

   if(Bars<=BarsUsed) return(0);
   ExtCountedBars=IndicatorCounted();
//---- check for possible errors
   if (ExtCountedBars<0) return(-1);
//---- last counted bar will be recounted
   if (ExtCountedBars>0) ExtCountedBars--;
//----
   AutocorrelationIndicator();
//---- done
   return(0);
}
//+------------------------------------------------------------------+
//| AutocorrelationIndicator Calculation                                        |
//+------------------------------------------------------------------+
void AutocorrelationIndicator(){
   int    i,pos = Bars-ExtCountedBars-1;
   int    per;
   
   Print(MaxPwr, " \n",
   R[9], " \n",   
   R[10], " \n",
   R[11], " \n",
   R[12], " \n",
   R[13], " \n",
   R[14], " \n");

//---- zero initial bars
   if(ExtCountedBars<1){
      for(i=1;i<BarsUsed;i++){
         HP[Bars-i]=0;
         RF[Bars-i]=0;
         ExtMapBuffer[Bars-i]=0; 
         pos--;
      }
   }
   else{ //this part is specific to this indicator. Took me exactly 2.5 hours on 15/4/2014 to figure out how to do this //v4.8
      pos = iBarShift(NULL,0,TimeOfLastCalculatedBar,true); //most of the time pos will be = 2, then will be reduced to 1 (see below)      
      if(pos == -1){ //have to recalculate the whole indicator since last stopping point cannot be found
                     //this is because the R[] array and MaxPwr are associated with bars
                     //if you cannot find the last bar at which the indicator was calculated, then
                     //R[] and MaxPwr will be incorrectly transferred to the wrong bar
         Print("Recalculating entire indicator... may take a while");
         
         pos = Bars - 1;
         for(i=1;i<BarsUsed;i++){
            HP[Bars-i]=0;
            RF[Bars-i]=0;
            ExtMapBuffer[Bars-i]=0; 
            pos--;
         }
      }
      else{
         pos--; //reduce by 1 bar, because bar [pos] has already been calculate
      }
   }

//---- calculate all other bars       
   double   X = 0;
   double   Y = 0;      
   double 	Sx = 0;
   double 	Sy = 0;
   double 	Sxx = 0;
   double 	Syy = 0;
   double 	Sxy = 0;     
   
   double CosPart = 0;
   double SinPart = 0;
   double SqSum = 0;     
   
   
   while(pos>0)   //not >=0 , but strictly >0
                  //this indicator can only be calculated for completed bars! //v4.8
                  //otherwise the transfer array R[] and transfer variable MaxPwr
                  //will get incorrect values!
   {
   //Highpass filter low-frequency cyclic components
   //Using two-pole Highpass filter to achieve 12db/octave attenuation roll-on (with falling frequency)   
      HP[pos] = MathPow((1-alpha/2),2) * (Close[pos]-2*Close[pos+1]+Close[pos+2]) + 2*(1-alpha)*HP[pos+1] - MathPow((1-alpha),2)*HP[pos+2];

   //Lowpass filter high-frequency cyclic components
      RF[pos] = 
                           c1 * (HP[pos] + HP[pos+1]) / 2 
                              + c2 * RF[pos+1] + c3 * RF[pos+2];
      
   //Autocorrelation indicator calcs start here
   //Pearson correlation for each value of lag
      for(int Lag = 3; Lag <= HPPeriod; Lag++){
         int M = AvgLength;
         if (AvgLength == 0){
            M=Lag;
         }
         X = 0;
         Y = 0;      
         Sx = 0;
      	Sy = 0;
      	Sxx = 0;
      	Syy = 0;
      	Sxy = 0;         
         for(int avcount = 0; avcount <= M-1; avcount++){
            X = RF[pos + avcount];
            Y = RF[pos + Lag + avcount];
            Sx = Sx + X;
            Sy = Sy + Y;
            Sxx = Sxx + X*X;
            Sxy = Sxy + X*Y;
            Syy = Syy + Y*Y;         
         }
         
         if((M*Sxx - Sx*Sx)*(M*Syy - Sy*Sy) > 0){
            Corr[Lag] = (M*Sxy - Sx*Sy)/MathSqrt((M*Sxx - Sx*Sx)*(M*Syy - Sy*Sy));
         }
      }

      for(per = LPPeriod; per <= HPPeriod; per++){
         CosPart = 0;
         SinPart = 0;
         SqSum = 0;
         
         for(int N = 3; N <= HPPeriod; N++){
            CosPart = CosPart + Corr[N]*MathCos(2*pi*N / per); 
            SinPart = SinPart + Corr[N]*MathSin(2*pi*N / per);
         }
         SqSum = MathPow(CosPart,2) + MathPow(SinPart,2);



         R[per] = 0.2*MathPow(SqSum,2) + 0.8*R[per];            
      }   

   //Find Maximum Power Level for Normalization
      MaxPwr =  0.995*MaxPwr;
      for(per = LPPeriod; per <= HPPeriod; per++){
         if(R[per] > MaxPwr){
            MaxPwr = R[per];
         }
      }
      for(per = LPPeriod; per <= HPPeriod; per++){
         if(MaxPwr != 0){
            Pwr[per] = R[per] / MaxPwr;      
         } else{  //so that prev value of Pwr[per] gets erased
            Pwr[per] = 0;
         }
      }

   //Compute the dominant cycle using the CG of the spectrum
      double Spx = 0;
      double Sp = 0;
      for(per = LPPeriod; per <= HPPeriod; per++){
         if(Pwr[per] >= MinPwrLimit){
            Spx = Spx + per * Pwr[per];
            Sp = Sp + Pwr[per];
         }
      }
      if(Sp != 0){
         ExtMapBuffer[pos] = Spx / Sp;
      } else{  //v3.0
         ExtMapBuffer[pos] = ExtMapBuffer[pos+1];
      }
       
      
   //Plot heatmap if requested
      if(HeatMap==true) {
         if((HeatMapBars > 0 && pos <= HeatMapBars))
            for(int n = LPPeriod; n<= HPPeriod; n++)
            {
               int Color1=0;
               int Color2=0;
      		   int Color3=0;
            
               if(Pwr[n] > 0.5)
               {
                  Color1 = 255;
                  Color2 = 255*(2*Pwr[n]-1);
               }
               else{
                  Color1 = 2*255*Pwr[n];
               }
  
               PlotRect((n+1),win+"AP"+(n+1)+"_"+pos+"_"+iTime(NULL,0,pos), RGB(Color1, Color2, Color3),1,pos);
            }
      }
   
      TimeOfLastCalculatedBar = iTime(NULL,0,pos);
         
      pos--;
   }
}

void PlotRect(double value,string name,color clr,double width,int bar)
{   
   ObjectDelete(name);
   ObjectCreate(name,OBJ_RECTANGLE,win,iTime(NULL,0,bar+1),value-0.5*width,iTime(NULL,0,bar),value+0.5*width);
   ObjectSet(name,OBJPROP_COLOR,clr);
}

int RGB(int R1, int G1, int B1)
{
   return (256*(256*B1 + G1) + R1);
}