// Downloaded From https://www.WiseStockTrader.com
// Center of Gravity (COG) indicator, original idea from El Mostafa Belkhayate
// Amibroker AFL code by E.M.Pottasch, 2011 
// Based on code by Fred Tonetti, 2006, n-th order Polynomial fit (see Amibroker Lib)
// JohnCW provided Gaussian_Eliminationsv function based on static variables.

//SetBarsRequired(sbrAll,sbrAll);
BI=BarIndex();
PF_EndBar=LastValue(BI);
PF_Y=(H+L)/2;
PF_Order=Param("nth Order",3,1,8,1);
PF_ExtraB=Param("Extrapolate Backwards",0,0,50,1);
PF_ExtraF=Param("Extrapolate Forwards",0,0,50,1);
Lookback=Param("Lookback Period",100,50,500,1);
sv=ParamToggle("Use Selected Value","Off|On",1);
norm=ParamToggle("Error Levels","Fibonacci|Normal",1);

if (sv)
{
	PF_EndBar=SelectedValue(bi);
	PF_BegBar=PF_EndBar-Lookback;
}
else
{
	PF_BegBar=PF_EndBar-Lookback;
}
function D2Set(L_value,i,j,L_name)
{
	local L_value,L_name,i,j;
	StaticVarSet(L_name + ":" + i + "," + j, L_value); 
}
function D2Get(i,j,L_name)
{
	local L_name,i,j;
	return(Nz(StaticVarGet(L_name + ":" + i + "," + j),0));
}
function Gaussian_Eliminationsv(GE_Order,GE_N,GE_SumXn,GE_SumYXn)
{
    w=0;Coeff=0;
    n=GE_Order+1;
	for(i=1;i<=n;i++ )
	{
		for(j=1;j<=n;j++)
		{
            if (i==1 AND j==1)
                D2Set(GE_N,i,j,"b");
            else
                D2Set(GE_SumXn[i+j-2],i,j,"b");
        }
        w[i]=GE_SumYXn[i];
    }
    n1=n-1;
    for(i=1;i<=n1;i++)
    {
        big=abs(D2Get(i,i,"b"));
        q=i;
        i1=i+1;
        for(j=i1;j<=n;j++)
        {
            ab=abs(D2Get(j,i,"b"));
            if(ab>=big)
            {
                big=ab;
                q=j;
            }
        }
        if (big!=0)
        {
            if (q!=i)
            {
                for (j=1;j<=n;j++)
                {
                    Temp=D2Get(q,j,"b");
                    D2Set(D2Get(i,j,"b"),q,j,"b");
                    D2Set(Temp,i,j,"b");
                }
                Temp=w[i];
                w[i]=w[q];
                w[q]=Temp;
            }
        }
        for(j=i1;j<=n;j++)
        {
            t=D2Get(j,i,"b")/D2Get(i,i,"b");
            for(k=i1;k<=n;k++)
            {
                D2Set(D2Get(j,k,"b")-t*D2Get(i,k,"b"),j,k,"b");
            }
            w[j]=w[j]-t*w[i];
        }
    }
    if(D2Get(n,n,"b")!=0)
    {
        Coeff[n]=w[n]/D2Get(n,n,"b");
        i=n-1;
        while(i>0)
        {
            SumY=0;
            i1=i+1;
            for(j=i1;j<=n;j++)
            {
                SumY=SumY+D2Get(i,j,"b")*Coeff[j];
            }
            Coeff[i]=(w[i]-SumY)/D2Get(i,i,"b");
            i=i-1;
        }
    }
    return Coeff;
}
function PolyFit(GE_Y,GE_BegBar,GE_EndBar,GE_Order,GE_ExtraB,GE_ExtraF)
{
	BI=BarIndex();
	GE_N=GE_EndBar-GE_BegBar+1;
	GE_XBegin=-(GE_N-1)/2;
	GE_X=IIf(BI<GE_BegBar,0,IIf(BI>GE_EndBar,0,(GE_XBegin+BI-GE_BegBar)));
	GE_X_Max=LastValue(Highest(GE_X));
	GE_X=GE_X/GE_X_Max;
	X1=GE_X;
	GE_Y=IIf(BI<GE_BegBar,0,IIf(BI>GE_EndBar,0,GE_Y));
	GE_SumXn=Cum(0);
	
	GE_SumXn[1]=LastValue(Cum(GE_X));
	GE_X2=GE_X*GE_X;GE_SumXn[2]=LastValue(Cum(GE_X2));
	GE_X3=GE_X*GE_X2;GE_SumXn[3]=LastValue(Cum(GE_X3));
	GE_X4=GE_X*GE_X3;GE_SumXn[4]=LastValue(Cum(GE_X4));
	GE_X5=GE_X*GE_X4;GE_SumXn[5]=LastValue(Cum(GE_X5));
	GE_X6=GE_X*GE_X5;GE_SumXn[6]=LastValue(Cum(GE_X6));
	GE_X7=GE_X*GE_X6;GE_SumXn[7]=LastValue(Cum(GE_X7));
	GE_X8=GE_X*GE_X7;GE_SumXn[8]=LastValue(Cum(GE_X8));
	GE_X9=GE_X*GE_X8;GE_SumXn[9]=LastValue(Cum(GE_X9));
	GE_X10=GE_X*GE_X9;GE_SumXn[10]=LastValue(Cum(GE_X10));
	GE_X11=GE_X*GE_X10;GE_SumXn[11]=LastValue(Cum(GE_X11));
	GE_X12=GE_X*GE_X11;GE_SumXn[12]=LastValue(Cum(GE_X12));
	GE_X13=GE_X*GE_X12;GE_SumXn[13]=LastValue(Cum(GE_X13));
	GE_X14=GE_X*GE_X13;GE_SumXn[14]=LastValue(Cum(GE_X14));
	GE_X15=GE_X*GE_X14;GE_SumXn[15]=LastValue(Cum(GE_X15));
	GE_X16=GE_X*GE_X15;GE_SumXn[16]=LastValue(Cum(GE_X16));
	
	GE_SumYXn=Cum(0);
	GE_SumYXn[1]=LastValue(Cum(GE_Y));
	GE_YX=GE_Y*GE_X;GE_SumYXn[2]=LastValue(Cum(GE_YX));
	GE_YX2=GE_YX*GE_X;GE_SumYXn[3]=LastValue(Cum(GE_YX2));
	GE_YX3=GE_YX2*GE_X;GE_SumYXn[4]=LastValue(Cum(GE_YX3));
	GE_YX4=GE_YX3*GE_X;GE_SumYXn[5]=LastValue(Cum(GE_YX4));
	GE_YX5=GE_YX4*GE_X;GE_SumYXn[6]=LastValue(Cum(GE_YX5));
	GE_YX6=GE_YX5*GE_X;GE_SumYXn[7]=LastValue(Cum(GE_YX6));
	GE_YX7=GE_YX6*GE_X;GE_SumYXn[8]=LastValue(Cum(GE_YX7));
	GE_YX8=GE_YX7*GE_X;GE_SumYXn[9]=LastValue(Cum(GE_YX8));
	
	GE_Coeff=Cum(0);
	
	GE_Coeff=Gaussian_Eliminationsv(GE_Order,GE_N,GE_SumXn,GE_SumYXn);
	
	for (i = 1; i <= GE_Order + 1; i++) printf(NumToStr(i, 1.0) + " = " + NumToStr(GE_Coeff[i], 1.9) + "\n");

	GE_X=IIf(BI<GE_BegBar-GE_ExtraB-GE_ExtraF,0,IIf(BI>GE_EndBar,0,(GE_XBegin+BI-GE_BegBar+GE_ExtraF)/GE_X_Max));
	
	GE_X2=GE_X*GE_X;GE_X3=GE_X2*GE_X;GE_X4=GE_X3*GE_X;GE_X5=GE_X4*GE_X;GE_X6=GE_X5*GE_X;
	GE_X7=GE_X6*GE_X;GE_X8=GE_X7*GE_X;GE_X9=GE_X8*GE_X;GE_X10=GE_X9*GE_X;GE_X11=GE_X10*GE_X; 
	GE_X12=GE_X11*GE_X;GE_X13=GE_X12*GE_X;GE_X14=GE_X13*GE_X;GE_X15=GE_X14*GE_X;GE_X16=GE_X15*GE_X; 
	
	GE_Yn=IIf(BI<GE_BegBar-GE_ExtraB-GE_ExtraF,-1e10,IIf(BI>GE_EndBar,-1e10,GE_Coeff[1]+
	GE_Coeff[2]*GE_X+GE_Coeff[3]*GE_X2+GE_Coeff[4]*GE_X3+GE_Coeff[5]*GE_X4+GE_Coeff[6]*GE_X5+
	GE_Coeff[7]*GE_X6+GE_Coeff[8]*GE_X7+GE_Coeff[9]*GE_X8));
	
	return GE_Yn;
}
	
Yn=PolyFit(PF_Y,PF_BegBar,PF_EndBar,PF_Order,PF_ExtraB,PF_ExtraF);

SetChartOptions(0, chartShowDates);
Title = "Symbol: "+ Name()+ "\nPoly Order: "+PF_Order;
Plot(C, "Close",colorLightGrey,styleCandle);
Plot(Yn,"",IIf(BI>PF_EndBar-PF_ExtraF,colorWhite,IIf(BI<PF_BegBar-PF_ExtraF,colorWhite,colorBlue)),styleThick,Null,Null,PF_ExtraF);

if(norm)
{
	se=StdErr((C-Yn),LookBack);se=se[PF_EndBar];
	//se=StDev(C,LookBack);se=se[PF_EndBar];
	seh2=Yn+ValueWhen(Yn,se*2);
	sel2=Yn-ValueWhen(Yn,se*2);
	seh1=Yn+ValueWhen(Yn,se*1);
	sel1=Yn-ValueWhen(Yn,se*1);
	Plot(seh2,"",IIf(BI>PF_EndBar-PF_ExtraF,colorWhite,IIf(BI<PF_BegBar-PF_ExtraF,colorWhite,ColorRGB(255,0,0))),styleThick,Null,Null,PF_ExtraF);
	Plot(sel2,"",IIf(BI>PF_EndBar-PF_ExtraF,colorWhite,IIf(BI<PF_BegBar-PF_ExtraF,colorWhite,ColorRGB(0,255,0))),styleThick,Null,Null,PF_ExtraF);
	Plot(seh1,"",IIf(BI>PF_EndBar-PF_ExtraF,colorWhite,IIf(BI<PF_BegBar-PF_ExtraF,colorWhite,ColorRGB(255,100,100))),styleDashed,Null,Null,PF_ExtraF);
	Plot(sel1,"",IIf(BI>PF_EndBar-PF_ExtraF,colorWhite,IIf(BI<PF_BegBar-PF_ExtraF,colorWhite,ColorRGB(100,255,100))),styleDashed,Null,Null,PF_ExtraF);
}
else
{
	se=StDev(C,LookBack);se=se[PF_EndBar];
	r1=(1+5^0.5)/2;
	se=se*r1;
	seh3=Yn+ValueWhen(Yn,se);
	sel3=Yn-ValueWhen(Yn,se);
	seh2=Yn+ValueWhen(Yn,se/(1.382));
	sel2=Yn-ValueWhen(Yn,se/(1.382));	
	seh1=Yn+ValueWhen(Yn,se/(1.382*1.618));
	sel1=Yn-ValueWhen(Yn,se/(1.382*1.618));	
	
	Plot(seh3,"",IIf(BI>PF_EndBar-PF_ExtraF,colorWhite,IIf(BI<PF_BegBar-PF_ExtraF,colorWhite,ColorRGB(255,0,0))),styleThick,Null,Null,PF_ExtraF);
	Plot(sel3,"",IIf(BI>PF_EndBar-PF_ExtraF,colorWhite,IIf(BI<PF_BegBar-PF_ExtraF,colorWhite,ColorRGB(0,255,0))),styleThick,Null,Null,PF_ExtraF);
	Plot(seh2,"",IIf(BI>PF_EndBar-PF_ExtraF,colorWhite,IIf(BI<PF_BegBar-PF_ExtraF,colorWhite,ColorRGB(255,100,100))),styleDashed,Null,Null,PF_ExtraF);
	Plot(sel2,"",IIf(BI>PF_EndBar-PF_ExtraF,colorWhite,IIf(BI<PF_BegBar-PF_ExtraF,colorWhite,ColorRGB(100,255,100))),styleDashed,Null,Null,PF_ExtraF);
	Plot(seh1,"",IIf(BI>PF_EndBar-PF_ExtraF,colorWhite,IIf(BI<PF_BegBar-PF_ExtraF,colorWhite,ColorRGB(255,200,200))),styleDashed,Null,Null,PF_ExtraF);
	Plot(sel1,"",IIf(BI>PF_EndBar-PF_ExtraF,colorWhite,IIf(BI<PF_BegBar-PF_ExtraF,colorWhite,ColorRGB(200,255,200))),styleDashed,Null,Null,PF_ExtraF);
}