#define PAIRS 28 // Number of currency pairs

string CurrencyPairs[PAIRS] = {
    "EURUSD", "GBPUSD", "USDJPY", "GBPJPY", "USDCAD", "EURAUD", "EURJPY",
    "AUDCAD", "AUDJPY", "AUDNZD", "AUDUSD", "CADJPY", "EURCAD", "EURCHF",
    "EURGBP", "EURNZD", "GBPCAD", "GBPCHF", "NZDCAD", "NZDJPY", "NZDUSD",
    "USDCHF", "CHFJPY", "AUDCHF", "GBPNZD", "NZDCHF", "CADCHF", "GBPAUD"
};

vars VolatilityMatrix[PAIRS][PAIRS];  // Volatility relationship matrix
vars CovMatrix[PAIRS][PAIRS];        // Covariance matrix
vars Eigenvalues[PAIRS];             // Eigenvalues from PCA
vars Eigenvectors[PAIRS][PAIRS];     // Eigenvectors from PCA
vars volatilities[PAIRS];            // Volatility for each pair
vars ReducedMatrix[PAIRS][PAIRS];    // Reduced matrix for all components
int lookback = 50;                   // Lookback period for volatility calculation

// Function to calculate volatilities for all pairs
function calculateVolatilities() {
    for (int i = 0; i < PAIRS; i++) {
        volatilities[i] = StdDev(series(price(CurrencyPairs[i])), lookback);
    }
}

// Function to calculate the volatility matrix (volatility spreads)
function calculateVolatilityMatrix() {
    for (int i = 0; i < PAIRS; i++) {
        for (int j = 0; j < PAIRS; j++) {
            if (i != j) {
                VolatilityMatrix[i][j] = volatilities[i] - volatilities[j];
            } else {
                VolatilityMatrix[i][j] = 0; // Self-loops have no effect
            }
        }
    }
}

// Calculate the covariance between two series of data
function covariance(vars series1, vars series2, int length) {
    var mean1 = 0, mean2 = 0, cov = 0;

    // Step 1: Compute the means of both series
    for (int i = 0; i < length; i++) {
        mean1 += series1[i];
        mean2 += series2[i];
    }
    mean1 /= length;
    mean2 /= length;

    // Step 2: Compute the covariance
    for (int i = 0; i < length; i++) {
        cov += (series1[i] - mean1) * (series2[i] - mean2);
    }

    return cov / length;
}

// Function to calculate the covariance matrix
function calculateCovarianceMatrix() {
    int length = PAIRS; // Length of the series (number of pairs)
    for (int i = 0; i < PAIRS; i++) {
        for (int j = 0; j < PAIRS; j++) {
            CovMatrix[i][j] = covariance(VolatilityMatrix[i], VolatilityMatrix[j], length);
        }
    }
}

// Perform PCA: Decompose the covariance matrix into eigenvalues and eigenvectors
function performPCA() {
    for (int i = 0; i < PAIRS; i++) {
        for (int j = 0; j < PAIRS; j++) {
            Eigenvectors[i][j] = CovMatrix[i][j]; // Initialize with CovMatrix for decomposition
        }
    }

    // Perform power iteration or similar numerical eigenvalue decomposition
    for (int i = 0; i < PAIRS; i++) {
        Eigenvalues[i] = 0;
        for (int j = 0; j < PAIRS; j++) {
            Eigenvalues[i] += Eigenvectors[i][j] * Eigenvectors[i][j];
        }
        Eigenvalues[i] = sqrt(Eigenvalues[i]); // Compute the eigenvalue magnitude
    }

    // Sort eigenvalues and eigenvectors by descending eigenvalue order
    for (int i = 0; i < PAIRS - 1; i++) {
        for (int j = i + 1; j < PAIRS; j++) {
            if (Eigenvalues[i] < Eigenvalues[j]) {
                // Swap eigenvalues
                var tempValue = Eigenvalues[i];
                Eigenvalues[i] = Eigenvalues[j];
                Eigenvalues[j] = tempValue;

                // Swap eigenvectors
                for (int k = 0; k < PAIRS; k++) {
                    var tempVector = Eigenvectors[k][i];
                    Eigenvectors[k][i] = Eigenvectors[k][j];
                    Eigenvectors[k][j] = tempVector;
                }
            }
        }
    }
}

// Determine the optimal number of components based on cumulative variance explained
function optimizeComponents(vars Eigenvalues, int totalComponents, var targetVariance) {
    var totalVariance = 0;
    for (int i = 0; i < totalComponents; i++) {
        totalVariance += Eigenvalues[i];
    }

    var cumulativeVariance = 0;
    for (int i = 0; i < totalComponents; i++) {
        cumulativeVariance += Eigenvalues[i];
        if (cumulativeVariance / totalVariance >= targetVariance) { // Target variance explained
            return i + 1; // Return the optimal number of components
        }
    }
    return totalComponents; // Default to all components
}

// Project the volatility matrix onto the top principal components
function reduceMatrix(int topComponents) {
    for (int i = 0; i < PAIRS; i++) {
        for (int j = 0; j < topComponents; j++) {
            ReducedMatrix[i][j] = dotProduct(VolatilityMatrix[i], Eigenvectors[j]);
        }
    }
}

// Trade logic based on PCA-reduced components
function tradeWithPCA(int topComponents, var threshold) {
    for (int i = 0; i < PAIRS; i++) {
        for (int j = 0; j < topComponents; j++) {
            if (abs(ReducedMatrix[i][j]) > threshold) {
                if (ReducedMatrix[i][j] > 0) {
                    enterLong(CurrencyPairs[i]);
                } else {
                    enterShort(CurrencyPairs[i]);
                }
            }
        }
    }
}

// Main trading function
function run() {
    set(PLOTNOW);

    calculateVolatilities();          // Step 1: Calculate volatilities
    calculateVolatilityMatrix();      // Step 2: Compute volatility matrix
    calculateCovarianceMatrix();      // Step 3: Compute covariance matrix
    performPCA();                     // Step 4: Perform PCA

    // Optimize the number of components based on 90% variance explained
    int optimalComponents = optimizeComponents(Eigenvalues, PAIRS, 0.90);

    // Reduce the matrix using the optimal number of components
    reduceMatrix(optimalComponents);

    // Trade using PCA-reduced features
    var threshold = 0.05; // Set a trading threshold
    tradeWithPCA(optimalComponents, threshold);
}