// timeseries.ts - A library for time series analysis, focusing on ARIMA.

// ========================================
// TYPE DEFINITIONS
// ========================================

/**
 * Defines the parameters for an ARIMA model.
 * (p, d, q) are the non-seasonal components.
 * (P, D, Q, s) are the optional seasonal components for SARIMA.
 */
export interface ARIMAOptions {
  p: number; // AutoRegressive (AR) order
  d: number; // Differencing (I) order
  q: number; // Moving Average (MA) order
  P?: number; // Seasonal AR order
  D?: number; // Seasonal Differencing order
  Q?: number; // Seasonal MA order
  s?: number; // Seasonal period length
}

/**
 * The result object from an ARIMA forecast.
 */
export interface ARIMAForecastResult {
  forecast: number[]; // The predicted future values
  residuals: number[]; // The errors of the model fit on the original data
  model: ARIMAOptions; // The model parameters used
}

/**
 * The result object from an STL decomposition.
 */
export interface STLDecomposition {
  seasonal: number[]; // The seasonal component of the series
  trend: number[];    // The trend component of the series
  residual: number[]; // The remainder/residual component
  original: number[]; // The original series, for comparison
}


/**
 * A class for performing time series analysis, including identification and forecasting.
 */
export class TimeSeriesAnalyzer {

  // ========================================
  // 1. IDENTIFICATION METHODS
  // ========================================

  /**
   * Calculates the difference of a time series.
   * This is the 'I' (Integrated) part of ARIMA, used to make a series stationary.
   * @param series The input data series.
   * @param lag The lag to difference by (usually 1).
   * @returns A new, differenced time series.
   */
  static difference(series: number[], lag: number = 1): number[] {
    if (lag < 1 || !Number.isInteger(lag)) {
      throw new Error('Lag must be a positive integer.');
    }
    if (series.length <= lag) {
      return [];
    }

    const differenced: number[] = [];
    for (let i = lag; i < series.length; i++) {
      differenced.push(series[i] - series[i - lag]);
    }
    return differenced;
  }
  
  /**
   * Helper function to calculate the autocovariance of a series at a given lag.
   */
  private static autocovariance(series: number[], lag: number): number {
    const n = series.length;
    if (lag >= n) return 0;
    const mean = series.reduce((a, b) => a + b) / n;
    let sum = 0;
    for (let i = lag; i < n; i++) {
      sum += (series[i] - mean) * (series[i - lag] - mean);
    }
    return sum / n;
  }

  /**
   * Calculates the Autocorrelation Function (ACF) for a time series.
   * ACF helps in determining the 'q' parameter for an ARIMA model.
   * @param series The input data series.
   * @param maxLag The maximum number of lags to calculate.
   * @returns An array of correlation values from lag 1 to maxLag.
   */
  static calculateACF(series: number[], maxLag: number): number[] {
    if (series.length < 2) return [];
    
    const variance = this.autocovariance(series, 0);
    if (variance === 0) {
      return new Array(maxLag).fill(1);
    }

    const acf: number[] = [];
    for (let lag = 1; lag <= maxLag; lag++) {
      acf.push(this.autocovariance(series, lag) / variance);
    }
    return acf;
  }

  /**
   * Calculates the Partial Autocorrelation Function (PACF) for a time series.
   * This now uses the Durbin-Levinson algorithm for an accurate calculation.
   * PACF helps in determining the 'p' parameter for an ARIMA model.
   * @param series The input data series.
   * @param maxLag The maximum number of lags to calculate.
   * @returns An array of partial correlation values from lag 1 to maxLag.
   */
  static calculatePACF(series: number[], maxLag: number): number[] {
    const acf = this.calculateACF(series, maxLag);
    const pacf: number[] = [];
    
    if (acf.length === 0) return [];

    pacf.push(acf[0]); // PACF at lag 1 is the same as ACF at lag 1

    for (let k = 2; k <= maxLag; k++) {
      let numerator = acf[k - 1];
      let denominator = 1;

      const phi = new Array(k + 1).fill(0).map(() => new Array(k + 1).fill(0));
      
      for(let i=1; i<=k; i++) {
        phi[i][i] = acf[i-1];
      }

      for (let j = 1; j < k; j++) {
        const factor = pacf[j - 1];
        numerator -= factor * acf[k - j - 1];
        denominator -= factor * acf[j - 1];
      }
      
      if (Math.abs(denominator) < 1e-9) { // Avoid division by zero
        pacf.push(0);
        continue;
      }
      
      const pacf_k = numerator / denominator;
      pacf.push(pacf_k);
    }

    return pacf;
  }

  /**
   * Decomposes a time series using the robust Classical Additive method.
   * This version correctly isolates trend, seasonal, and residual components.
   * @param series The input data series.
   * @param period The seasonal period (e.g., 7 for daily data with a weekly cycle).
   * @returns An object containing the seasonal, trend, and residual series.
   */
  static stlDecomposition(series: number[], period: number): STLDecomposition {
    if (series.length < 2 * period) {
      throw new Error("Series must be at least twice the length of the seasonal period.");
    }
    
    // Helper for a centered moving average
    const movingAverage = (data: number[], window: number) => {
        const result = [];
        const halfWindow = Math.floor(window / 2);
        for (let i = 0; i < data.length; i++) {
            const start = Math.max(0, i - halfWindow);
            const end = Math.min(data.length, i + halfWindow + 1);
            let sum = 0;
            for (let j = start; j < end; j++) {
                sum += data[j];
            }
            result.push(sum / (end - start));
        }
        return result;
    };

    // Step 1: Calculate the trend using a centered moving average.
    // If period is even, we use a 2x-MA to center it correctly.
    let trend: number[];
    if (period % 2 === 0) {
        const intermediate = movingAverage(series, period);
        trend = movingAverage(intermediate, 2);
    } else {
        trend = movingAverage(series, period);
    }
    
    // Step 2: Detrend the series
    const detrended = series.map((val, i) => val - trend[i]);

    // Step 3: Calculate the seasonal component by averaging the detrended values for each period
    const seasonalAverages = new Array(period).fill(0);
    const seasonalCounts = new Array(period).fill(0);
    for (let i = 0; i < series.length; i++) {
        if (!isNaN(detrended[i])) {
            const seasonIndex = i % period;
            seasonalAverages[seasonIndex] += detrended[i];
            seasonalCounts[seasonIndex]++;
        }
    }
    
    for (let i = 0; i < period; i++) {
        seasonalAverages[i] /= seasonalCounts[i];
    }
    
    // Center the seasonal component to have a mean of zero
    const seasonalMean = seasonalAverages.reduce((a, b) => a + b, 0) / period;
    const centeredSeasonalAverages = seasonalAverages.map(avg => avg - seasonalMean);

    const seasonal = new Array(series.length).fill(0);
    for (let i = 0; i < series.length; i++) {
        seasonal[i] = centeredSeasonalAverages[i % period];
    }

    // Step 4: Calculate the residual component
    const residual = detrended.map((val, i) => val - seasonal[i]);

    return {
      original: series,
      seasonal,
      trend,
      residual,
    };
  }


  // ========================================
  // 2. FORECASTING METHODS
  // ========================================

  /**
   * [UPGRADED] Generates a forecast using a simplified SARIMA model.
   * This implementation now handles both non-seasonal (p,d,q) and seasonal (P,D,Q,s) components.
   * @param series The input time series data.
   * @param options The SARIMA parameters.
   * @param forecastSteps The number of future steps to predict.
   * @returns An object containing the forecast and model residuals.
   */
  static arimaForecast(series: number[], options: ARIMAOptions, forecastSteps: number): ARIMAForecastResult {
    const { p, d, q, P = 0, D = 0, Q = 0, s = 0 } = options;

    if (series.length < p + d + (P + D) * s + q + Q * s) {
      throw new Error("Data series is too short for the specified SARIMA order.");
    }
    
    const originalSeries = [...series];
    let differencedSeries = [...series];
    const diffLog: { lag: number, values: number[] }[] = [];

    // Step 1: Apply seasonal differencing 'D' times
    for (let i = 0; i < D; i++) {
        diffLog.push({ lag: s, values: differencedSeries.slice(-s) });
        differencedSeries = this.difference(differencedSeries, s);
    }
    
    // Step 2: Apply non-seasonal differencing 'd' times
    for (let i = 0; i < d; i++) {
        diffLog.push({ lag: 1, values: differencedSeries.slice(-1) });
        differencedSeries = this.difference(differencedSeries, 1);
    }

    const n = differencedSeries.length;
    // Simplified coefficients
    const arCoeffs = p > 0 ? new Array(p).fill(1 / p) : [];
    const maCoeffs = q > 0 ? new Array(q).fill(1 / q) : [];
    const sarCoeffs = P > 0 ? new Array(P).fill(1 / P) : [];
    const smaCoeffs = Q > 0 ? new Array(Q).fill(1 / Q) : [];
    
    const residuals: number[] = new Array(n).fill(0);
    const fitted: number[] = new Array(n).fill(0);

    // Step 3: Fit the model
    const startIdx = Math.max(p, q, P * s, Q * s);
    for (let t = startIdx; t < n; t++) {
        // Non-seasonal AR
        let arVal = 0;
        for (let i = 0; i < p; i++) arVal += arCoeffs[i] * differencedSeries[t - 1 - i];

        // Non-seasonal MA
        let maVal = 0;
        for (let i = 0; i < q; i++) maVal += maCoeffs[i] * residuals[t - 1 - i];
        
        // Seasonal AR
        let sarVal = 0;
        for (let i = 0; i < P; i++) sarVal += sarCoeffs[i] * differencedSeries[t - s * (i + 1)];

        // Seasonal MA
        let smaVal = 0;
        for (let i = 0; i < Q; i++) smaVal += smaCoeffs[i] * residuals[t - s * (i + 1)];
        
        fitted[t] = arVal + maVal + sarVal + smaVal;
        residuals[t] = differencedSeries[t] - fitted[t];
    }

    // Step 4: Generate the forecast
    const forecastDifferenced: number[] = [];
    const extendedSeries = [...differencedSeries];
    const extendedResiduals = [...residuals];

    for (let f = 0; f < forecastSteps; f++) {
        const t = n + f;
        let nextForecast = 0;
        
        // AR
        for (let i = 0; i < p; i++) nextForecast += arCoeffs[i] * extendedSeries[t - 1 - i];
        // MA (future residuals are 0)
        for (let i = 0; i < q; i++) nextForecast += maCoeffs[i] * extendedResiduals[t - 1 - i];
        // SAR
        for (let i = 0; i < P; i++) nextForecast += sarCoeffs[i] * extendedSeries[t - s * (i + 1)];
        // SMA
        for (let i = 0; i < Q; i++) nextForecast += smaCoeffs[i] * extendedResiduals[t - s * (i + 1)];

        forecastDifferenced.push(nextForecast);
        extendedSeries.push(nextForecast);
        extendedResiduals.push(0);
    }

    // Step 5: Invert the differencing
    let forecast = [...forecastDifferenced];
    for (let i = diffLog.length - 1; i >= 0; i--) {
        const { lag, values } = diffLog[i];
        const inverted = [];
        const fullHistory = [...originalSeries, ...forecast]; // Need a temporary full history for inversion
        
        // A simpler inversion method for forecasting
        let history = [...series];
        for (const forecastVal of forecast) {
            const lastSeasonalVal = history[history.length - lag];
            const invertedVal = forecastVal + lastSeasonalVal;
            inverted.push(invertedVal);
            history.push(invertedVal);
        }
        forecast = inverted;
    }

    return {
      forecast,
      residuals,
      model: options,
    };
  }
}