ARIMA_PLUS Forecasting Framework

• ARIMA_PLUS is a modular decomposition framework that separates trend, holiday, seasonal, and anomalous components for interpretable time series analysis.
• It preprocesses data through regularization, interpolation, holiday extraction, and spike cleaning before applying ARIMA to model the residual trend.
• The framework automates model selection and scaling in BigQuery ML, achieving competitive accuracy metrics (MASE, RMSE, sMAPE) on large datasets.

ARIMA_PLUS is a modular decomposition framework for univariate, and in ARIMA_PLUS_XREG form ARIMAX-style multivariate, time series forecasting and anomaly detection, implemented as a fully managed service inside Google BigQuery. In the formulation introduced for BigQuery ML, it combines a sequential pipeline for preprocessing, holiday extraction, spikes and dips cleaning, multi-seasonality extraction, step-change adjustment, and residual trend modeling by non-seasonal ARIMA with in-database automation for model selection, forecasting, and anomaly detection. Its defining emphasis is the joint pursuit of large-scale automation, interpretability, and accuracy: the same fitted object exposes decomposed components, prediction intervals, anomaly scores, and SQL-native interfaces for large collections of time series (Cheng et al., 28 Oct 2025).

1. Definition and formal structure

ARIMA_PLUS treats a time series as a sum of interpretable components estimated sequentially. For historical data, the framework writes

$$\begin{aligned} \text{history\_time\_series} &= \text{trend} + \text{holiday\_effect} \ &\quad + \text{spikes\_and\_dips} + \text{step\_changes} + \text{residual} \ &\quad + \sum_{m\in\{\text{daily, weekly, monthly, quarterly, yearly}\}} \text{seasonality}_m \ &\quad + \sum_{i=1}^N \text{attribution\_feature}_i . \end{aligned}$$

For forecasts, the framework writes

$$\begin{aligned} \text{forecasted\_time\_series} &= \text{trend} + \text{holiday\_effect} \ &\quad + \sum_{m\in\{\text{daily, weekly, monthly, quarterly, yearly}\}} \text{seasonality}_m \ &\quad + \sum_{i=1}^N \text{attribution\_feature}_i . \end{aligned}$$

This decomposition is implemented through nine main steps: preprocessing or regularization, missing value interpolation, holiday extraction, spikes and dips cleaning, seasonality extraction, step-change adjustment, trend modeling by ARIMA, seasonality re-addition, and holiday re-addition. The sequential structure is central: each module fits a specific component, removes it, and passes a cleaned series to the next stage, so ARIMA is applied only after other systematic effects have been extracted (Cheng et al., 28 Oct 2025).

Stage	Operation	Primary output
1	Preprocessing / regularization	Equally spaced series
2	Missing value interpolation	Completed series
3	Holiday extraction	holiday_effect
4	Spikes and dips cleaning	spikes_and_dips
5	Seasonality extraction	seasonality_m
6	Step-change adjustment	step_changes
7	Trend modeling	ARIMA trend
8	Seasonality re-addition	Forecast seasonality
9	Holiday re-addition	Final forecast

The modularization is not merely an implementation convenience. It is the basis of the framework’s interpretability, since holidays, anomalies, seasonalities, structural changes, trend, and covariate attributions are represented as distinct outputs rather than latent effects folded into one monolithic forecasting model.

2. Data processing, holidays, and transient anomalies

The initial stages regularize timestamped observations into an equally spaced series. Timestamps are converted to UTC, malformed timestamps are cleaned, the series is sorted and deduplicated, and the frequency is detected from the mode of inter-arrival times. A regular grid is then constructed, with special handling for calendar frequencies such as monthly and quarterly data by mapping them to a month index so that unequal month lengths do not break equal-spacing assumptions. Time zones and daylight saving time are handled by working in a local time index and mapping back to UTC afterward (Cheng et al., 28 Oct 2025).

Missingness is treated as a first-class modeling issue. ARIMA_PLUS distinguishes native missing values, special numeric codes treated as missing, and periods explicitly marked as outliers. The interpolation layer supports linear interpolation, splines, LOESS, and Kriging, and introduces two additional principles. First, interpolation is seasonality-aware so that later seasonal estimation is not distorted. Second, for outlier periods it is noise-aware: if an outlier period $O$ is replaced, the framework computes a local smooth $\widetilde{Y}_i$ and restores local noise through

$$Y_i^{\text{final}} = (Y_i - \widetilde{Y}_i) + \hat{Y}_i, \quad i\in O.$$

Holiday modeling is performed before spike detection so that holiday-driven deviations are not misclassified as anomalies. The system supports fixed holidays, moving holidays, multi-region holiday tables, and user-defined custom holidays with names, dates, and pre- and post-holiday windows. A notable design feature is the use of “sub-holidays”: each position inside a holiday window is treated as a distinct sub-holiday, allowing effects to vary sharply within a holiday span rather than being forced into one aggregate effect. Overlapping holidays are reconciled by taking at most one positive and one negative effect, specifically the maximum positive and minimum negative estimated contributions, rather than summing all overlapping effects (Cheng et al., 28 Oct 2025).

Temporary anomalies are modeled through the spikes-and-dips module. After holiday removal, the series is passed through forward and backward double exponential smoothing. Points are retained as spikes or dips only when they are flagged in both directions and have the same deviation sign in the forward and backward smoothers, with additional edge handling near the beginning and end of the series. Repeated seasonal spikes occurring every year at the same position are excluded from the anomaly set. Detected spike and dip points are then set to missing and reconstructed via the interpolation layer, and the difference between original and reconstructed values becomes the spikes_and_dips component (Cheng et al., 28 Oct 2025).

3. Multi-seasonality, structural changes, and exogenous regressors

ARIMA_PLUS does not use seasonal ARIMA. Instead, seasonality is handled by an enhanced STL-based decomposition designed for multiple seasonal periods, irregular seasonal periods, dynamic seasonal amplitudes, and seasonality-over-seasonality effects. For daily-frequency data, the candidate periods include weekly, monthly, quarterly, and yearly effects; more generally, the system fits STL iteratively over candidate periods, tests significance, deduplicates overlapping seasonal patterns, and retains significant components as separate $\text{seasonality}_m$ outputs (Cheng et al., 28 Oct 2025).

A distinctive case is irregular monthly behavior, such as end-of-month effects in daily or hourly data. In this setting, each calendar month is rescaled to a standardized month of length 30 days, or $30 \times 24$ for hourly data, seasonal structure is estimated in the standardized domain, and the result is interpolated back to original timestamps. This permits monthly seasonality extraction even though actual months have 28 to 31 days.

Dynamic seasonality is forecast through what the framework calls “forecasting by leaping.” If $S_t$ is a seasonal component with period $K$, the method forms $K$ sub-series,

$$\begin{aligned} \text{forecasted\_time\_series} &= \text{trend} + \text{holiday\_effect} \ &\quad + \sum_{m\in\{\text{daily, weekly, monthly, quarterly, yearly}\}} \text{seasonality}_m \ &\quad + \sum_{i=1}^N \text{attribution\_feature}_i . \end{aligned}$$0

forecasts each sub-series separately using a trend model such as double exponential smoothing, and then stitches the sub-series forecasts back together. This allows, for example, a weekly seasonal pattern to evolve gradually over time instead of remaining fixed.

Longer-lived regime changes are handled separately from transient spikes. Change-point detection is based on local forward forecasts, local backward forecasts on the reversed series, and residual-based z-scores. Candidate forward change starts and backward change ends are paired into change periods, merged when necessary, and then passed to a Chow-style adjustment stage. Within each identified change interval, values are set to missing, stable periods before and after the interval are compared, and a change-adjustment series is constructed so that

$$\begin{aligned} \text{forecasted\_time\_series} &= \text{trend} + \text{holiday\_effect} \ &\quad + \sum_{m\in\{\text{daily, weekly, monthly, quarterly, yearly}\}} \text{seasonality}_m \ &\quad + \sum_{i=1}^N \text{attribution\_feature}_i . \end{aligned}$$1

The result is a cleaned history in which historical level and slope changes are neutralized before ARIMA fitting, while the estimated change adjustments are retained as an interpretable step_changes component (Cheng et al., 28 Oct 2025).

For exogenous regressors, ARIMA_PLUS_XREG uses regression with ARIMA errors. The covariate stage fits ridge regression,

$$\begin{aligned} \text{forecasted\_time\_series} &= \text{trend} + \text{holiday\_effect} \ &\quad + \sum_{m\in\{\text{daily, weekly, monthly, quarterly, yearly}\}} \text{seasonality}_m \ &\quad + \sum_{i=1}^N \text{attribution\_feature}_i . \end{aligned}$$2

computes fitted covariate contribution, subtracts it from the target, and then applies the univariate ARIMA_PLUS pipeline to the residual series. During forecasting, future covariate values are mapped through the fitted ridge model and added to the forecasted residual component. Categorical regressors are dummy-encoded automatically (Cheng et al., 28 Oct 2025).

4. ARIMA trend layer, forecasting, and anomaly detection

After holidays, spikes and dips, seasonalities, step changes, and—when present—covariate contributions are removed, the remaining signal is modeled by non-seasonal ARIMA:

$$\begin{aligned} \text{forecasted\_time\_series} &= \text{trend} + \text{holiday\_effect} \ &\quad + \sum_{m\in\{\text{daily, weekly, monthly, quarterly, yearly}\}} \text{seasonality}_m \ &\quad + \sum_{i=1}^N \text{attribution\_feature}_i . \end{aligned}$$3

where $$\begin{aligned} \text{forecasted\_time\_series} &= \text{trend} + \text{holiday\_effect} \ &\quad + \sum_{m\in\{\text{daily, weekly, monthly, quarterly, yearly}\}} \text{seasonality}_m \ &\quad + \sum_{i=1}^N \text{attribution\_feature}_i . \end{aligned}$$4 is the cleaned residual series. The auto-ARIMA procedure first applies repeated KPSS tests to determine the differencing order $$\begin{aligned} \text{forecasted\_time\_series} &= \text{trend} + \text{holiday\_effect} \ &\quad + \sum_{m\in\{\text{daily, weekly, monthly, quarterly, yearly}\}} \text{seasonality}_m \ &\quad + \sum_{i=1}^N \text{attribution\_feature}_i . \end{aligned}$$5, then searches over $$\begin{aligned} \text{forecasted\_time\_series} &= \text{trend} + \text{holiday\_effect} \ &\quad + \sum_{m\in\{\text{daily, weekly, monthly, quarterly, yearly}\}} \text{seasonality}_m \ &\quad + \sum_{i=1}^N \text{attribution\_feature}_i . \end{aligned}$$6 and $$\begin{aligned} \text{forecasted\_time\_series} &= \text{trend} + \text{holiday\_effect} \ &\quad + \sum_{m\in\{\text{daily, weekly, monthly, quarterly, yearly}\}} \text{seasonality}_m \ &\quad + \sum_{i=1}^N \text{attribution\_feature}_i . \end{aligned}$$7 subject to $$\begin{aligned} \text{forecasted\_time\_series} &= \text{trend} + \text{holiday\_effect} \ &\quad + \sum_{m\in\{\text{daily, weekly, monthly, quarterly, yearly}\}} \text{seasonality}_m \ &\quad + \sum_{i=1}^N \text{attribution\_feature}_i . \end{aligned}$$8, optionally includes linear drift when $$\begin{aligned} \text{forecasted\_time\_series} &= \text{trend} + \text{holiday\_effect} \ &\quad + \sum_{m\in\{\text{daily, weekly, monthly, quarterly, yearly}\}} \text{seasonality}_m \ &\quad + \sum_{i=1}^N \text{attribution\_feature}_i . \end{aligned}$$9, chooses the model minimizing AIC, and fits candidates by maximum likelihood with L-BFGS. By design, the framework restricts ARIMA to trend and residual modeling, leaving seasonality to STL and event effects to dedicated modules (Cheng et al., 28 Oct 2025).

Forecasts are then reconstructed by re-adding future seasonal and holiday components to the ARIMA trend forecast, and, in the XREG variant, the future covariate contribution. The same decomposition also underlies anomaly detection. The expected historical value is defined as the sum of trend, holiday effect, seasonality, and step changes, excluding spikes and dips and residual noise. Residuals are treated as Gaussian with standard deviation $O$0 estimated from the ARIMA stage, and anomaly probability is defined by

$O$1

where $O$2 is the standard normal CDF. For a chosen threshold $O$3, the interval size is

$O$4

yielding the pointwise interval

$O$5

A point is flagged if its observed value lies outside this interval or equivalently if its anomaly probability exceeds the threshold. This unifies forecasting and anomaly detection in a single component-based model rather than treating them as separate tasks (Cheng et al., 28 Oct 2025).

5. In-database implementation, scale, and empirical performance

ARIMA_PLUS is exposed in BigQuery ML as MODEL_TYPE = 'ARIMA_PLUS' for univariate series and MODEL_TYPE = 'ARIMA_PLUS_XREG' for series with exogenous regressors. Forecasting is accessed through ML.FORECAST; decomposition through ML.EXPLAIN_FORECAST; anomaly detection through ML.DETECT_ANOMALIES; evaluation through ML.EVALUATE; and ARIMA coefficients through ML.ARIMA_COEFFICIENTS. A single CREATE MODEL statement can train many series simultaneously, with TIME_SERIES_ID_COL defining independent series-specific pipelines executed in parallel (Cheng et al., 28 Oct 2025).

The system automates frequency detection, data regularization, missing value handling, holiday lookup through HOLIDAY_REGION, spike and dip cleaning, change-point detection, seasonality selection, and auto-ARIMA configuration. It also exposes user controls such as AUTO_ARIMA_MAX_ORDER, NON_SEASONAL_ORDER, INCLUDE_DRIFT, SEASONALITIES, CLEAN_SPIKES_AND_DIPS, ADJUST_STEP_CHANGES, FORECAST_LIMIT_LOWER_BOUND, FORECAST_LIMIT_UPPER_BOUND, and time-series-length controls including MIN_TIME_SERIES_LENGTH and MAX_TIME_SERIES_LENGTH. Defaults were tuned empirically; for example, AUTO_ARIMA_MAX_ORDER = 2 and MAX_TIME_SERIES_LENGTH = 1024 are reported as effective accuracy–speed compromises (Cheng et al., 28 Oct 2025).

The infrastructure claims are unusually central to the identity of ARIMA_PLUS. The framework is directly built into the query engine of BigQuery in Google Cloud, uses a simple SQL interface, and automates technical steps such as data cleaning and model selection. The paper reports that it is possible to forecast 100 million time series using only 1.5 hours with a throughput of more than 18000 time series per second (Cheng et al., 28 Oct 2025).

Empirically, the framework is benchmarked on the Monash forecasting repository. On 41 datasets with available metrics, it achieves the lowest geometric mean MASE, beating the second-best model by 17%, and the best average-rank MASE, approximately 8.5 versus approximately 9.6 for the second-best model. It is reported as top-1 on 11 datasets and top-5 on 19 datasets. The comparison set includes established statistical models such as ETS, ARIMA, TBATS, and Prophet as well as newer neural models such as DeepAR, N-BEATS, PatchTST, and TimeMixer. In the ARIMA_PLUS_XREG case study for PM2.5 forecasting, adding wind speed and temperature as regressors improves MAE by 8.2%, MSE by 32%, RMSE by 14.6%, and sMAPE by 7.5% relative to the univariate version (Cheng et al., 28 Oct 2025).

6. Interpretability, limitations, and terminological context

Interpretability follows directly from the decomposition. ARIMA_PLUS can display historical and forecast trend, each seasonal component, holiday effects, spikes and dips, step changes, residuals, and covariate contributions. The paper’s Google I/O Wikipedia pageviews example shows how this componentization changes model behavior. With default US holidays, the model captures a declining long-run trend, a yearly seasonal hump around May, and event-day surges as spikes and dips. When explicit custom event holidays are supplied for Google I/O dates, those surges are reclassified as holiday effects rather than anomalies, and forecast accuracy over the three-day event window improves materially: MAE is reduced by 58%, RMSE by 62.3%, and sMAPE by 62.1%. The corresponding sum-of-squares attribution shifts from 2.3% to 50.9% for holiday effect and from 51.4% to 10.1% for spikes and dips, making the decomposition more aligned with domain knowledge (Cheng et al., 28 Oct 2025).

The framework also has clear limits. Each series is modeled independently, so it remains a local rather than global model; it assumes Gaussian residuals for intervals and anomaly detection; trend is restricted to what can be captured by segmented, difference-stationary ARIMA after preprocessing; and holiday or seasonality extraction requires sufficient history. The paper also notes that increasing AUTO_ARIMA_MAX_ORDER beyond 2 slightly worsens MASE while increasing training time, and that increasing MAX_TIME_SERIES_LENGTH from 128 to 2048 improves MASE by only about 1.86% while increasing training time by about 58%, underscoring that much of the framework’s accuracy comes from decomposition and cleaning rather than from a more complex ARIMA search alone (Cheng et al., 28 Oct 2025).

Related work suggests that the string “ARIMA+” is also used more loosely in the literature to denote ARIMA-inspired or ARIMA-augmented methods rather than the specific BigQuery ML system. Examples include a CNN-based direct multi-step ARMA or ARIMA-style module combined with RevIN for non-stationarity handling (Kim et al., 12 Sep 2025), PACF-based parameterization of ARUMA models presented as a more general algebraic form underlying ARIMA and SARIMA classes relevant to ARIMA_PLUS-type systems (Halliday et al., 2022), conditional continuation simulation of fitted SARIMA models as an ARIMA-plus workflow for scenario generation (Sak et al., 2012), and a parallel hybrid of ARIMA with polynomial classifiers combined by an analytically optimized weight (Nguyen et al., 11 May 2025). This suggests a useful distinction: in BigQuery ML, ARIMA_PLUS is the proper name of a concrete modular forecasting and anomaly-detection framework, whereas in adjacent research “ARIMA+” often functions as a broader label for extensions that retain ARIMA as a core component while adding decomposition, simulation, nonlinear modeling, or alternative parameterization.