Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
68,240 characters · 11 sections · 24 citation commands
Model selection confidence sets for time series models with applications to electricity load data
\doublespacing
{\it Keywords:} model uncertainty, model selection confidence set, lower boundary models, electricity load forecasting, time series.
Model selection is an important tool for practitioners and policymakers as it enhances the interpretability and practical relevance of statistical models. This is particularly evident in time-series analysis, where competing models can deliver similar in-sample fits while producing substantially different forecasts stock,Rossi. More often than not, however, the data do not clearly support a unique model, so model selection uncertainty becomes a central concern.
This issue is important in commodity markets, particularly for partially storable commodities such as electricity, where country-level load modeling has gained prominence. Electricity load exhibits pronounced intra-day variation and responds to a broad set of drivers, including weather conditions station,recency,lowvoltage,ireland, economic indicators FEZZI,LiddleHasanov2024Correction, and seasonal patterns taylor2,DORDONNAT2008566,seasonal,bookseason, making the identification of a single best specification difficult hong2019gefcom2017,dong2025stelfsurvey. Moreover, the optimal model structure may change over time as demand patterns evolve (concept drift) and regulatory reforms reshape consumption behavior Enrich2024TOUSpain,azeem2024tla,cao2025opdf. Understanding and controlling this variability matters because even modest forecasting improvements can reduce operating costs by improving unit commitment and dispatch, informing reserve requirements under net-load uncertainty, and lowering congestion-management needs pierro2020residualload,wang2023valueprob,bernecker2025congestion,hasan2025lfreview. This becomes even more prominent during extreme load conditions (e.g., heat- or cold-driven demand spikes), which are important to predict, and model comparisons can depend on how performance is evaluated when attention is placed on extremes lerch2017forecasters.
Our work is motivated by a specific empirical problem: modeling hourly electricity load in Italy over the period from 2019 to 2023. This type of data, augmented with temperature and renewable generation measures, naturally gives rise to many competing specifications, making model selection uncertainty unavoidable. Italy provides a particularly informative case, as the sample includes major structural disruptions, particularly the COVID-19 pandemic and the 2022 energy crisis, that substantially affect load dynamics and complicate model choice. As an illustration of selection uncertainty, Figure (ref) shows the logarithmic root mean squared error (RMSE) of 28,672 linear time series model specifications fitted with electricity load data described in Section (ref). The $-\log_{10}(p\text{-value})$ from likelihood ratio (LR) tests comparing a large benchmark model to all candidate submodels are also reported. Notably, there are a number of models with similar prediction performance that cannot be discarded by the LR test. All these models provide statistically indistinguishable explanations of the true data generating process at a given level of accuracy.
Existing approaches to model selection typically focus on selecting a single optimal model. They can be broadly categorized into three classes: hypothesis testing diebold2002,Romano2005,Giacomini,CLARK2007291, information criteria Billah01102005,EBIC,AIC,ding2018model, and penalized methods tibshirani1996regression,fan2001variable,zou2005regularization,penalized,MCP,UematsuTanaka2019, Kock2020. Information criteria such as the Akaike Information Criterion (AIC) and the Bayesian Information Criterion (BIC) often lead to substantially different model choices, particularly in finite samples hansen2007least, while regularization methods such as LASSO exhibit considerable selection variability in noisy or high-dimensional settings BuhlmannVanDeGeer2011HighDim,Kock2020. When multiple models provide comparable explanatory power, conditioning inference on a single selected specification may overlook meaningful structural differences and understate the uncertainty inherent in the model selection process.
In this paper we address model selection uncertainty by introducing the Model Selection Confidence Set (MSCS) for univariate time series models involving autoregressive and moving average components with exogenous inputs (ARMAX). The MSCS is defined as the set of all model specifications that are statistically indistinguishable from the true data-generating process at confidence level $(1-\alpha)\%$. Analogous to confidence intervals in parameter estimation, the MSCS size serves as a diagnostic for selection uncertainty: large sets reflect higher noise, with many models equally supported by the data while concentration around fewer models indicates lower uncertainty on the underlying data generating process. A key feature of the MSCS is its subset of Lower Boundary Models (LBMs), i.e., the minimal models in the MSCS that cannot be further simplified without significant information loss. Studying LBMs helps identify key model terms: predictors and lags that appear across many LBMs are likely essential, whereas those that occur only sporadically may be less influential. In this way, the MSCS provides deeper insight than traditional single-model selection or standard model-averaging techniques.
Several studies have been devoted in developing sets of models, each statistically justified by the data, to better account for model uncertainty. MCS propose constructing a model confidence set (MCS) from a given set of candidate models by a sequence of equivalence tests on the currently remaining models, followed by a rule that removes the worst performers. Their method obtains a subset of the original models that is meant to contain the set of models with the best empirical performance under some given loss function. While sharing a similar motivation, our methodology builds on the confidence set framework of VSCS for variable selection in linear models, which employs screening approaches based on the F-test, and later extended by Ferrari2019 and MCBounds to more general settings using the likelihood ratio test. The crucial distinction lies in their objective: the MCS seeks to identify the best predictive models from the candidate pool, whereas the latter framework aims to identify the true underlying process that generated the data.
This paper makes three main contributions. First, we develop a methodological framework for constructing confidence sets of univariate time series models, extending variable selection confidence set approaches to accommodate temporal dependencies in time series data. Second, we evaluate finite-sample performance through Monte Carlo experiments, confirming that the MSCS attains the desired asymptotic coverage probability and documenting how the sizes of the MSCS and the LBM vary with sample size and across different data-generating processes. Third, we apply the methodology to Italian hourly electricity load data, showing that accounting for model selection uncertainty substantially enhances forecasting robustness and reveals multiple plausible predictor combinations across different hours of the day.
The remainder of the paper is as follows. Section (ref) introduces the ARMAX modeling framework and develops the Model Selection Confidence Set (MSCS) methodology in this context. In the same section, we formally define Lower Boundary Models within the ARMAX framework. Section (ref) investigates the finite-sample behaviour of the procedure through Monte Carlo experiments under a range of data-generating processes and model misspecification scenarios. The analysis focuses on coverage probabilities, the size and composition of the MSCS, and the ability of the method to recover structurally relevant specifications. Section (ref) applies the proposed methodology to hourly Italian electricity load data, examining both the empirical composition of the MSCS and the predictive performance of the models contained in the set. Particular attention is given to forecast accuracy and the stability of model terms inclusion across competing specifications. Section (ref) concludes with a discussion of the main findings and directions for future research.
Let $y_t$ be the electricity load at time $t= 1,...,T$. The autoregressive moving average model of order $(p,q)$ that incorporates $s$ exogenous variables observed with $r$ lags, denoted as ARMAX$(p,q,r)$, is defined as:
where $x_t=(x_{t1}, \dots, x_{ts})^\top $, is a $s\times 1$ vector of predictors, $\varepsilon_1, \dots ,\varepsilon_t$ are independently and identically distributed errors with zero mean and variance $\sigma^2$. Let $\boldsymbol{\theta }= (\alpha, \boldsymbol\beta^\top, \boldsymbol\vartheta^\top , \boldsymbol\eta^\top)^\top$ be the overall parameter vector where $\boldsymbol{\beta}= (\beta_1, \dots, \beta_p)^\top \in \mathbb{R}^p$, $\boldsymbol{\vartheta} = (\vartheta_1, \dots, \vartheta_q)^\top \in \mathbb{R}^q$ and $\boldsymbol{\eta}= (\eta^\top_1, \dots, \eta^\top_r)^\top \in \mathbb{R}^{sr}$ are the coefficients associated with lagged response variables, lagged error terms, and exogenous variables, respectively. We consider the situation where multiple elements of $\boldsymbol{\theta}$ might be equal to zero. A model \( m \) is defined as a tuple \( m = \{S_{\text{AR}}, S_{\text{MA}}, S_{\text{X}}\} \), where \( S_{\text{AR}} = \{ j \in \{1, \dots, p\} \mid \beta_j \neq 0 \} \) represents the set of autoregressive lags with nonzero coefficients in the parameter vector $\boldsymbol{\beta}$, \( S_{\text{MA}} = \{ k \in \{1, \dots, q\} \mid \vartheta_k \neq 0 \} \) is the set of moving average lags with nonzero coefficients in $\boldsymbol{\vartheta }$, and \( S_{\text{X}} = \{ (i, l) \mid i \in \{1, \dots, s\}, l \in \{1, \dots, r\}, \eta_{i,l} \neq 0 \} \) corresponds to the set of exogenous variables and their lags that have nonzero coefficients. The dimension of model $m$ is denoted by $d_m = |m|$; while the space of all feasible models is denoted by $\mathcal{M}$ and has cardinality \( |\mathcal{M}| = 2^{p + q + rs } \) , when all lags are considered. The true model is denoted by \( m^* = (S_{\text{AR}}^*, S_{\text{MA}}^*, S_{\text{X}}^*) \), where each subset contains the indices of terms with nonzero coefficients in the true data-generating process.
For a given confidence level $(1-\alpha)$, where $0 < \alpha < 1$ (e.g., 0.95 or 0.99), we aim to construct a confidence set of ARMAX models, $\widehat{\mathcal{M}}_\alpha$, satisfying
This property, referred to as asymptotic coverage, ensures that the true data generating process $m^\ast$ is included in the MSCS with large probability in large samples. In contrast to classical model selection procedures that return a single “best” specification, the MSCS acknowledges the intrinsic uncertainty of the selection step and delivers a set of models that are statistically indistinguishable at the prescribed confidence level. As a result, inference and prediction are no longer conditioned on a single selected model but instead accounts for the sampling variability inherent in model comparison. The coverage property therefore provides a formal frequentist guarantee that model uncertainty is explicitly controlled, rather than implicitly ignored.
In the context of ARMAX models, we rely on the likelihood ratio (LR) test to decide whether a given candidate model $m$ is statistically plausible compared to a full model, $m_f$, which includes all potentially useful predictors and lag terms. The full model serves as a comprehensive benchmark, capturing the maximum possible structure in the data. In order for the LR test to be meaningful, we assume that $T> d_f = |m_f|$; then the candidate model $m$ is nested in $m_f$ and both models ($m, m_f$) are uniquely estimated via maximum likelihood under standard regularity conditions. Model $m$ is rejected when
where $L(m)$ and $L(m_f)$ are the maximized likelihood functions for models $m$ and $m_f$, respectively, and $q_{1-\alpha}(\delta)$ is the $(1-\alpha)$-quantile for the central chi-square random variable with $\delta$ degrees of freedom, denoted by $\chi^2_{\delta}$. The MSCS $\widehat{\mathcal{M}}_\alpha$ is defined by all the models not rejected by the LR test, i.e.,
and by default, the full model is included in $\widehat{\mathcal{M}}_\alpha$. By construction, if the true model is not the full model, the confidence set $\widehat{\mathcal{M}}_\alpha$ has asymptotic coverage probability of \( 1-\alpha \) as $T \rightarrow \infty$, as stated in Eq. ((ref)).
This result follows trivially from the fact that when \( m = m^{\ast} \), under the null hypothesis that the restricted model \( m \) is correctly specified, we have that \( \Lambda(m) \xrightarrow{d} \chi^2_{d_f - d_m} \), as \( T \to \infty \). A proof for general likelihoods and diverging number of parameters is available in Theorem 1 of Ferrari2019. Although we do not provide formal results in this paper, following similar steps one can show that the same behavior holds in the current setting for fixed number of parameters, $d < T$, under standard regularity conditions for the LR test in ARMAX models involving: (i) the correct specification of the true model, (ii) identifiability of parameters, (iii) finite and nonsingular Fisher information matrix, (iv) weak dependence conditions such as strong mixing or ergodicity, and (v) regularity of the innovation process, where the error terms \( \varepsilon_t \) are assumed to be i.i.d. or weakly dependent with finite variance. Our procedure works also when alternative tests, such as the Wald or the Lagrange multiplier (score) tests (e.g., see WaldLR), are employed for model screening since they are asymptotically equivalent under standard regularity conditions.
The confidence set \( \widehat{\mathcal{M}}_\alpha \) allows us to determine whether a model \( m \) is overly parsimonious, meaning that it may exclude important variables. If \( m \notin \widehat{\mathcal{M}}_\alpha \), such a model likely overlooks some key terms. Unlike traditional model selection methods that force a single best model, the MSCS explicitly quantifies selection uncertainty by identifying a set of plausible models. This prevents overconfidence in a single specification and mitigates the risk of misleading conclusions when multiple models fit the data similarly well. Moreover, the structure of \( \widehat{\mathcal{M}}_\alpha \) provides insights into predictor stability, distinguishing between core explanatory variables and those whose inclusion is selection-dependent. The size of the MSCS is also informative: a small confidence set suggests strong support for a particular specification, while a larger one signals greater uncertainty, revealing multiple competing models.
The MSCS described in the previous section could have a large size mainly due to the possibility of including irrelevant parameters that survive the likelihood ratio screening. To manage the potential size of the MSCS, we focus on a smaller, yet informative subset, the set of lower boundary models (LBMs), denoted by $\widehat{\mathcal{B}}_\alpha$, defined by the set of models in $\widehat{\mathcal{M}}_\alpha$ without any nested sub-models in $\widehat{\mathcal{M}}_\alpha$. Intuitively, the LBMs is the collection of the most parsimonious models that remain well-supported by the data at a given confidence level.
VSCS study the LBMs in the context of linear models finding that both their composition and size of $\widehat{\mathcal{M}}_\alpha$ provide valuable insights into overall variable selection uncertainty. Their findings on LBMs can be extended to the ARMAX setting when the total number of parameters in the full model, $d$, is fixed, and the sample size, $T$, is large under the following scenarios.
As \( T \to \infty \), the lower boundary models increasingly include the terms of the true data-generating process with increasing probability, effectively identifying all relevant model terms. Particularly, as \( T \to \infty \), if \( m^\ast \) is not the full model, then $P(\widehat{\mathcal{B}}_\alpha = \{m^\ast\}) \to 1 - \alpha$, see VSCS for this result in linear models. A similar result may be established for ARMAX processes by following Theorem 2 of Ferrari2019, provided that standard regularity conditions hold to ensure the asymptotic validity of the likelihood ratio test.
In this paper, we observe that the cardinality and the composition of LBMs provide valuable insights into model terms selection uncertainty in time-dependent settings. Furthermore, by examining how these lower boundary models relate to more complex models in the MSCS, researchers can gain insights into the essential versus potentially superfluous terms, thus offering a principled approach to understand the trade-off between model simplicity and explanatory power.
The finite-sample behavior of $\widehat{\mathcal{B}}_\alpha$ is evaluated in Section (ref) via Monte Carlo experiments, while in Section (ref) we illustrate how the cardinality and composition of LBMs inform model selection uncertainty in electricity load forecasting.
Moreover, to reduce the risk of being excessively parsimonious by omitting important terms that might appear in only some models we explore the utility of a combination of LBMs, called union model. This specification is defined as $m^{\text{max}} = \cup_{m \in \widehat{\mathcal{B}_\alpha}} m$. The union model is especially valuable when dealing with noisy temporal data: while individual models might capture specific patterns, the union model $m^{\text{max}}$ ensures all relevant temporal dynamics are represented. This is confirmed in Monte Carlo experiments in Section (ref) showing that $m^{\text{max}}$ consistently contains the true model across all settings, with a probability of nearly 1, underscoring its robustness in covering the correct model, even in small samples.
To quantify the relevance of specific terms, we define the inclusion importance metrics, denoted as \( {II}(\theta_j) \), defined as the frequency with which a given parameter $ \theta_j \in \Theta$ appears as nonzero in the LBMs of the selected ARMAX confidence set. This metric provides insight into the relevance and stability of individual predictors across competing models and it is defined as
where $K(\theta_j)$ is the number of times $\theta_j$ appears in $\widehat{\mathcal{B}}_\alpha $ and \( | \widehat{\mathcal{B}}_\alpha | \), represents the cardinality of the lower boundary $\widehat{\mathcal{B}}_\alpha $. In the context of ARMAX modeling, the inclusion importance determines which exogenous variables and lagged terms are most important when sampling variability is considered. A large value of \( {II}(\theta_j) \) indicates that $\theta_j$ is a core model term, whereas a small value suggests that $\theta_j$ has limited importance.
We extend the inclusion importance measure to the entire confidence set, \( \widehat{\mathcal{M}}_\alpha \). Specifically, we define the inclusion importance for the entire MSCS as $ \widetilde{II}(\theta_j) = {\widetilde{K}(\theta_j)}/{|\widehat{\mathcal{M}}_\alpha|}$, where \( \widetilde{K}(\theta_j) \) denotes the number of times the parameter \( \theta_j \) is nonzero across all models in \( \widehat{\mathcal{M}}_\alpha \). Since the LR test tends to retain models that are larger than the true model, irrelevant predictors appear randomly in approximately half of the models within \( \widehat{\mathcal{M}}_\alpha \). To account for this behavior, we define a normalized inclusion importance metric that corrects for the expected inclusion rate of spurious predictors:
The normalized inclusion importance emphasizes the terms with strong inclusion signals while down-weights those that appear inconsistently due to overfitting. A large value of $\widetilde{II}_s(\theta_j)$ indicates that a predictor is included in all models, while a value near zero suggests weak or spurious inclusion.
To study the finite-sample properties of the MSCS and LBMs, we conduct several Monte Carlo experiments by generating data from six ARMAX models specified by various configurations of the autoregressive, moving average and exogenous components as reported in Table (ref). For all the models, we set $\sigma =0.8$ and consider both independent and correlated predictors generated from a multivariate normal distribution taking, uniform correlation between each predictor pair with values $\rho \in \{0, 0.7\}$. All the considered models do not include the intercept, specify a zero-lag for the exogenous variables and have normally distributed error terms. In each experiment, we consider 500 Monte Carlo samples of size $T \in \{ 250, 1000, 10000\}$. The MSCS is computed using exhaustive enumeration of the sample space consisting of all feasible submodels of the full model ARMAX($3$, $2$, $6$) and consider nominal confidence levels $(1-\alpha) \in \{0.95,0.99\}$.
For each setting, we study the key properties of the MSCS and related LBMs in finite samples. Particularly, we compute Monte Carlo estimates of the following quantities: 1) MSCS size $|\widehat{\mathcal{M}}_{\alpha}|$; 2) average size of models in the lower boundary $\widehat{\mathcal{B}}_{\alpha}$; 3) coverage probability $P_\alpha^* = P(m^* \in \widehat{\mathcal{M}}_{\alpha})$; 4) expected distance $ d^{*, max} = E[d(m^\ast, m^{\text{max}})]$, where $d(m_1, m_2)$ is the Hamming distance counting the number of different elements in $m_1$ compared to model $m_2$; 5) probability that the true model is included in the LBM union model, $P_\alpha^{max} = P(m^\ast \subseteq m^{\text{max}})$, where $m^{\text{max}}$ is defined in Section (ref).
Table (ref) shows Monte Carlo estimates for the MSCS and LBM set size. The number of LBMs is relatively small, suggesting that even when the MSCS contains many models, only a few of them are essential for capturing the underlying process. The presence of correlated exogenous predictors (\(\rho = 0.7\)) does not significantly alter these trends, indicating that the selection procedure remains robust under moderate correlation.
When comparing different confidence levels, both LBM set and MSCS size increase with the confidence level, given the stricter threshold for model rejection.
The behavior of the MSCS and LBM set exhibits distinct dependencies on sample size. The LBM set size consistently decreases as the sample size increases, converging to 1, as expected, but the rate of convergence varies depending on model complexity. For the ARMA(1,0) model (A--B), the LBM size remains close to 1 even in small samples, indicating precise model selection. In contrast, the ARMA(1,1) model (C--D) initially shows larger LBM sizes of about $2$ due to greater uncertainty in small samples caused by estimation fluctuations but gradually converges to 1 as the sample size increases. The ARMA(2,1) case (E--F) exhibits the slowest convergence, with an initial over-selection of parameters ranging from about $2.1$ to $2.2$, as the model selection process struggles with the interaction of autoregressive and moving average components. The presence of MA terms introduces greater model uncertainty, delaying MSCS stabilization. However, as the sample size grows, estimation variance decreases, leading to more accurate model identification, and the LBM size approaches 1.
The behavior of the full MSCS is less consistent and depends on whether the true model is ARMA(1,1) or ARMA(1,0). Since ARMA(1,0) is structurally farther from the full ARMA(3,2) model, more models survive the LR test, making rejection less frequent. In small samples, the LR test struggles to distinguish between an AR(1) model with short-term noise and an ARMA(1,1) model, leading to frequent misclassification. With increasing sample size, reduced estimation variance allows for a clearer identification of the MA(1) component, shrinking the MSCS. Conversely, when the true model is ARMA(1,0), small samples make it difficult to detect spurious MA components, resulting in a more conservative selection. As sample size grows, small but systematic autocorrelations in the residuals may mimic a MA structure, increasing the likelihood of selecting models with unnecessary MA terms, thereby expanding the MSCS.
Table (ref) shows estimates of the coverage probability of the MSCS and the distance of the true model from the LBMs. While the coverage approaches the nominal level with large $T$, in small samples it is not the case due to the finite-sample behavior of the likelihood ratio (LR) test, which tends to over-reject simpler models when the true model is nested within a larger one. Thus, in relatively small samples the LR screening somewhat struggles to distinguish between alternative specifications due to increased estimation variability and bias. The Hamming distance between $m^{\text{max}}$ and $m^\ast$ displays high values in small samples, thus indicating pronounced model selection uncertainty, but decreases sharply with increasing sample size. Notably, the union of LBMs consistently contains the true model across all settings, with a probability of nearly 1, underscoring the robustness of the LBM approach in covering the correct model, even in small samples. These results confirm that while small-sample variability introduces uncertainty in model selection, the union of LBMs provides a safe model that ensures the presence of the true term, and increasing the sample size leads to more accurate model identification.
Recent empirical work on electricity load forecasting lacks consensus on a comprehensive set of load drivers. Temperature and calendar effects are the most acknowledged determinants bashiri2023, with load responding nonlinearly to temperature henley1997non, damm2017effects. Given this uncertainty, we adopt the MSCS methodology to assess the model uncertainty and study the importance of individual model terms.
We analyze hourly electricity load data (MWh) for Italy obtained from the London Stock Exchange Group (LSEG) (\url{https://www.lseg.com/en/data-analytics}), from January 1 2019 to September 29 2023 consisting of a total of 41592 observations. Figure (ref) shows the hourly load (MWh) time series for Italy, which exhibits moderate annual seasonality, with peaks during winter and summer due to increased heating and cooling loads. We collect a set of key predictors to explain variations in electricity load; these include forecasted solar (Solar) and wind (Wind) generation in MWh as renewable energy sources, as well as temperature (Temperature) and its squared term (Temperature$^2$) to capture nonlinear effects on load \footnote{Hourly temperature data were retrieved from the Milan Lambrate meteorological station through the Lombardy Regional Agency for Environmental Protection (\href{https://www.arpalombardia.it/temi-ambientali/meteo-e-clima/form-richiesta-dati/}{ARPA Lombardia}), specifically form sensor ID 2001.}. Calendar effects are modeled through indicator variables for months from December to February (Winter) and from June to August (Summer), for weekends and public holidays (\texttt{Week/Hol.}), and for the first COVID-19 lockdown in Italy (\texttt{Lockdown}), which had a noticeable impact on electricity consumption patterns. To account for the inherent persistence in load, we consider intra-day lagged values of electricity load of one (\texttt{Lag$_{1}$}), two (\texttt{Lag$_{2}$}) and three (\texttt{Lag$_{3}$}) hours.
To construct the MSCS, we compare the full ARMAX($7$,$2$) model incorporating 11 exogenous predictors detailed in Section (ref) against all possible submodels \footnote{To work on a manageable model space we do not allow for gaps in auto-regressive and moving-average orders and we do not consider lagged values of external predictors. Cardinality of $\mathcal{M}$ is then restricted to 28672 different ARMAX structures.} based on the LR statistic as defined in Equation ((ref)). Setting $p=7$ ensures coverage of weekly consumption cycles, capturing recurring electricity load at the same hour across different days. Similarly, $q=2$ provides sufficient coverage of short-term autocorrelations in forecasting errors, particularly those influenced by daily fluctuations. Given the hourly nature of the data, we split the dataset into 24 separate time series, each corresponding to a specific hour of the day, and focus mainly on working hours (8-17).
In Table (ref), we present the size of the MSCS and its corresponding LBMs set, revealing distinct patterns in model selection uncertainty across working hours. The MSCS size remains consistently low throughout the morning and early afternoon, suggesting moderate uncertainty. This stability ends at hour 15, when the MSCS size increases by several orders of magnitude, indicating a rise in uncertainty, and then vanishing the hour after.
The LBM set size remains small across all hours, yet exhibits meaningful variation. In the morning, the set size is minimal, suggesting only moderate uncertainty in identifying the simplest adequate model. Moving into the afternoon, the set size increases, corresponding to increased uncertainty about which minimal specification best captures load dynamics. Crucially, this expansion in LBM diversity peaks at hour 14, one hour before the explosion in MSCS size. This temporal sequence, combined with a substantial reduction in the model's autoregressive order $p$ around the same time, provides evidence of a significant structural shift in the underlying uncertainty pertaining the data-generating process during the afternoon hours.
Table (ref) presents the inclusion importance of different model terms of the lower boundary models corresponding to the MSCS for working hours. Despite the varying uncertainty reflected in LBM set sizes, the core explanatory variables such as Week/Hol., Temperature, Temperature$^2$ and Solar and key lags like Lag$_{1}$, Lag$_{2}$ remain largely consistent across all hours. Consistent selection of such predictors confirms the dominant influence of intraday hourly lags, calendar effects, and weather conditions on electricity load.
Further insights emerge when comparing the structures of LBMs selected at the 99% and 95% confidence levels. For instance, the 99% confidence level produces greater diversity in the model's parameters specification, admitting various autoregressive structures across the boundary models. Conversely, for the same hours, the 95% confidence level often results in a consensus on a single, specific structure, indicating a clearer model preference.
Beyond these core terms, the inclusion of environmental and event-based predictors reveals more nuanced patterns in electricity load. Solar appears in at least one boundary model for all daylight hours, aligning with its direct impact on the electricity grid through photovoltaic generation and its correlation with cooling load. Wind, however, is included far less frequently, suggesting its influence is more sporadic or less significant than other factors. The seasonal and event variables---Summer, Winter and Lockdown---also appear in specific LBMs rather than forming a consistent pattern, which implies their impact is concentrated in certain hours of the day rather than exerting uniform effects.
The autoregressive and moving average components exhibit a structured but highly flexible pattern. While many hours favor a weekly autoregressive dependency with $p=7$, reflecting the strong weekly cyclicality of electricity consumption, the models show significant adaptation throughout the day. Shorter autoregressive orders ($p=1,2,4$) are selected at various hours, indicating that while weekly patterns are important, shorter-term dynamics become more relevant at specific times. For the moving average component, both $q=1$ and $q=2$ are frequently included reflecting flexibility in capturing short-term forecast errors.
The results highlighted in Table (ref) are confirmed by the findings in Figure (ref), where we show the normalized inclusion importance introduced in Equation ((ref)) for key predictor variables at different hours. The contribution of each variable to explaining load variation is clearly time-dependent. Moreover, the close agreement between inclusion importance at the confidence set level ($\widetilde{II}_s$) and at the LBM level ($II$) indicates that the main findings on intra-day hourly loads, calendar effects, and temperature are driven by features common to the whole confidence set, rather than by the choice of any particular parsimonious specification.
Next, we use the MSCS to assess whether models selected by common information criteria are trustworthy. Table (ref) compares predictor inclusion patterns between our method and the single best model chosen by three conventional criteria: the Akaike Information Criterion (AIC), the Bayesian Information Criterion (BIC), and the model minimizing out-of-sample Root Mean Squared Error (RMSE). We examine hours 8, 12, and 16 to capture morning, midday, and evening load dynamics. The analysis reveals substantial temporal variation in whether criterion-selected models fall within the MSCS. At hours 8 and 12, all three criteria—AIC, BIC, and RMSE—select models that are included within the MSCS at the 99% confidence level, indicating that these models are statistically plausible during morning and midday periods. However, at hour 16, both BIC and RMSE identify models outside the MSCS bounds, while only AIC selects a model within our confidence set. This result reflects the fundamental trade-offs embedded in these selection criteria. AIC's weaker penalty for model complexity tends to favor models with more parameters, which increases the likelihood that AIC-selected models fall within the MSCS when the true data-generating process exhibits substantial structure. Conversely, BIC's stronger parsimony penalty favors more parsimonious specifications, which may fall outside the MSCS when the underlying load dynamics require additional complexity to achieve statistical plausibility. The divergence at hour 16 suggests that evening load patterns exhibit greater complexity than earlier periods, requiring a richer model structure that BIC's penalty excludes but that the MSCS deems statistically necessary. On the other hand, RMSE's focus on minimizing prediction error, rather than balancing fit and complexity through an information-theoretic penalty, can similarly lead to underfitted models that lack the statistical support required for MSCS inclusion. These findings adds to the value of the MSCS in providing a rigorous benchmark for evaluating whether conventional criteria strike an appropriate balance between parsimony and model adequacy. Furthermore, examining individual predictors reveal nuanced temporal patterns. Among intraday lags, the one-hour lag (Lag_1) appears consistently across all methods and hours, while the two-hour and three-hour lags (Lag_2, Lag_3) display more variable selection patterns, with inclusion rates changing across different hours and selection criteria.
Temperature and its squared term exhibit time-varying importance: largely excluded at hour 8 but frequently selected at hours 12 and 16, consistent with the thermal sensitivity of energy consumption during warmer parts of the day. Solar shows strong inclusion at hour 8 but diminishing importance at hours 12 and 16, reflecting its natural availability pattern. Week./Hol. indicator maintains high inclusion during morning and midday but loses significance at hour 16, possibly due to convergence between weekday and weekend load patterns during evening hours. Seasonal indicators (Summer, Winter) show moderate inclusion primarily at hour 12, while the Lockdown indicator appears only sporadically, selected solely by RMSE models at hours 12 and 16. The ARMA structure parameters provide further insight into temporal modeling requirements. The weekly lag $y_{t-7}$ is selected across most methods and hours, confirming strong weekly seasonality in energy load. The MA structure varies by hour: $\varepsilon_{t-2}$ is required at hour 8, whereas $\varepsilon_{t-1}$ alone suffices at later hours.
Given the ability of the MSCS to identify models statistically equivalent to the true data-generating process, we study whether our methodology also ensures the selection of well-performing forecasting models. We evaluate this by comparing the predictive performance of models inside the MSCS with the ones left outside. Accuracy is measured using out-of-sample point forecast metrics for hourly load, specifically the root-mean-square error and the mean absolute error:
where $T$ is the number of observations, $R$ is the length of the rolling window and $\hat{y}_{h,t+1|t}$ are the individual hourly load forecasts. Specifically, results are based on a one-step ahead rolling forecasting process with a window size of three years. The initial estimation sample goes from 1 January 2019 to 31 December 2021 and then the forecasting evaluation period starts on 1 January 2023 and ends on 29 September 2023; for a total of 637 out-of-sample point forecasts.
Figure (ref) compares the out-of-sample forecast error distributions (RMSE and MAE) at hours 8, 12, and 16 for models that are included in the MSCS and for those that are excluded. Both metrics convey a consistent message: excluded models exhibit a wide dispersion of forecast errors, including many specifications with substantially poorer performance, whereas models in the MSCS are tightly concentrated in the low-error region. This clear separation indicates that the MSCS isolates a set of model specifications that deliver competitive short-term forecasts of electricity load, and thus provides a useful diagnostic for model choice in forecasting applications. The contrast in the distributions mirrors the MSCS selection mechanism, which is designed to retain models whose losses are statistically indistinguishable from the best-performing specification while discarding models with significantly worse performance.
The results for the LBMs further illustrate the trade-off between parsimony and predictive accuracy. At hour 8, the LBMs rank among the most accurate models, while at hours 12 and 16 they remain highly competitive but lie closer to the right tail of the MSCS distribution. This suggests that, although not always close to the empirical optimum, LBMs offer a good compromise between model parsimony and forecasting performance.
Overall, these findings provide empirical support that the models retained in the confidence set also perform well out of sample, reinforcing the practical usefulness of the MSCS framework for prediction tasks.
Quantifying uncertainty in model selection is becoming a central topic in time series analysis and forecasting. In this paper, we introduced and applied the Model Selection Confidence Set (MSCS) methodology for ARMAX models. The key idea is to move beyond the traditional single-model paradigm and provide, instead, a statistically valid set of admissible specifications. In doing so, the MSCS makes selection variability explicit and provides a principled way to assess structural uncertainty in dynamic environments where misspecification risks are high.
Monte Carlo experiments illustrated the finite-sample behavior of the approach under different ARMAX structures. As sample size increases, the MSCS attains appropriate coverage probabilities and the size of the corresponding LBM sets contracts. Moreover, simpler AR structures tend to be identified more quickly than models that include additional moving-average and exogenous components.
To summarize stability across admissible models, we proposed inclusion importance statistics that quantify how frequently terms appear within the LBMs and across the entire confidence set. In the application to Italian electricity load data we find that selection uncertainty varies systematically over the day: intraday lags, temperature, calendar effects, and solar generation emerge consistently as key drivers, while the overall degree of structural uncertainty is time-varying with morning hours yielding smaller, more stable confidence sets, whereas afternoons admitting a broader range of plausible specifications. This highlights the practical value of acknowledging selection uncertainty when modeling and forecasting electricity load.
Looking ahead, applications to other commodity markets or energy sectors will likely involve much larger predictor pools and richer lag structures, making exhaustive enumeration of the model space infeasible. A natural direction is therefore to combine the MSCS with a conservative screening step that reduces the candidate pool to a manageable subset, treat that reduced specification as the “full model,” and then construct the confidence set over nested models within the restricted space. In addition, because many applications operate at sample sizes where asymptotic approximations may be inaccurate, future work should investigate bootstrap likelihood ratios or other resampling-based calibrations to improve finite-sample coverage and sharpen uncertainty quantification.