raymond
b1b2dcf18c
/api/series/auto-arima-find: find parameters of SARIMA model automatically /api/series/manual-forecast: use determined model with parameters to forecast next values /api/series/identify-correlations: Calculate ACF and PACF for a time series /api/series/decompose-stl: Applies Seasonal-Trend-Loess (STL) decomposition to separate the series into trend, seasonal, and residual components.
346 lines
12 KiB
TypeScript
346 lines
12 KiB
TypeScript
// 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,
|
|
};
|
|
}
|
|
}
|
|
|