// Least Squares Polynomial Fit test
// comparing a least squares fit using matrix calculations
// with one that does not
 
SetBarsRequired(500,0);

order = Param( "n-th Order", 3, 1, 10, 1 );
nbar = Param( "Lookback nbars", 50, 10, 500, 1 );
mbck = Param( "Move Back", 0, 0, 100, 1 );
clevel = Param( "Confidence Level", 2, 1, 3, 0.1 );
extend = Param( "Extend Fit (Bars)", 10, 0, 20, 1 );

eb = BarCount - 1 - mbck;
bars = nbar;
prc = C;

SetChartOptions( 0, chartShowDates );
SetBarFillColor( IIf( C > O, ColorRGB( 0, 75, 0 ), IIf( C <= O, ColorRGB( 75, 0, 0 ), colorLightGrey ) ) );
Plot( C, "", IIf( C > O, ColorRGB( 0, 255, 0 ), IIf( C <= O, ColorRGB( 255, 0, 0 ), colorLightGrey ) ), styleCandle, Null, Null, 0, 0, 1 );

////////////////////////////////////// START FIT USING MATRIX CALC //////////////////////////////////////
function MxInvert( mm )
{
    n = rows = MxGetSize( mm, 0 );
    iden = MxIdentity( rows );
    inv = Matrix( rows, rows, 0 );

    for( i = 0; i < n; i++ )
    {
        for( j = 0; j < n; j++ )
        {
            if( i != j )
            {
                ratio = mm[j][i] / mm[i][i];

                for( k = 0; k < 2 * n; k++ )
                {
                    if( k < n )
                    {
                        mm[j][k] -= ratio * mm[i][k];
                    }
                    else
                        if( k >= n )
                        {
                            iden[j][k - n] -= ratio * iden[i][k - n];
                        }
                }
            }
        }
    }

    for( i = 0; i < n; i++ )
    {
        a = mm[i][i];

        for( j = 0; j < 2 * n; j++ )
        {
            if( j < n )
            {
                mm[i][j] /= a;
            }
            else
                if( j >= n )
                {
                    iden[i][j - n] /= a;
                }
        }
    }

    for( i = 0; i < n; i++ )
    {
        for( j = n; j < 2 * n; j++ )
        {
            if( j < n )
            {
                inv[i][j] = mm[i][j];
            }
            else
                if( j >= n )
                {
                    inv[i][j - n] = iden[i][j - n];
                }
        }
    }

    return inv;
}
function get_y( eb, bars )
{
    lb = eb;
    fb = eb - bars;
    rows = ( lb - fb + 1 );
    yy = Matrix( rows, 1, 0 );
    cnt = 0;

    for( i = fb; i <= lb; i++ )
    {
        yy[cnt][ 0 ] = prc[ i ];
        cnt = cnt + 1;
    }

    return yy;
}
function get_x( eb, bars )
{
    lb = eb;
    fb = eb - bars;
    rows = ( lb - fb + 1 );
    xx = Matrix( rows, order + 1, 1 );

    for( j = 1; j <= order; j++ )
    {
        cnt = 0;

        for( i = fb; i <= lb; i++ )
        {
            xx[cnt][j] = ( i - ( lb + fb ) / 2 ) ^ j;
            cnt = cnt + 1;
        }
    }

    return xx;
}
function calculateCoefficients( eb, bars )
{
    xx = get_x( eb, bars );
    yy = get_y( eb, bars );
    xxt = MxTranspose( xx );
    aa = MxInvert( xxt @ xx ) @ xxt @ yy;

    return aa;
}
function calculateFit( eb, bars )
{
    global rr;
    global xxa;
    global zza;
    global rrext;
    global rrextxxa;
    global rrextzxa;

    lb = eb;
    fb = eb - bars;

    aa = calculateCoefficients( eb, bars );

    rr = rrext = Null; // store the fit in rr

    for( i = fb; i <= lb; i++ )
    {
        rr[i] = aa[0][0];

        for( j = 1; j <= order; j++ )
        {
            rr[i] = rr[i] + aa[j][0] * ( i - ( lb + fb ) / 2 ) ^ j;
        }
    }

    // extended fit
    if( lb == eb )
    {
        for( i = lb + 1; i <= lb + extend; i++ )
        {
            rrext[i - extend] = aa[0][0];

            for( j = 1; j <= order; j++ )
            {
                vv = ( i - ( lb + fb ) / 2 );
                rrext[i - extend] = rrext[i - extend] + aa[j][0] * vv ^ j;
            }
        }
    }

    // calculate standard deviation
    sdp = 0;

    for( i = fb; i <= lb; i++ )
    {
        sdp = sdp + ( prc[i] - rr[i] ) ^ 2;
    }

    sd = sqrt( sdp / ( bars - 2 ) ); // devide by ( bars - 2 ) corresponding to StdErr function

    xxa = rr + sd * clevel;
    zza = rr - sd * clevel;
    rrextxxa = rrext + sd * clevel;
    rrextzxa = rrext - sd * clevel;
}

calculateFit( eb, bars );
Plot( rr, "", colorYellow, styleLine, Null, Null, 0, 0, 2 );
Plot( rrext, "", colorYellow, styleDashed | styleNoLabel | styleNoRescale, Null, Null, extend, 0, 2 );
Plot( xxa, "", colorYellow, styleline, null, null, 0, 0, 2 );
Plot( zza, "", colorYellow, styleline, null, null, 0, 0, 2 );
////////////////////////////////////// END FIT USING MATRIX CALC //////////////////////////////////////


////////////////////////////////////// START FIT USING OTHER METHOD //////////////////////////////////////
function CalculateCoefficients1( aa, bb )
{
    n = MxGetSize( aa, 0 );
    ll = uu = Matrix( n, n, 0 );
    xx = yy = 0;

    for( j = 0; j < n; j++ )
    {
        for( i = 0; i < n; i++ )
        {
            if( i <= j )
            {
                uu[i][j] = aa[i][j];

                for( k = 0; k <= i - 1; k++ )
                    uu[i][j] -= ll[i][k] * uu[k][j];

                if( i == j )
                    ll[i][j] = 1;
                else
                    ll[i][j] = 0;
            }
            else
            {
                ll[i][j] = aa[i][j];

                for( k = 0; k <= j - 1; k++ )
                    ll[i][j] -= ll[i][k] * uu[k][j];

                ll[i][j] /= uu[j][j];
                uu[i][j] = 0;
            }
        }
    }

    for( i = 0; i < n; i++ )
    {
        yy[i] = bb[i];

        for( j = 0; j < i; j++ )
        {
            yy[i] -= ll[i][j] * yy[j];
        }
    }

    for( i = n - 1; i >= 0; i-- )
    {
        xx[i] = yy[i];

        for( j = i + 1; j < n; j++ )
        {
            xx[i] -= uu[i][j] * xx[j];
        }

        xx[i] /= uu[i][i];
    }

    return xx;
}

function CalculateFit1( eb, bars )
{
    global reg;
    global x1a;
    global z1a;
    global regext;
    global regextx1a;
    global regextz1a;

    reg = x1a = z1a = Null;
    regext = regextx1a = regextz1a = Null;

    lb = eb;
    fb = eb - bars;
    nb = lb - fb;

    if( eb > bars )
    {
        aa = Matrix( order + 1, order + 1, 0 );
        bb = 0;

        // fill matrix A
        for( i = 0; i <= order; i++ )
        {
            for( j = 0; j <= order; j++ )
            {
                cnt = 0;

                for( k = fb; k <= lb; k++ )
                {
                    vv = ( k - ( lb + fb ) / 2 );
                    aa[i][j] = aa[i][j] + ( vv ^ ( i + j ) );
                    //aa[i][j] = aa[i][j] + ( cnt ^ ( i + j ) );
                    cnt = cnt + 1;
                }
            }
        }

        aa[0][0] = cnt;

        // fill matrix B
        for( i = 0; i <= order; i++ )
        {
            cnt = 0;

            for( j = fb; j <= lb; j++ )
            {
                vv = ( j - ( lb + fb ) / 2 );
                bb[i] = bb[i] + prc[j] * ( vv ^ i );
                //bb[i] = bb[i] + prc[j] * ( cnt ^ i );
                cnt = cnt + 1;
            }
        }

        // calculate coefficients
        xx = CalculateCoefficients1( aa, bb );

        // store the fit in reg
        cnt = 0;

        for( i = fb; i <= lb; i++ )
        {
            reg[i] = xx[0];

            for( j = 1; j <= order; j++ )
            {
                vv = ( i - ( lb + fb ) / 2 );
                reg[i] = reg[i] + xx[j] * vv ^ j;
                //reg[i] = reg[i] + xx[j] * cnt ^ j;
            }

            cnt = cnt + 1;
        }

        // extended fit (only when channel is active at last bar)
        if( lb == eb )
        {
            for( i = lb + 1; i <= lb + extend; i++ )
            {
                regext[i - extend] = xx[0];

                for( j = 1; j <= order; j++ )
                {
                    vv = ( i - ( lb + fb ) / 2 );
                    regext[i - extend] = regext[i - extend] + xx[j] * vv ^ j;
                    //regext[i - extend] = regext[i - extend] + xx[j] * cnt ^ j;
                }

                cnt = cnt + 1;
            }
        }

        // calculate standard deviation
        sdp = 0;

        for( i = fb; i <= lb; i++ )
        {
            sdp = sdp + ( prc[i] - reg[i] ) ^ 2;
        }

        sd = sqrt( sdp / ( bars - 2 ) ); // devide by ( bars - 2 ) corresponding to StdErr function

        x1a = reg + sd * clevel;
        z1a = reg - sd * clevel;
        regextx1a = regext + sd * clevel;
        regextz1a = regext - sd * clevel;
    }
}

calculateFit1( eb, bars );

Plot( reg, "", colorBlue, styleline, null, null, 0, 0, 8 );
Plot( regext, "", colorBlue, styleDashed | styleNoLabel | styleNoRescale, Null, Null, extend, 0, 8 );
Plot( x1a, "", colorBlue, styleline, null, null, 0, 0, 8 );
Plot( z1a, "", colorBlue, styleline, null, null, 0, 0, 8 );
////////////////////////////////////// END FIT USING OTHER METHOD //////////////////////////////////////