//+------------------------------------------------------------------+
//|                                                     polyMA.mq4 |
//|                                                   |
//|                                                                  |
//+------------------------------------------------------------------+
#property copyright ""
#property link      ""
// Include Neural Network package
#property indicator_chart_window
#property indicator_buffers 1
#property indicator_color1 Red
double Buffer1[];
extern int smoothing = 10;
int       endpos       = 0;               // Last value position
double    multStdDev   =1.96;            // Bands separation
bool      comments     =false;            // Comment on/off
// ---- buffers
double RegBfr[];
// ---- global variables
int      c_handle,p_handle,
current_day,previous_day;
double   errorstddev[4],
data[1],
regression[1],
regressionerror[1];
int init()
{
    IndicatorBuffers(1);
    SetIndexStyle(0,DRAW_LINE);
    SetIndexBuffer(0,Buffer1);
    ArrayResize(RegBfr, smoothing+1);
    //---- indicators
    //----
    return(0);
}
int start()
{
    int k = 0;
    int nCountedBars;
    nCountedBars=IndicatorCounted();
    double prevPrice = 0.0;
    bool   bFound;
    double dCurrent;
    int limit;
    int    counted_bars=IndicatorCounted();
    //---- check for possible errors
    if(counted_bars<0) return(-1);
    //---- last counted bar will be recounted
    if(counted_bars>0) counted_bars--;
    limit=Bars-counted_bars;
    //---- main loop
    if(limit > 2000)
      limit = 2000;
    for(int i=0; i<limit; i++)
    {
        ArrayResize(data,smoothing);
        ArrayResize(regression,smoothing);
        ArrayResize(regressionerror,smoothing);
        // ----
        ArrayInitialize(data,EMPTY_VALUE);
        for(k=0;k<smoothing;k++)
        {
         if(i+k == 0)
            data[k]=Open[i+k];
         else
            data[k]=(High[i+k]+Low[i+k])/2;
        }
        // ----
        ArrayInitialize(regression,EMPTY_VALUE);
        iPolReg(smoothing, 1, data, regression);
        ArrayInitialize(regressionerror,EMPTY_VALUE);
        for(k=0;k<smoothing;k++)
        regressionerror[k]=MathAbs(data[k]-regression[k]);
        errorstddev[0]=iStdDevOnArray(regressionerror,0,smoothing,0,MODE_SMA,0);
        // ----
        ArrayInitialize(regression,EMPTY_VALUE);
        iPolReg(smoothing, 2, data, regression);
        for(k=0;k<smoothing;k++)
        regressionerror[k]=MathAbs(data[k]-regression[k]);
        errorstddev[1]=iStdDevOnArray(regressionerror,0,smoothing,0,MODE_SMA,0);
        // ----
        ArrayInitialize(regression,EMPTY_VALUE);
        iLogReg(smoothing, data, regression);
        for(k=0;k<smoothing;k++)
        regressionerror[k]=MathAbs(data[k]-regression[k]);
        errorstddev[2]=iStdDevOnArray(regressionerror,0,smoothing,0,MODE_SMA,0);
        // ----
        ArrayInitialize(regression,EMPTY_VALUE);
        iExpReg(smoothing, data, regression);
        for(k=0;k<smoothing;k++)
        regressionerror[k]=MathAbs(data[k]-regression[k]);
        errorstddev[3]=iStdDevOnArray(regressionerror,0,smoothing,0,MODE_SMA,0);
        // ----
        BestReg();
        Buffer1[i] = RegBfr[0];
    }
    //----
    return(0);
}
void iPolReg(int datapoints, int grade, double& array1[], double& array2[])
{
    // y=c0+c1x+c2x^2+c3x^3...cnx^n
    // grade+1    n?mero de funciones y de coeficientes a calcular
    // 2*grade    valor del mayor exponente
    // datapoints n?mero de valores a considerar
    // pos        first value position (left to right direction)
    int     pos,exp,k,row,col,
    initialrow,initialcol,
    loop,i,j;
    double  sumx,sumyx,sum;
    double  sumxvalues[1],
    sumyxvalues[1],
    matrix[1][10],
    constant[1];
    pos=datapoints-1;
    ArrayResize(sumxvalues,2*grade+1);
    ArrayResize(sumyxvalues,grade+1);
    ArrayInitialize(sumxvalues,0.0);
    ArrayInitialize(sumyxvalues,0.0);
    sumxvalues[0]=datapoints;
    for(exp=1;exp<=2*grade;exp++)
    {
        sumx=0.0;
        sumyx=0.0;
        for(k=1;k<=datapoints;k++)
        {
            sumx+=MathPow(k,exp);
            if(exp==1)
            sumyx+=array1[pos-k+1];
            else if(exp<=grade+1)
            sumyx+=array1[pos-k+1]*MathPow(k,exp-1);
        }
        sumxvalues[exp]=sumx;
        if(sumyx!=0.0)
        sumyxvalues[exp-1]=sumyx;
    }
    ArrayResize(matrix,grade+1);
    ArrayInitialize(matrix,0.0);
    for(row=0;row<=grade;row++)
    for(col=0;col<=grade;col++)
    matrix[row][col]=sumxvalues[row+col];
    initialrow=1;
    initialcol=1;
    for(loop=1;loop<=grade;loop++)
    {
        for(row=initialrow;row<=grade;row++)
        {
            sumyxvalues[row]=sumyxvalues[row]-(matrix[row][loop-1]/matrix[loop-1][loop-1])*sumyxvalues[loop-1];
            for(col=initialcol;col<=grade;col++)
            matrix[row][col]=matrix[row][col]-(matrix[row][loop-1]/matrix[loop-1][loop-1])*matrix[loop-1][col];
        }
        initialrow++;
        initialcol++;
    }
    ArrayResize(constant,grade+1);
    ArrayInitialize(constant,0.0);
    j=0;
    for(i=grade;i>=0;i--)
    {
        if(j==0)
        constant[i]=sumyxvalues[i]/matrix[i][i];
        else
        {
            sum=0.0;
            for(k=j;k>0;k--)
            sum+=constant[i+k]*matrix[i][i+k];
            constant[i]=(sumyxvalues[i]-sum)/matrix[i][i];
        }
        j++;
    }
    ArrayResize(array2,datapoints);
    ArrayInitialize(array2,EMPTY_VALUE);
    k=1;
    for(i=datapoints-1;i>=0;i--)
    {
        sum=0.0;
        for(j=0;j<=grade;j++)
        sum+=constant[j]*MathPow(k,j);
        array2[i]=sum;
        k++;
    }
}
void iLogReg(int datapoints, double& array1[], double& array2[])
{
    // y=c0*x^c1
    // lny=lnc0+c1lnx <=> y=a+bx
    // c0=e^a
    // c1=b
    int     pos,exp,k,row,col,
    initialrow,initialcol,
    loop,i,j;
    double  sumx,sumyx,
    lnx,a;
    double  sumxvalues[1],
    sumyxvalues[1],
    matrix[1][10],
    constant[1];
    pos=datapoints-1;
    ArrayResize(sumxvalues,3);
    ArrayResize(sumyxvalues,2);
    ArrayInitialize(sumxvalues,0.0);
    ArrayInitialize(sumyxvalues,0.0);
    sumxvalues[0]=datapoints;
    for(exp=1;exp<=2;exp++)
    {
        sumx=0.0;
        sumyx=0.0;
        for(k=1;k<=datapoints;k++)
        {
            lnx=MathLog(k);
            sumx+=MathPow(lnx,exp);
            if(exp==1)
            sumyx+=MathLog(array1[pos-k+1]);
            else
            sumyx+=MathLog(array1[pos-k+1])*MathPow(lnx,exp-1);
        }
        sumxvalues[exp]=sumx;
        if(sumyx!=0.0)
        sumyxvalues[exp-1]=sumyx;
    }
    ArrayResize(matrix,2);
    ArrayInitialize(matrix,0.0);
    for(row=0;row<=1;row++)
    for(col=0;col<=1;col++)
    matrix[row][col]=sumxvalues[row+col];
    sumyxvalues[1]=sumyxvalues[1]-(matrix[1][0]/matrix[0][0])*sumyxvalues[0];
    matrix[1][1]=matrix[1][1]-(matrix[1][0]/matrix[0][0])*matrix[0][1];
    ArrayResize(constant,2);
    ArrayInitialize(constant,0.0);
    constant[1]=sumyxvalues[1]/matrix[1][1];
    a=(sumyxvalues[0]-(constant[1]*matrix[0][1]))/matrix[0][0];
    constant[0]=MathExp(a);
    ArrayResize(array2,datapoints);
    ArrayInitialize(array2,EMPTY_VALUE);
    k=1;
    for(i=datapoints-1;i>=0;i--)
    {
        array2[i]=constant[0]*MathPow(k,constant[1]);
        k++;
    }
}
void iExpReg(int datapoints, double& array1[], double& array2[])
{
    // y=c0*e^(xc1)
    // lny=lnc0+c1x <=> y=a+bx
    // c0=e^a
    // c1=b
    int     pos,exp,k,row,col,
    initialrow,initialcol,
    loop,i,j;
    double  sumx,sumyx,a;
    double  sumxvalues[1],
    sumyxvalues[1],
    matrix[1][10],
    constant[1];
    pos=datapoints-1;
    ArrayResize(sumxvalues,3);
    ArrayResize(sumyxvalues,2);
    ArrayInitialize(sumxvalues,0.0);
    ArrayInitialize(sumyxvalues,0.0);
    sumxvalues[0]=datapoints;
    for(exp=1;exp<=2;exp++)
    {
        sumx=0.0;
        sumyx=0.0;
        for(k=1;k<=datapoints;k++)
        {
            sumx+=MathPow(k,exp);
            if(exp==1)
            sumyx+=MathLog(array1[pos-k+1]);
            else
            sumyx+=MathLog(array1[pos-k+1])*MathPow(k,exp-1);
        }
        sumxvalues[exp]=sumx;
        if(sumyx!=0.0)
        sumyxvalues[exp-1]=sumyx;
    }
    ArrayResize(matrix,2);
    ArrayInitialize(matrix,0.0);
    for(row=0;row<=1;row++)
    for(col=0;col<=1;col++)
    matrix[row][col]=sumxvalues[row+col];
    sumyxvalues[1]=sumyxvalues[1]-(matrix[1][0]/matrix[0][0])*sumyxvalues[0];
    matrix[1][1]=matrix[1][1]-(matrix[1][0]/matrix[0][0])*matrix[0][1];
    ArrayResize(constant,2);
    ArrayInitialize(constant,0.0);
    constant[1]=sumyxvalues[1]/matrix[1][1];
    a=(sumyxvalues[0]-(constant[1]*matrix[0][1]))/matrix[0][0];
    constant[0]=MathExp(a);
    ArrayResize(array2,datapoints);
    ArrayInitialize(array2,EMPTY_VALUE);
    k=1;
    for(i=datapoints-1;i>=0;i--)
    {
        array2[i]=constant[0]*MathExp(k*constant[1]);
        k++;
    }
}
void BestReg()
{
    string  short_name;
    double  currentstddev;
    int     c_stddevindex,i;
    if(comments)
    Comment("LinReg=",errorstddev[0]," || QuadReg=",errorstddev[1],
    " || LogReg=",errorstddev[2]," || ExpReg=",errorstddev[3]);
    currentstddev=errorstddev[0];
    c_stddevindex=0;
    for(i=1;i<=3;i++)
    {
        if(errorstddev[i]<currentstddev && i == 5)
        {
            currentstddev=errorstddev[i];
            c_stddevindex=i;
        }
    }
    switch(c_stddevindex)
    {
        case 0 : iPolReg(smoothing, 1, data, regression); short_name="LinReg";  break;
        case 1 : iPolReg(smoothing, 2, data, regression); short_name="QuadReg"; break;
        case 2 : iLogReg(smoothing, data, regression);    short_name="LogReg";  break;
        case 3 : iExpReg(smoothing, data, regression);    short_name="ExpReg";  break;
        default :iPolReg(smoothing, 1, data, regression); short_name="LinReg";
    }
    SetIndexLabel(0,short_name);
    IndicatorShortName(short_name+"("+smoothing+")");
    ArrayInitialize(RegBfr,EMPTY_VALUE);
    for(i=smoothing-1;i>=0;i--)
    RegBfr[i+endpos]=regression[i];
}