//+------------------------------------------------------------------+
//|                                              PolyFitBands_v1.mq4 |
//|                                Copyright � 2007, TrendLaboratory |
//|            http://finance.groups.yahoo.com/group/TrendLaboratory |
//|                                   E-mail: igorad2003@yahoo.co.uk |
//+------------------------------------------------------------------+
#property copyright "Copyright � 2007, TrendLaboratory"
#property link      "http://finance.groups.yahoo.com/group/TrendLaboratory"

#property indicator_chart_window
#property indicator_buffers   6
#property indicator_color1    Orange
#property indicator_width1    2
#property indicator_color2    Orange
#property indicator_width2    1
#property indicator_style2    2       
#property indicator_color3    Orange
#property indicator_width3    2
#property indicator_color4    Orange
#property indicator_width4    1
#property indicator_style4    2 
#property indicator_color5    Orange
#property indicator_width5    2
#property indicator_color6    Orange
#property indicator_width6    1
#property indicator_style6    2  

//---- external variables
extern int     Price       =   0;   //Apply to Price(0-Close;1-Open;2-High;3-Low;4-Median price;5-Typical price;6-Weighted Close)
extern int     Length      =  30;   //Number of bars for an evaluation
extern int     Degree      =   3;   //Degree of a Polynomial(no more 12)
extern int     MA_Length   =   1;   //Period of Preliminary Smoothing
extern int     MA_Mode     =   0;   //Mode of MA:0-SMA,1-EMA,2-Wilders(SMMA),3-LWMA
extern int     ProjLength  =   5;   //Length of Projection 
extern double  K_Sigma     =   1;   //Multiplier of Sigma(Standard Deviation)
extern int     FitMode     =   0;   //Fitting Mode: 0-Fitting,1-Moving
extern int     CenterMode  =   0;   //Center Line Plot: 0-without,1-with 
extern int     CountBars   = 500;   //Maximum Number of Bars(0-all Bars) 
//---- indicator buffers
double PolyFitMA[];
double Projection[];
double UpBand[];
double UpProject[];
double DnBand[];
double DnProject[];
double Polynom[];
double Project[];
//---- global variables
double PriceArray[];
int cBars;
//+---------------------------------------------------------------------------+
//| Custom indicator initialization function                                  |
//+---------------------------------------------------------------------------+
int init()
{
   IndicatorDigits(MarketInfo(Symbol(), MODE_DIGITS )+2);
   IndicatorBuffers(8);
   SetIndexStyle(0, DRAW_LINE);
   SetIndexStyle(1, DRAW_LINE);
   SetIndexStyle(2, DRAW_LINE);
   SetIndexStyle(3, DRAW_LINE);
   SetIndexStyle(4, DRAW_LINE);
   SetIndexStyle(5, DRAW_LINE);
   
   SetIndexBuffer(0, PolyFitMA);
   SetIndexBuffer(1, Projection);
   SetIndexBuffer(2, UpBand);
   SetIndexBuffer(3, UpProject);
   SetIndexBuffer(4, DnBand);
   SetIndexBuffer(5, DnProject);
   SetIndexBuffer(6, Polynom);
   SetIndexBuffer(7, Project);
   
   if(CountBars == 0) cBars = Bars;
   else cBars = CountBars;
   
   SetIndexDrawBegin( 0,Bars-cBars+(Length+MA_Length));
   SetIndexDrawBegin( 2,Bars-cBars+2*(Length+MA_Length));
   SetIndexDrawBegin( 4,Bars-cBars+2*(Length+MA_Length));
         
   SetIndexShift(1,ProjLength);
   SetIndexShift(3,ProjLength);
   SetIndexShift(5,ProjLength);
   
   string short_name="PolyFitBands("+Price+","+Length+","+Degree+")";
   IndicatorShortName(short_name);
   SetIndexLabel(0,"PolyFitCenter");
   SetIndexLabel(1,"CenterProject");
   SetIndexLabel(2,"UpPolyBand");
   SetIndexLabel(3,"UpProjection");
   SetIndexLabel(4,"DnPolyBand");
   SetIndexLabel(5,"DnProjection");

   ArrayResize(PriceArray, Length);
   
   return(0);
}

double PolyFit(double price[],int deg,int len,int bar)
{
double   result=0,Sum=0;
double   AX[12,12],BX[12],CX[12],ZX[12],Pow[24];
int      Row[12];
  
// FILL MATRIX FIRST
   if (len<=1)	Sum=price[0];
   else
   {	
   
   for (int j=1;j<=deg+1;j++) BX[j]=0;
		
		for (int k=1;k<=len;k++) 
  	   {
  	   double YK=price[len-k];
  	   double XK=k;
  	   double Prod=1;
  	      for (j=1;j<=deg+1;j++)
    	   {
    	   BX[j]+= YK*Prod;
    	   Prod*=XK;
  		   }
      }
      
   for (j=1;j<=2*deg;j++) Pow[j]=0;
	
	Pow[0]=len;
	
      for (k=1;k<=len;k++)
	   {
  	   XK=k;
  	   Prod=k;
  	      for (j=1;j<=2*deg;j++)
  	      {
    	   Pow[j]+=Prod;
    	   Prod*=XK;
  		   }
      }	
	
	   for (j=1;j<=deg+1;j++)
	   {	
         for (int l=1;l<=deg+1;l++) 
  	      AX[j,l]=Pow[j+l-2];
      }	
//NOW SOLVE FOR COEFFICIENTS

	for (j=1;j<=deg+1;j++) Row[j]=j;
  	
      for (int i=1;i<=deg;i++)
  	   {
  	      for (k=i+1;k<=deg+1;k++)
         {
            if(MathAbs(AX[Row[k],i]) > MathAbs(AX[Row[i],i]))
      	   {	
      	   int Temp=Row[i];
      	   Row[i]=Row[k];
      	   Row[k]=Temp;
            }
         }

         for (k=i+1;k<=deg+1;k++)
    	   {
    	   if(AX[Row[i],i]!=0) AX[Row[k],i]=AX[Row[k],i]/AX[Row[i],i];
    		   for (l=i+1;l<=deg+1;l++)
      	   AX[Row[k],l]=AX[Row[k],l]-AX[Row[k],i]*AX[Row[i],l];
    		}
		}

   ZX[1]=BX[Row[1]];
      
      for (k=2;k<=deg+1;k++)
  	   {
  	   Sum=0;
      for (l=1;l<=k-1;l++) Sum+=AX[Row[k],l]*ZX[l];
  	   
  	   ZX[k]=BX[Row[k]]-Sum;
	   }
     
	CX[deg+1]=ZX[deg+1]/AX[Row[deg+1],deg+1];
      	
      for (k=deg;k>=1;k--)
  	   {
  	   Sum=0;
  	   l=k+1;  		
  	   for (l=k+1;l<=deg+1;l++) Sum += AX[Row[k],l]*CX[l];
  	   CX[k]=(ZX[k]-Sum)/AX[Row[k],k];
	   }
    
//NOW COMPUTE NEXT POINT IN SERIES AND RETURN}

	Sum=CX[deg+1];
	for (k=deg;k>=1;k--) Sum=CX[k]+Sum*(len+bar);
   
   }
   return (Sum);
}	
//+------------------------------------------------------------------+
//| PolyFitBands_v1                                                     |
//+------------------------------------------------------------------+
int start()
{
   int i, j, k, m, shift, counted_bars=IndicatorCounted(), limit;
        
   
   if ( counted_bars < 0 )  return(0);
   if ( counted_bars ==0 )  limit=cBars-1; 
         
   if (FitMode==1 && counted_bars > 0) {limit = Bars - counted_bars; int len = Length+MA_Length;}
   else
   if (FitMode==0) {limit = Length; len = cBars - Length+1;} 
   
   if ( counted_bars < 1 || FitMode==0)
   { 
      for(i=1;i<len;i++)
      { 
      Polynom[cBars-i]=0; 
      PolyFitMA[cBars-i]=EMPTY_VALUE;    
      UpBand[cBars-i]=EMPTY_VALUE;
      DnBand[cBars-i]=EMPTY_VALUE;
      }
   }   
   
   if ( ProjLength > 0 )
   {  
      for(i=1;i<Length;i++)
      {
      Projection[i]=EMPTY_VALUE;  
      UpProject[i]=EMPTY_VALUE;
      DnProject[i]=EMPTY_VALUE;
      }
   } 
     
   for(shift=0;shift<limit;shift++ )
   {
      if((shift==0 && FitMode==0) || (FitMode==1))
      {   
         for(j=0;j<Length;j++) 
         PriceArray[j] = iMA(NULL,0,MA_Length,0,MA_Mode,Price,shift*FitMode + j);
      
         double Sum=0;
         for(j=Length-1;j>=0;j--)
         {
         double poly = PolyFit(PriceArray,Degree,Length,-j); 
         double del = PriceArray[j] - poly; 
         Sum += del*del;
         if (FitMode==0) Polynom[j] = poly;  
         }
      
      if(FitMode==1) Polynom[shift] = poly; 
      
         if (Length-1 > 0 && Sum > 0)
         {
            if(Length < 32) double StdDev = MathSqrt(Sum/(Length-1));   
            else
            StdDev = MathSqrt(Sum/Length);       
         }
      }
      
   UpBand[shift] = Polynom[shift] + K_Sigma * StdDev;    
   DnBand[shift] = Polynom[shift] - K_Sigma * StdDev;
      
      if(ProjLength > 0 && shift == 0)
      {
         for (j=0;j<=ProjLength;j++)
         Project[j] = PolyFit(PriceArray,Degree,Length,ProjLength - j);
      }  
   }
      
   for(shift=0;shift<limit;shift++ )
   {
   if (CenterMode == 1) PolyFitMA[shift] = Polynom[shift]; 
     
      if(ProjLength > 0 && shift == 0)
      {
         for (j=0;j<=ProjLength;j++)
         { 
         if (CenterMode == 1) Projection[j] = Project[j];
         UpProject[j] = Project[j] + K_Sigma * StdDev;   
         DnProject[j] = Project[j] - K_Sigma * StdDev;
         }
      }
   }     
   return(0);
}
//+---------------------------------------------------------------------------+