Evaluating the statistical efficiency of robust and classical estimators for ARIMA models under contaminated time series

Introduction

Autoregressive Integrated Moving Average (ARIMA) models are among the most prevalent and efficient standard tools for time series analysis and forecasting, widely applied across various economic, financial, and social sciences. The development of this integrated methodology is credited to Box and Jenkins (1976) [14] , who established an inferential framework based on combining two primary components: the Autoregressive (AR) component, which relies on the historical values of the series, and the Moving Average (MA) component, which tracks past forecast errors, while addressing non-stationary behavior through the mechanism of differencing (𝑑). Despite the high statistical efficiency of ARIMA models under ideal standard conditions, they face severe challenges when time series are contaminated by outliers. Outliers represent a natural and inherent phenomenon in real-world data, arising from sudden structural shifts, external shocks, or measurement errors. Extreme observations cause fundamental distortions in parameter estimates and inflate error variance when relying on conventional estimators such as the Maximum Likelihood (ML) method, due to its high sensitivity and breakdown point in the presence of minor statistical contamination, ultimately leading to misleading inferences and forecasts [24] . Consequently, recent research has shifted toward developing and adopting robust estimation methods as an effective alternative. These estimators possess the capacity to downweight contaminated observations and yield stable, highly durable parameters compared to classical approaches [21] .

Research Problem

Numerous applied studies face fundamental challenges regarding the accuracy of statistical modeling and forecasting for future time series trends, negatively impacting the reliability of decisions derived from these estimates. These challenges are particularly exacerbated when dealing with time series affected by statistical contamination-specifically the presence of outliers-which causes violations of underlying classical assumptions, such as the normality of residual distributions. Consequently, the efficiency of classical estimation methods degrades significantly due to their extreme sensitivity to contaminated data, leading to biased parameter estimates and inflated standard errors. Hence, the core research problem centers on how to establish robust predictive methodologies within such contexts. This necessitates the implementation of robust estimation techniques capable of downweighting the influence of contaminated observations, thereby refining forecasting models and enhancing their performance under complex data conditions.

Research Objectives

This research aims to achieve several primary objectives, formulated as follows:  Testing and diagnosing the impact of statistical contamination and non-normality resulting from outliers on the behavior and stability of ARIMA model parameters.  Comparing the efficiency of ARIMA model estimators using the classical Maximum Likelihood (ML) estimation method against the performance of the Robust Iteratively Reweighted Least Squares (RLS) estimation method to determine the optimal model and estimator.  Evaluating the goodness-of-fit and forecasting accuracy of the estimated models across different estimators using comparative criteria, specifically the Mean Absolute Error (MAE), Akaike Information Criterion (AIC), and Bayesian Information Criterion (BIC).  Applying the estimated classical and robust models to real-world time series contaminated with outliers to validate their practical performance and empirical efficiency.

Literature Review

 Mohammed and Abd (2007) [23] investigated robust estimation methods for the mixed autoregressive moving average model ARMA(1,1) and compared them with conventional approaches through simulation experiments. Their findings revealed that the classical Least Squares (LS) method remains the optimal and most efficient choice in the absence of extreme observations. However, in the presence of outliers within the series, the robust estimation method (RAMB) demonstrated superior stability and parameter durability, highlighting the necessity of adopting robust estimators to preserve modeling accuracy.  Muler, Peña, and Yohai (2009) [24] introduced a mathematical framework developing a novel class of robust estimators for ARMA models. The proposed approach was based on a specialized formulation of M-estimators combined with an innovative residual calculation mechanism designed to constrain the influence of severe outliers, ensuring that their statistical impact remains confined solely to the period of occurrence rather than propagating throughout the series. Monte Carlo simulations demonstrated that these newly developed estimators significantly outperformed standard M-estimators and traditional diagnostic procedures, owing to their consistent asymptotic properties.  Flimmel et al. (2017) [17] addressed the econometric issues resulting from data contamination by outliers, which induce severe biases leading to inaccurate order selection and misleading ARMA models. Utilizing both simulation studies and empirical applications, the authors emphasized the critical role of robust treatments. Their results established that employing robust techniques serves as a resilient foundation for preventing parameter bias and ensuring accurate model identification and order determination in outlier-contaminated data.  Arai (2020) [8] examined the implementation of Iteratively Reweighted Least Squares (RLS) algorithms for forecasting complex time series subject to statistical gaps (missing data). The author introduced two learning and updating algorithms incorporating the Kalman Filter with RLS via time and measurement update mechanisms. The algorithmic outputs indicated that robust iterative estimators with adaptive weighting structures offer rapid convergence and smooth tracking in forecasting and fitting, underscoring the role of iterative (RLS) estimators in stabilizing dynamic forecasting models.  Zhao et al. (2022) [26] conducted a cross-sectional comparative study evaluating the predictive efficiency of classical time series models, specifically ARIMA and the Grey model GM(1,1), against Long Short-Term Memory (LSTM) deep neural networks. Predictive performance was assessed using MAE, RMSE, and MAPE criteria. The empirical findings demonstrated that advanced non-linear frameworks like LSTM exhibit superior forecasting ability in modeling complex dynamic trends compared to classical structures.  Bianco et al. (2024) [13] directly addressed outlier management in ARIMA time series models, specifically distinguishing between Innovation Outliers (IO) and Additive Outliers (AO). The authors developed robust estimation procedures as an alternative to classical Maximum Likelihood (ML) estimation, which suffers efficiency losses under data contamination. Monte Carlo simulation results confirmed that utilizing residuals derived from robust estimation enhances outlier detection capabilities and prevents the inflation of error variance compared to sensitive conventional methods.  Ardia et al. (2025) [9] focused on assessing the forecasting performance of standard models when time series experience severe shocks causing deviations from residual normality (e.g., market shocks and periods of high uncertainty). The authors implemented advanced robust estimators, such as the Huber M-estimator, comparing them with classical approaches using Mean Absolute Error (MAE) and information criteria. Their empirical results demonstrated that integrating robust methods mitigates parameter drift and achieves superior forecasting accuracy and statistical stability in contaminated series.  Otopeleanu et al. (2025) [25] introduced an advanced algorithmic solution based on a robust regularized RLS algorithm using third-order tensor decomposition. The study aimed to reduce computational complexity while maintaining parameter estimation efficiency in high-dimensional dynamic systems. Experimental results verified that the proposed robust algorithm exhibits high durability and clear superiority in downweighting distortions caused by signal contamination relative to standard RLS algorithms.

 Mano et al. (2026) [20] evaluated the impact of stochastic outlier characteristics on parameter accuracy in ARIMA models. They conducted a comparative analysis between traditional models reliant on ideal normality assumptions and robust models incorporating error-correction factors and weighting adjustments for extreme observations. The study concluded that robust estimators guarantee parameter unbiasedness and prevent misleading inferences resulting from abrupt shocks in the data.

Theoretical Framework 2.1 Time Series

A time series is defined as a sequence of observations ordered chronologically over time. These observations are typically dependent rather than statistically independent, where each observation relies on preceding ones. This temporal dependency is effectively exploited to generate reliable forecasts that reflect the sequential ordering of the data. Thus, 𝑦 𝑡 represents the observation at time 𝑡, 𝑦 𝑡-1 denotes the preceding observation, and 𝑦 𝑡+1 represents the subsequent observation [1] . Time series analysis is a crucial statistical methodology that has undergone significant developments, establishing itself as an effective tool for forecasting. This approach relies on analyzing the behavior of a given phenomenon over a specified time horizon and leveraging the extracted historical patterns to project future trends [2] .

Stationary Time Series

A time series representation over a given period  

may sometimes exhibit structural consistency with its representation over a

, indicating homogeneity in the behavior of the series. This property is referred to as stationarity. In general, a time series is visually identified as stationary through its time plot, or formally defined as stationary if it possesses a constant mean and variance over time, remaining unaffected by systematic trends or seasonal shifts. Specifically, a time series is said to be weakly stationary if it satisfies the following conditions [4] :

Where:  𝑌: The population mean of the time series.  𝑦 𝑡: The observation of the series at time 𝑡.  𝑛: The total number of observations.

Variance Stationarity

𝑉𝑎𝑟(𝑌 𝑡 ) = 𝐸(𝑌 𝑡 -𝜇) 2 = 𝜎 𝑌 2 = 𝛾 0 (3) where 𝛾 0 represents the variance of the time series, which remains constant over time and independent of 𝑡.

Absence of Trend:

There is no systematic increase or decrease in the values.

Equal relationship between values over time (Autocorrelation):

Meaning that the correlation between values depends only on the time interval between them.

Autocorrelation Function (ACF)

It is a measure that determines the degree of correlation between observations of the same variable at different time periods. It also serves as a tool to identify the characteristics of time series. The concept of autocorrelation is defined as the correlation of a random variable with its preceding or subsequent variable over a specified time period, i.e., the correlation of the series with itself or shifted by a certain period. This correlation is given by the following formula [3] :

Where:  𝜌 𝑘: Represents the autocorrelation of 𝑦 𝑡 values at lag 𝑘.  𝛾 𝑘: Represents the autocovariance at lag 𝑘.  𝜎 𝑌 2: Represents the variance of 𝑦 𝑡 values.

Partial Autocorrelation Function (PACF)

It calculates the correlation between 𝑦 𝑡 and 𝑦 𝑡-𝑘 while holding 𝑦 constant across the other periods, i.e., calculating the degree of correlation between 𝑦 𝑡 and 𝑦 𝑡-𝑘 controlling for the effect of other values. From the autoregressive model equation AR of order 𝑝, it is obtained as follows [5] :

That is, the first partial autocorrelation is equal to the first autocorrelation, and by varying the number of terms, it will yield partial autocorrelations for different values of 𝑘. Its general formula is:

Where: 𝜙 ˆ𝑘𝑘 is the estimated value of 𝜙 𝑘𝑘 .

Where: 𝜙 𝑘𝑗 is the partial correlation at period 𝑘𝑗.

Box-Jenkins Models

The Box-Jenkins methodology is a systematic approach for time series analysis and forecasting, primarily focusing on identifying and estimating forecasting models for temporal data, particularly those that can be expressed as Autoregressive Integrated Moving Average (ARIMA) models [3] . Box-Jenkins models were developed between the years (1970-1976) and are abbreviated as (B-JM). They are used in analyzing temporal data within the time domain and are divided into two categories: the first is Seasonal Models, which will not be detailed as the study relies on the second category, Non-Seasonal Models, which are used to represent two types of series [6] :

Stationary Time Series Models

These time series models are defined as possessing the property of stationarity, meaning they do not contain a general trend in their data and have a constant mean around which they fluctuate. These models include three types:

Autoregressive (AR) Models

Autoregression is defined as the correlation between the current observations of a time series and previous observations of the same series. The model order is denoted by (𝑝), where (𝑝) represents the order of autoregression as an integer. Autoregressive models of lag order AR(𝑝) can be represented as follows [7] :

For example, the first-order autoregressive model AR(1) is:

And if 𝑦 𝑡 represents the deviation of the observations from the mean 𝜇, the equation takes the following form: (10) Where the stationarity condition for the autoregressive parameter 𝜙 1 is:

The value of 𝜙 1 must be less than one in absolute value for the series to be stationary:

Where the autocorrelation coefficient 𝜌 𝑘 for this model is:

And the partial autocorrelation coefficient is equal to the parameter of the autoregressive model AR(1), while the other partial autocorrelation coefficients are equal to zero [10] . As for the second-order autoregressive model AR(2), it is:

And the autocorrelation function for this model is:

And the stationarity condition for this function is:

While the partial autocorrelation parameters (PACF) for the AR(2) model diminish, as current values depend only on values at time (𝑡 -1) and (𝑡 -2). The model will be clear at (lag 1, lag 2), and the coefficients begin to decay rapidly toward zero, meaning that the effect of subsequent lags (𝑠) is weak or insignificant.

Where 𝜙 22 is at the second lag, and each symbol can be expressed as follows:  𝑦 𝑡: Observations of the time series at time 𝑡.  𝜙 𝑖: Model parameters.

 𝑝: Order of the autoregression.  𝜇: Mean of the observations.

Moving Average (MA) Models

The studied phenomenon can be represented by another model other than the autoregressive model, which is the moving average model. The difference between this model and the previous one is that moving averages depend on error terms at different time points. Moving average models of order (𝑞)-where (𝑞) represents the order of the moving average, denoted as MA(𝑞)-can be represented by the following mathematical model [12] :

When moving average parameters are added to the autocorrelation function (ACF) equation, the equation becomes as follows:

When 𝐾 ≤ 𝑞 𝜌 𝑘 = 0 ,𝐾 > 𝑞 And the parameters of the partial autocorrelation function (PACF) never decay absolute to zero (as they decrease exponentially). For example, the first-order moving average model MA( 1) is [11] :

The behavior of the "Autocorrelation Function and Partial Autocorrelation Function" in the MA(1) model is similar to that of the AR(1) model. The second-order moving average model MA( 2) is:

Where the Autocorrelation Function (ACF) equation for the MA(2) model can be formulated as follows:

While the Partial Autocorrelation Function (PACF) equation is as follows:

The symbols above can be expressed as:  𝑦 𝑡: Observations of the series at time 𝑡.  𝜃 𝑖: Model parameters.  𝑞: Order of the moving average.

Non-Stationary Time Series Models

These models apply to represent series in which a general trend is one of their components, resulting in multiple means around which the data fluctuate. These models can be represented using the first type of models, but after removing non-stationarity from the original series by using the differencing method and taking an appropriate number of differences of the original data for time periods representing the length of the season, adding the letter (𝑑) to the notation [4] . When introducing the difference operator to the ARMA(𝑝, 𝑞) model, the model transforms into an ARIMA(𝑝, 𝑑, 𝑞) model. These models include the same stationary time series models; however, to distinguish between this type and the first type of models, the word "Integrated" is added to the model name to indicate the application of stationary time series models to non-stationary time series after converting them into stationary time series to represent this type of model. It is preferable to use the mixed model after adding the letter (I), making the model abbreviation (ARIMA), and the model becomes of order (𝑝, 𝑑, 𝑞). Converting non-stationary time series into stationary ones avoids many problems, as stationary temporal properties are independent of time. Furthermore, it becomes easier to develop stationarity theory to represent stochastic time series, meaning that the new forecasted series will also be stationary. From a practical standpoint, the calculations required to estimate the parameters of stationary models require less work. The model is defined as follows:

To make the series stationary in this model, and assuming that 𝑦 𝑡 denotes the original time series data, and 𝑤 𝑡 denotes the time series data after taking the first difference, the equation for calculating the first difference is as follows:

That is, the value of the stationary time series 𝑤 𝑡 at time 𝑡 equals the difference between the value of the original series 𝑦 𝑡 at time 𝑡 and its value 𝑦 𝑡-1 at the previous time 𝑡 -1.

That is, where 𝑤 𝑡 denotes the time series data after taking the first difference, the equation for calculating second-order differences is as follows: [16] 𝑍 𝑡 = 𝛻 2 𝑌 𝑡 = 𝛻(∇𝑌 𝑡 ) = (𝑌 𝑡 -𝑌 𝑡-1 ) -(𝑌 𝑡-1 -𝑌 𝑡-2 ) = 𝑌 𝑡 -2𝑌 𝑡-1 + 𝑌 𝑡-2 = 𝑌 𝑡 -2𝛽𝑌 𝑡 + 𝛽 2 𝑌 𝑡 = (1 -2𝛽 + 𝛽 2 )𝑌 𝑡 (22) If 𝑑 is the number of differences required to obtain a stationary time series, then the general equation for the stationary series is as follows:

Where:  𝑤 𝑡: Denotes the stationary series after taking 𝑑 differences.  𝑦 𝑡: Denotes the original series.  Δ 𝑑 𝑦 𝑡 (or (1 -𝐵) 𝑑 𝑦 𝑡 ): Represents the operation of taking 𝑑 differences of the original series.

 (1 -𝐵) 𝑑: Represents the 𝑑-th differencing using the Lag Operator (𝐵).

Accordingly, the series is expressed as 𝑦 𝑡 , and to render it stationary, 𝑑 differences are taken from it to obtain the series Z t , where 𝑑 represents the minimum number of differences required to obtain a time series that is stationary in the mean.

Stages of Applying Box-Jenkins Models in Forecasting 2.4.1 Identification & Diagnostic Phase

In this phase, the orders (𝑝, 𝑑, 𝑞) of the ARIMA model are identified, which involves monitoring the "Autocorrelation Function (ACF) and Partial Autocorrelation Function (PACF)" as a fundamental step in time series analysis [15] .  If the Autocorrelation Function (ACF) values fall within the 95% confidence interval bounds, this indicates that the autocorrelation coefficient 𝜌 𝑘 does not significantly differ from zero, demonstrating the stationarity of the series and that it is integrated of order zero (𝐼(0)). In this case, analyses can be performed using the original values of the variable 𝑦 𝑡 without requiring any transformations.

 If the Autocorrelation Function (ACF) does not fall within the 95% confidence interval bounds for several time lags, then 𝜌 𝑘 differs significantly from zero, indicating that the time series is non-stationary and requires differencing to transform it into a stationary series. When taking the first difference, 𝑑 = 1, and when taking the second difference, 𝑑 = 2. In the ARIMA(𝑝, 𝑑, 𝑞) model, 𝑝 and 𝑞 can be identified by inspecting the pattern of the Autocorrelation Function (ACF) and Partial Autocorrelation Function (PACF): 1.

When autocorrelations decay exponentially toward zero, it indicates the presence of an Autoregressive (AR) model, whose order is determined by the number of partial autocorrelations that differ significantly from zero.

2.

If partial autocorrelations decay exponentially toward zero, the model is a Moving Average (MA) model, whose order is determined by the number of statistically significant autocorrelations.

3.

If both autocorrelations and partial autocorrelations decay exponentially toward zero, this indicates an ARMA model.

The order of both AR and MA is determined by observing the behavior of the autocorrelation and partial autocorrelation functions.

If the autocorrelation function does not decay rapidly with increasing lag orders, it indicates that the time series is non-stationary and requires differencing. The diagnostic process includes reviewing these two functions to determine the stationarity of the time series and the parameters (𝑝, 𝑑, 𝑞), as shown in the table below [4] :

Estimation & Testing Phase

After identifying the model and determining its order, the parameters of this model are estimated, representing the parameters 𝜙 𝑖 in the AR(𝑝) model, 𝜃 𝑗 in the MA(𝑞) model, and (𝜙 𝑖 , 𝜃 𝑗 ) in the ARMA(𝑝, 𝑑, 𝑞) model.

First: Classical Maximum Likelihood Estimation (MLE)

Maximum Likelihood Estimation is considered the most common standard tool for estimating the parameters of time series models (Berenguer-Rico et al., 2021) [5] .This method derives its strength from the assumption that the random errors 𝑒 𝑡 are normally distributed with full normality. Classically, it is characterized by satisfying the conditions of efficiency and asymptotic consistency with large, clean sample sizes free from statistical contamination, where we write the probability density function for the AR(1) model as:

Taking the natural logarithm (ln)

and set it ∅ 𝟏 equal to zero

From which we obtain the estimator:

The probability density function of the errors for the MA(1) model, which depends on all previous errors, is written as follows:

Assuming (𝒆 𝒕 = 𝟎), the function becomes:

Taking the logarithm, we obtain:

Since 𝑒 𝑡 is a complex function in 𝜃 1 due to continuous substitution (𝑒 𝑡 = 𝑍 𝑡 -𝜃 1 𝑧 𝑡-1 + 𝜃 2 𝑧 𝑡-2 … ), the derivative 𝜕𝑒 𝑡 𝜕𝜃 1 is complex and does not yield a simple algebraic solution. Therefore, a numerical algorithm such as Newton-Raphson is used to find the value that minimizes ∑ 𝑒 𝑡 2 .

Second: Robust Recursive Least Squares (RLS) Method

The Recursive Least Squares (RLS) method represents an advanced temporal extension for model parameters step-by-step without the need to recompute calculations when new time series observations stream in [8] . Given the standard algorithm's sensitivity to breakdown when exposed to outliers, it has been robustified in recent literature by introducing time-varying regularization parameters that protect the dynamic weight matrix from inflation and computational instability (Mahadi et al., 2022). Additionally, it has been algorithmically formulated in an advanced manner based on time series matrix decomposition to reduce computational complexity during the continuous updating of the model parameters [25] .

Then we differentiate with respect to [∅ 𝟏 ] and set it equal to zero:

After simplification, the equation becomes:

And we can write the MA(1) model, where the objective function is:

Then we differentiate the function using the weighted Gauss-Newton algorithm:

Outliers and Their Statistical Detection via Boxplot

Outliers are defined in the time series literature as inconsistent observations that sharply depart from the general pattern of the series. They typically arise from measurement errors, sudden shocks, or geopolitical and security crises. Outlier contamination in time series is divided into two main types [16]:  Innovation Outliers (IO): Shocks that affect the error term at a specific moment and continue to propagate and influence future observations through the dynamic model mechanism of the time series.  Additive Outliers (AO): Exogenous incidental shocks that strike a single time observation without leaving an effect on the future memory of the series.

To detect these outliers and accurately identify contamination periods in the time series, the Boxplot is employed as a non-parametric diagnostic statistical tool based on dividing data into dynamic quartiles summarizing the geometric distribution of the series:  First Quartile (Q1): Represents the point below which 25% of the time series data falls when arranged in ascending order.  Median or Second Quartile (Q2): Represents the midpoint of the time series (50%).  Third Quartile (Q3): Represents the point below which 75% of the time series data falls.

The Interquartile Range (IQR), which measures the dispersion of the core of the series, is calculated:

Based on Tukey's rules, upper and lower detection bounds called "fences" are drawn. That is, any time observation in the series that exceeds the upper fence or falls below the lower fence is immediately classified as an extreme outlier. These outliers violate the assumptions of normality and homoscedasticity (constant variance) in time series. Given that the Maximum Likelihood method has a breakdown point approaching zero, a single outlier is sufficient to render traditional estimates of the time series model completely biased [12] . Herein lies the rationale for using robust methods, as they possess a high breakdown point reaching up to 50%, thereby mitigating the impact of these shocks [16] .

Criteria for Model Selection 1. Mean Absolute Deviation (MAD)

It is defined as the mean of the absolute differences between the actual and estimated values, serving as a measure of accuracy and quality of the estimator, given as follows: (19, 5)

Akaike Information Criterion (AIC)

Akaike (1974) proposed selecting the order 𝑝 such that the model variance remains small relative to the increase in the number of estimated parameters. Its formula is given by: [22] 𝑨𝑰𝑪 (𝒑 , 𝒒) = 𝒏 𝑳𝒏𝝈 ̂𝒂 𝟐 + 𝟐(𝒑 + 𝒒) (33)

Where:  𝑛: Number of observations.  (𝑝 + 𝑞): Total number of estimated model parameters.

 𝜎 ˆ2: Estimated error variance.

The AIC criterion is considered both a diagnostic tool and a criterion for determining the model order. The appropriate model for the time series data is selected corresponding to the minimum AIC value.

Bayesian Information Criterion (BIC) / Schwarz

Criterion: Shibata (1976) pointed out that the AIC test sometimes tends to select a higher order than the true autoregressive order. Consequently, Akaike (1979) extended the AIC criterion into the Bayesian Information Criterion (BIC) to estimate 𝑝 for the ARMA(𝑝) model according to the following formula: [3] The figure above illustrates the data distribution pattern, showing that the series of observations (blue points) is widely dispersed.

The core concentration of the data (50%) lies within the blue box body between the first and third quartiles (𝑄 1 = 141.6 and 𝑄 3 = 290.3). Regarding the skewness of the distribution, the proximity of the median line (the horizontal red line at 162.015) to the bottom of the box demonstrates that the data exhibits positive skewness (right-skewed). This is a result of the substantial upward production surges during recent years (2019-2025). Furthermore, regarding outlier identification, the boxplot reveals the presence of outliers in the series. The graph above illustrates that the original time series, represented by the upper blue plot, exhibits a long-term upward secular trend reflecting the expansion of oil facilities and increased production capacity. Quantities rose from 103.28 in 2000 to surpass 300, reaching an all-time high of 327.89 at the end of the period in 2025. This upward behavior clearly confirms that the raw series is non-stationary in the mean. To achieve stationarity due to the non-stationarity of the raw data, the first difference (𝑑 = 1) was applied to the time series, indicated by the lower green plot. Upon applying the first difference, the series eliminated the drawbacks of drift growth and began dynamically oscillating around the zero line (a constant mean), providing clear evidence of attaining stationarity. The graph above illustrates the Autocorrelation Function (ACF) and Partial Autocorrelation Function (PACF) after taking the first difference, confirming that the series has attained stationarity. We observe that the values of the Autocorrelation Function, represented by the upper blue plot (ACF), lie entirely within the critical confidence bounds (the orange dashed lines) across all time lags without exceeding them. This provides conclusive statistical evidence of time series stationarity after applying the first difference, thereby achieving full stationarity. This is further corroborated by the Partial Autocorrelation Function, represented by the lower red plot (PACF), where the function values fluctuate weakly and randomly within the safe range of statistical confidence limits (the orange dashed lines), recording no significant breaches across most lags. This joint behavior of both functions confirms standard stationarity, demonstrating that the series is well-controlled and that the residuals have been transformed into random white noise suitable for adoption in forecasting models. The table above presents the estimation results and statistical comparison for three integrated models: ARIMA(1,1,0), ARIMA(0,1,1), and ARIMA(1,1,1). Their autoregressive parameters (𝜙) and moving average parameters (𝜃) were estimated using three different methods, where the results indicate a complete superiority of the Robust Recursive Least Squares (RLS) method over the classical Maximum Likelihood (ML) method. Specifically, the Robust RLS method achieved a substantial and tangible reduction in the Mean Absolute Error (MAE) values; the error in the ARIMA(0,1,1) model dropped from 10.55432 under the ML method to its lowest level at 4.490651 using the RLS method, reflecting high parameter precision and superior capability in minimizing deviations. Thus, the Robust RLS method successfully controlled the estimates with complete econometric efficiency. The scientific rationale for this behavior lies in the fact that the RLS method relies on flexible weighting mechanisms that automatically downweight extreme observations, conferring higher econometric robustness when dealing with real-world time series contaminated by outliers. Therefore, according to the Box-Jenkins methodology, the optimal model is the one that achieves the joint minimum across all accuracy and information criteria. Upon examining the models estimated by the superior method (RLS), we observe that the ARIMA(0,1,1) model recorded the absolute lowest values among all models, achieving the lowest value for the Mean Absolute Error (MAE = 4.490651). This aligns with the other criteria, as the model recorded the lowest value for the Akaike Information Criterion (AIC = 90.38945) and the lowest value for the Schwarz-Bayesian Information Criterion (BIC = 92.8272). Furthermore, the estimated value of the moving average parameter reached (𝜃), implying that current natural gas production quantities in Iraq are not directly linearly linked to their own past values as much as they are sharply and clearly influenced by the "memory of dynamic random shocks" and accumulated errors from the previous period (𝑒 𝑡-1 ). From an economic perspective, this behavior reflects the actual reality of the country, where gas production quantities are suddenly and abruptly affected by geopolitical changes, emergency security situations, or forced technical shutdowns in oil and gas fields and facilities-all of which enter as random shocks. Consequently, the moving average component MA(1) captures their mathematical impact and dynamically corrects it with utmost precision. Accordingly, the ARIMA(0,1,1) model estimated by the Robust RLS method is adopted as the sole optimal model, serving as the basis for forecasting future years. that the time series data for natural gas production in Iraq departs from a normal distribution. This sharp deviation from normality is attributed to the contamination of the series with outliers resulting from fluctuations in production operations and the memory of geopolitical, economic, and emergency structural shocks in Iraq. This was further reflected in the geometric analyses of the Boxplot through sharp positive skewness and the presence of accelerated production surges in the second half of the study period (2019-2025).

2. The results demonstrated that examining the Autocorrelation Function (ACF) and Partial Autocorrelation Function (PACF) of the raw series exhibited clear non-stationary behavior in the mean. Specifically, the ACF values were characterized by a slow and gradual decay extending across time lags, whereas the PACF recorded a sharp and abrupt cutoff at the first lag. This necessitated treating the series using first-order differencing (𝑑 = 1), which successfully transformed the series into a fully stationary process oscillating around a constant mean with zero-variance behavior. 3. The statistical evaluation results showed complete superiority of the Robust Recursive Least Squares (RLS) method over the classical Maximum Likelihood (ML) method across all candidate models. The reason for this lies in the inability of the classical method to withstand the effect of outliers due to its extreme sensitivity to contamination (its breakdown point approaches zero). Conversely, the RLS algorithm, through its flexible weight matrix and dynamic regularization penalties, succeeded in isolating the impact of those outliers, substantially reducing the Mean Absolute Error (MAE) with an improvement rate exceeding expectation in some models. 4. The results of accuracy criteria and comparison measures showed that the optimal model for representing the series is the robust Integrated Autoregressive Moving Average model, ARIMA(0,1,1), estimated via the RLS method. It achieved the absolute lowest values for error and information criteria. 5. The results revealed that the moving average parameter of the adopted model recorded an estimated value of 𝜃, alongside a complete absence of autoregressive parameters (𝜙 = 0). Economically and practically, this result proves that natural gas production in Iraq is not governed by a strong, continuous linear link with its historical production values. Instead, it is entirely and impulsively influenced by the memory of instantaneous random shocks, technical shutdowns, and sudden security and political fluctuations affecting oil facilities-where the moving average component MA(1) captures and dynamically corrects this impact.

Recommendations

1. The researcher recommends the adoption of the robust model derived in this study, ARIMA(0,1,1) estimated by the RLS method, by the Iraqi Ministry of Oil and the Ministry of Planning as a primary benchmark tool for forecasting natural gas production quantities for upcoming years. This will contribute to formulating realistic development plans built on scientific readings free from outlier distortion. 2. We recommend that service and production departments and institutions avoid using traditional estimation methods (such as Maximum Likelihood, ML, or Ordinary Least Squares, OLS) when dealing with real-world time series affected by sharp exogenous shocks and emergency conditions. Instead, they should utilize robust algorithms (such as RLS), which have proven their capability in reducing bias and enhancing modeling efficiency. 3. We recommend developing hybrid models that combine robust time series models (RLS-ARIMA) with artificial intelligence techniques, neural networks, or grey systems to evaluate forecasting efficiency under small sample sizes and high-uncertainty conditions.

References

1. Ahmed D, Atta GM. Forecasting consumer price indices in Iraq. Proceedings of the 3rd Scientific Conference, University of Kerbala. 2018;2.
2. Jaber AG, Ankoush AD. Using Box-Jenkins models to forecast the volume of outgoing and incoming postal parcels. Proceedings of the 3rd Scientific Conference, University of Kerbala. 2018;2.
3. Shamal AH, Hamid YK. A study on the number of divorce cases in the city of Hilla using time series in neural networks. Proceedings of the 3rd Scientific Conference, University of Kerbala. 2018;2.
4. Abd Al-Kadhim IJ. Comparing the efficiency of black-box models with grey system models for time series forecasting with an application. Kut Journal of Economics and Administrative Sciences. 2025;17(59):26-53.
5. Al-Ujeili SSM, Ali FA. Forecasting inflation rates in Iraq using integrated Box-Jenkins models for the period 1971-2015. Proceedings of the 3rd Scientific Conference, University of Kerbala. 2018;2.
6. Ministry of Oil -Iraq. Quantity of natural gas produced in Iraq for the period (2000-2025). Baghdad, Iraq; 2026.Ahmed BK, Rahim SA, Maaroof BB, Taher HA. Comparison between ARIMA and Fourier ARIMA model to forecast the demand of electricity in Sulaimani Governorate. Qalaai Zanist Journal. 2020;5(3):908-940.
7. Al-Nasser AMH, Al-Saffawi SY. Comparison between ordinary and robust estimators for low-order bivariate combined time series models. Journal of Applied Statistical Science. 2005;14(3):215-228.
8. Arai K. Recursive least square: RLS method-based time series data prediction for many missing data. International Journal of Advanced Computer Science and Applications. 2020;11(11):1-7.
9. Ardia D, Hoogerheide L, van Dijk HK. Macroeconomic forecasting with outlier-robust VAR and ARIMA models. Journal of Applied Econometrics. 2025;40(3):312-329. DOI:10.1002/jae.3088.
10. Askari M, Askari H. Time series grey system prediction-based models: gold price forecasting. Trends in Applied Sciences Research. 2011;6(11):1287.
11. Barreto H, Maharry D. Least median of squares and regression through the origin. Computational Statistics & Data Analysis. 2006;50(6):1391-1397.
12. Berenguer-Rico V, Johansen S, Nielsen B. A model where the least trimmed squares estimator is maximum likelihood. Technical Report. Department of Economics, University of Oxford; 2021.
13. Bianco AM, Boente G, Rodrigues IM. Outlier detection in regression models with ARIMA errors using robust estimates. Journal of Time Series Analysis. 2024;45(2):184-201. DOI:10.1111/jtsa.12718.
14. Box GEP, Jenkins GM. Time series analysis: forecasting and control. Holden-Day; 1976.
15. Chatfield C. The analysis of time series: an introduction. 6th ed. Chapman and Hall/CRC Press; 2003.
16. Donatos GS, Meintanis SG. Robust estimators of AR-models: a comparison. Journal of Statistical Computation and Simulation. 2003;73(3):157-167.
17. Flimmel S, Čamaj M, Malá I, Procházka J. Robust methods for ARMA model order identification in the presence of outliers. Journal of Financial and Actuarial Metrics. 2017;12(1):45-59.
18. Kwon S, Oh M, Park B. ARIMA model selection and prediction intervals under contaminated distributions. Journal of Forecasting. 2026;45(1):115-130. DOI:10.1002/for.3195.
19. Mahadi M, Ballal T, Moinuddin M, Al-Saggaf UM. A recursive least-squares with a time-varying regularization parameter. Applied Sciences. 2022;12(4).
20. Mano R, Wijaya A, Setyawan I. Evaluating the validity of time series models under outliers characteristics. Communications in Statistics -Simulation and Computation. 2026;55(4):402-418. DOI:10.1080/03610918.2024.2311504.
21. Maronna RA, Martin RD, Yohai VJ. Robust statistics: theory and methods (with R). 2nd ed. Wiley; 2019.
22. Mills TC. Applied time series analysis: a practical guide to modeling and forecasting. Academic Press; 2019.
23. Mohammed SB, Abd NA. Robust estimators of ARMA(1,1) model (empirical study). Al-Rafidain Journal of Computer Sciences and Mathematics. 2007;4(2):89-104.
24. Muler N, Peña D, Yohai VJ. Robust estimation for ARMA models. The Annals of Statistics. 2009;37(2):816-840. DOI:10.1214/07-AOS570.
25. Otopeleanu RA, Paleologu C, Benesty J, Stanciu CL, Dogariu LM, Costea RL. Robust regularized recursive least-squares algorithm based on third-order tensor decomposition. Algorithms. 2025;18(3):768.
26. Zhao D, Zhang H, Cao Q, Wang Z, He S, Zhou M, et al. The research of ARIMA, GM (1,1), and LSTM models for prediction of TB cases in China. PLoS ONE. 2022;17(2):e0262734. DOI:10.1371/journal.pone.0262734.