//+------------------------------------------------------------------+
//|                                                 PolyFitMA_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 2
#property indicator_color1  Orange
#property indicator_width1  2
#property indicator_color2  Orange
#property indicator_width2  1
#property indicator_style2  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    =  50;   //Number of bars for an evaluation
extern int Degree    =   4;   //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 =  0;   //Length of Projection
extern int FitMode   =   0;   //Fitting Mode: 0-Fitting,1-Moving
extern int CountBars =1000;   //Maximum Number of Bars(0-all Bars)
//---- indicator buffers
double PolyFitMA[];
double Project[];
//---- global variables
double PriceArray[];
int cBars;
//+---------------------------------------------------------------------------+
//| Custom indicator initialization function                                  |
//+---------------------------------------------------------------------------+
int init()
{
   IndicatorDigits(MarketInfo(Symbol(), MODE_DIGITS )+2);
   
   SetIndexStyle(0, DRAW_LINE);
   SetIndexStyle(1, DRAW_LINE);
   
   SetIndexBuffer(0, PolyFitMA);
   SetIndexBuffer(1, Project);
   
   if(CountBars == 0) cBars = Bars;
   else cBars = CountBars;
   
   SetIndexDrawBegin( 0,Bars-cBars+(Length+MA_Length));
      
   SetIndexShift(1,ProjLength);
   string short_name="PolyFitMA("+Price+","+Length+","+Degree+")";
   IndicatorShortName(short_name);
   SetIndexLabel(0,"PolyFitMA");
   SetIndexLabel(1,"Projection");
   
   ArrayResize(PriceArray, Length);
   
   return(0);
}

double PolyFit(double price[],int deg,int len,int bar)
{
double   result=0,Sum=0;
double   AX[13,13],BX[13],CX[13],ZX[13],Pow[26];
int      Row[13];  
// 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=0;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;
	   }
     
	if(AX[Row[deg+1],deg+1] != 0) 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]; //AX[Row[KX],Col]*CX[Col];
  	   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);
}	
//+------------------------------------------------------------------+
//| PolyFitMA_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++) 
      PolyFitMA[cBars-i]=EMPTY_VALUE;    
   }   
   
   if ( ProjLength > 0 )
   {  
      for(i=1;i<Length;i++)
      Project[i]=EMPTY_VALUE;  
   } 
      
   for(shift = 0;shift<limit;shift++ )
   {
      for(j = 0; j < Length; j++) 
      PriceArray[j] = iMA(NULL,0,MA_Length,0,MA_Mode,Price,shift*FitMode + j);
      
   PolyFitMA[shift] = PolyFit(PriceArray,Degree,Length,shift*(FitMode-1));
   
      if(ProjLength > 0 && shift == 0)
      {
         for (j=0;j<=ProjLength;j++)
         Project[j] = PolyFit(PriceArray,Degree,Length,ProjLength - j);
      }
   }
   return(0);
}
//+---------------------------------------------------------------------------+


