//+------------------------------------------------------------------+
//| 
//| 
//+------------------------------------------------------------------+
#property copyright "PowerLawMBK_MACD Matt Kennel. 2010"
#property link      ""

#property indicator_separate_window
#property indicator_buffers 2
#property indicator_color1 Yellow
#property indicator_color2 Blue

//---- input parameters


extern int CountBars=4000;
extern double PowerExponent1 = 1.0; 
extern double PowerExponent2 = 1.05;
extern int KernelLength = 250; 
extern int offset = 1; 
extern double EMAperiod = 1.5; 

int counter=0;
//---- buffers

double diffbuffer[], emabuffer[]; 
double kernel1[],kernel2[]; 

//+------------------------------------------------------------------+
//| Custom indicator initialization function                         |
//+------------------------------------------------------------------+
int init()
  {
   string short_name;
//---- indicator line
   IndicatorBuffers(2);
   SetIndexStyle(0,DRAW_LINE);
   SetIndexBuffer(0,diffbuffer);

   SetIndexStyle(1,DRAW_LINE);
   SetIndexBuffer(1,emabuffer);
  
   initialize_kernel(kernel1, PowerExponent1);
   initialize_kernel(kernel2, PowerExponent2);
   
   return(0);
  }
  
void initialize_kernel(double& kernel[], double PowerExponent) {
   double kernelsum; 
   ArrayResize(kernel,KernelLength); 
   
   kernelsum = 0.0; 
   for (int i=0; i<KernelLength; i++) {
      kernel[i] = MathPow( (i+1.0)*0.01, -PowerExponent); 
      kernelsum += kernel[i];
   }
   for (i = 0; i<KernelLength; i++) {
      kernel[i] = kernel[i] / kernelsum; 
   }
}  

double convolve(double array[], double kernel[], int n) {
// return sum(i=0..n-1) array[i]*kernel[i]
// conventionally kernel[*] sums to 1, but this is not enforced here.
   double sum = 0.0;
   for (int i=0; i<n; i++) 
      sum += array[i]*kernel[i];

   return(sum); 
}
  
//+------------------------------------------------------------------+
//| SilverTrend_PowerLaw                                             |
//+------------------------------------------------------------------+
int start()
  {   
   CountBars = MathMin(Bars-KernelLength-1,CountBars); 

   int i,shift,counted_bars=IndicatorCounted();
   int i1,i2;
   int limit; 
//----
   double buffer[]; 
   double pema = 2.0/(EMAperiod + 1.0); 
//----
   if (counted_bars<0) return(-1);
   if (counted_bars>0) counted_bars--;
   
   ArrayResize(buffer,KernelLength);

  
   limit = CountBars-counted_bars;
   if (limit < 1) limit = 1; 
   
   for (shift = limit; shift>=0; shift--) 
   { 
      double plma1, plma2; 
              
      // for offset = 1 start computing using previous bar's data only, this may be 
      // needed for causality if using a trading system
      for (i=0; i<KernelLength; i++) {
         buffer[i] = Close[shift+offset+i];
      }

      plma1 = convolve(buffer, kernel1, KernelLength);
      plma2 = convolve(buffer, kernel2, KernelLength); 
      diffbuffer[shift] = (plma2-plma1)/((PowerExponent2-PowerExponent1)*Point);
 
      if (shift < CountBars) {
         emabuffer[shift] = pema *diffbuffer[shift] + (1-pema)*emabuffer[shift+1];
      } else {
         emabuffer[CountBars] = diffbuffer[CountBars]; 
      } 
   } // end for
   return(0);
}
//+------------------------------------------------------------------+