// correlation-matrix.afl
// Part 36 - Matrices and Advanced Data Handling
//
// Builds the correlation matrix of daily returns for a list of symbols using
// AFL's matrix type and the @ matrix-product operator.
//
// WHY A MATRIX
//   Correlation between every pair in a list is a two-dimensional answer.
//   With ordinary arrays you would need one named variable per pair, or a
//   nest of loops. Standardise each symbol's returns, stack them as the
//   columns of one Lookback x K matrix R, and the entire answer is a single
//   product: MxTranspose( R ) @ R.
//
// ASSUMPTIONS
//   - Daily end-of-day bars. No real-time feed, no Professional edition.
//   - Bars are matched by position in the window, so all symbols must share
//     the database's trading calendar. Mixing exchanges with different
//     holidays would silently pair up different dates.
//   - Correlation is measured over the last Lookback bars only. It describes
//     that window. It is not a stable property of the pair and says nothing
//     about what the next window will look like.
//   - A symbol with any Null return inside the window is dropped and named
//     in the report rather than being quietly treated as a zero return.
//   - Matrix support needs AmiBroker 6.00; MxSetBlock and MxTranspose need
//     6.10. Verified against AmiBroker 7.00.1.

_SECTION_BEGIN("Correlation Matrix");

SymbolList = "AAPL,MSFT,KO,XOM,JNJ";
Lookback   = Param( "Lookback bars", 120, 20, 500, 10 );

LastBar  = BarCount - 1;
FirstBar = BarCount - Lookback;   // first bar of the measurement window

Usable   = "";                    // comma-separated list of symbols that qualify
Rejected = "";
Cols     = 0;

if( BarCount > Lookback + 1 )
{
    // -----------------------------------------------------------------------
    // Pass 1 - standardise each symbol's returns over the window and keep the
    //          result under a name built from the ticker.
    // -----------------------------------------------------------------------
    for( i = 0; ( Sym = StrExtract( SymbolList, i ) ) != ""; i++ )
    {
        SetForeign( Sym );

        Ret   = ROC( Close, 1 );
        Mean  = MA( Ret, Lookback );
        Sigma = StDev( Ret, Lookback, False );   // sample SD, divides by n-1

        // Count Nulls inside the window with array processing rather than a
        // loop: Sum() of a Boolean array is a rolling count.
        NullsInWindow = Sum( IsNull( Ret ), Lookback );

        // Pre-dividing by sqrt(n-1) means the matrix product below IS the
        // correlation matrix, with no scalar division of a matrix afterwards.
        Standardised = ( Ret - Mean ) / ( Sigma * sqrt( Lookback - 1 ) );

        BadData = NullsInWindow[ LastBar ] > 0 OR Sigma[ LastBar ] <= 0;

        RestorePriceArrays();

        if( BadData )
        {
            Rejected = Rejected + Sym + " ";
        }
        else
        {
            VarSet( "z_" + Sym, Standardised );
            Usable = Usable + Sym + ",";
            Cols   = Cols + 1;
        }
    }
}

if( Cols >= 2 )
{
    // -----------------------------------------------------------------------
    // Pass 2 - load the standardised windows into the columns of a matrix.
    // -----------------------------------------------------------------------
    R = Matrix( Lookback, Cols, 0 );

    for( j = 0; ( Sym = StrExtract( Usable, j ) ) != ""; j++ )
    {
        Standardised = VarGet( "z_" + Sym );

        // MxSetBlock takes its values starting at element 0 of the array, so
        // the window has to be moved to the front of a scratch array first.
        // Copying Close is simply a way to allocate an array of the right
        // length; every element used is overwritten below.
        Block = Close;

        for( k = 0; k < Lookback; k++ )
            Block[ k ] = Standardised[ FirstBar + k ];

        R = MxSetBlock( R, 0, Lookback - 1, j, j, Block );
    }

    // The only true matrix product in AFL. Every other operator is
    // element-wise. R is Lookback x Cols, so the transpose is Cols x Lookback
    // and the product is Cols x Cols.
    Corr = MxTranspose( R ) @ R;

    Report = "Correlation of daily returns, last "
             + NumToStr( Lookback, 1.0 ) + " bars\n"
             + "Order: " + Usable + "\n"
             + "Matrix: " + NumToStr( MxGetSize( Corr, 0 ), 1.0 ) + " x "
             + NumToStr( MxGetSize( Corr, 1 ), 1.0 ) + "\n";

    if( Rejected != "" )
        Report = Report + "Dropped: " + Rejected + "\n";

    printf( "%s\n", Report );
    printf( "%s\n", MxToString( Corr ) );
}
else
{
    Report = "Not enough usable symbols. BarCount = "
             + NumToStr( BarCount, 1.0 )
             + ", window = " + NumToStr( Lookback, 1.0 ) + " bars.\n";

    if( Rejected != "" )
        Report = Report + "Dropped: " + Rejected + "\n";

    printf( "%s\n", Report );
}

Plot( Close, "Close", colorDefault, styleCandle );

Title = Report;

_SECTION_END();