Big time series–large collections of observations recorded sequentially over time–have become increasingly common across a wide range of disciplines, including science, business, and industry. The time series could be univariate, or multivariate where several correlated series are analyzed together. Examples of univariate big time series include the business-related complete NASDAQ daily returns over more than 12,000 days, or the engineering-related appliances energy consumption time series of length almost 20,000 in [1]. An example of big multivariate time series is the weekly initial unemployment claims for several U.S. states of length over 1,200 weeks in [3].

Time series analysis typically involves preliminary steps (such as transformations, removing trends or seasonality through differencing), followed by model specification, parameter estimation, diagnostic checking, and forecasting. The focus of this paper is on order identification, a crucial step in model specification. While the method developed here is applicable to both univariate and multivariate time series, this paper focuses exclusively on multivariate time series due to their complexity and to a lack of well-established order identification methods when compared to the univariate case.

From a forecasting perspective, we aim to choose autoregressive (AR) orders that are large enough to capture the underlying structure of the data, yet not so large that estimation error and overfitting become problematic. To discourage overfitting, a penalty term is introduced, as is standard in model selection problems, which are well-studied in the context of regression. We restrict our attention to fitting vector autoregressive models of nonnegative order p,  or VAR(p), using the final prediction error (FPE) criterion. The motivations for this setting include: prediction is often the main goal of time series analysis; any stationary time series can be approximated by an AR(p) process for sufficiently large p; AR(p) processes are simpler and more intuitive than mixed autoregressive-moving average (ARMA) models due to their resemblance to classical regression models, which is especially useful in the multivariate case.

For small time series, identifying the optimal order p is straightforward: one computes the FPE for each candidate p from 1 up to some maximum M, and selects the value that minimizes the criterion. However, this exhaustive approach becomes computationally prohibitive in large time series at the estimation stage. To address this, we introduce an emulator-based optimization method inspired by the design of computer experiments literature [2, 8, 9]. The method proceeds as follows:

1. 1.

   A small subset of equally-spaced AR orders is selected from the finite set of candidate values, providing the initialization set.

2. 2.

   A fast kriging emulator is fit to emulate (i.e. predict) FPE values over the full range of AR orders.

3. 3.

   An efficient global optimization (EGO) algorithm [4] incorporating the fast kriging emulator is used to find the order that minimizes the FPE, thereby identifying a suitable model while drastically reducing computational cost.

Kriging methodology has previously been applied in spatial statistics and computer experiments, and more recently in the context of information criteria (e.g., AIC, BIC) for time series models, as detailed in [10] and [11]. The latter focused on univariate model selection using emulator-based optimization, while the former extended these ideas to multivariate time series. In both cases, the emphasis was on minimizing criteria related to model fit. In contrast, the present work adopts an emulator-based optimization framework targeting predictive performance, as measured by the one-step-ahead mean squared prediction error through the FPE criterion.

Alternative approaches exist for VAR order selection. Classical information criteria such as AIC are commonly used as baseline benchmarks (e.g. [5]) and may have predictive interpretations, but as it will be illustrated in this paper they could overfit the model in some cases, with negative implications for prediction. Penalized likelihood methods, such as LASSO-based VAR (e.g. [6]), introduce sparsity and perform variable selection, but require tuning of regularization parameters and do not directly target a finite autoregressive order. In contrast, the method proposed in this paper directly targets the prediction-oriented FPE criterion, while achieving computational efficiency via kriging-based emulation and efficient global optimization. Maximum likelihood estimation (MLE) would be the ideal estimation method, but it becomes itself computationally prohibitive in large-dimension problems. Instead, we use ordinary least squares (OLS) since the sample sizes considered are large enough for this method to achieve optimal asymptotic properties.

The present work is closely related to previous contributions on emulator-based order selection for autoregressive moving-average models [10, 11], but differs in several important respects. First, whereas prior work focused on univariate or relatively low-dimensional ARMA models and model-fit-based criteria such as BIC, the present paper considers VAR models where the computational burden increases substantially with the dimension K,  sample size T and the candidate order p. Second, we focus on the FPE criterion, which provides a finite-sample measure of prediction error that depends explicitly on K, p, and T, and is particularly relevant when Kp/T is not negligible. Third, we develop an emulator-based optimization framework tailored to this setting, using Gaussian-process emulation and efficient global optimization (EGO) to reduce the number of model fits required for order identification.

Taken together, these contributions extend emulator-based order selection methods to a more challenging multivariate setting and demonstrate that substantial computational gains can be achieved without sacrificing predictive performance.

In Sect. 2, we review the prediction-based criterion used for autoregressive model selection. Section 3 introduces the kriging-based optimization algorithm and Sect. 4 applies it to simulated and real-world time series data. Section 5 concludes with a summary and discussion of potential extensions.

Following results from [5], let \(y_t = (y_{1t}, \ldots , y_{Kt})'\) be a vector of K stationary time series, obtained after applying appropriate transformations and/or differencing to remove trend and seasonal components. A zero-mean VAR(p) model for the transformed data is given by

where \(u_t = (u_{1t}, \ldots , u_{Kt})'\) is a white noise process with mean \(\textbf{0}\) and covariance matrix \(\Sigma _u\).

In practice, the order p of the VAR model is unknown and must be identified from data. This work focuses on selecting the order p that yields good predictions of future values. Let T be the sample size, K the number of component series, and \(\tilde{\Sigma }_u(p)\) the maximum likelihood estimator of the covariance matrix for a VAR(p) model.

Following [5, Eq. 4.3.23], the one-step-ahead prediction error variance for a VAR(p) model yields the Final Prediction Error (FPE) criterion,

Although FPE and AIC are asymptotically equivalent as \(T \rightarrow \infty \), the equivalence relies on the ratio Kp/T being negligible. As shown in [5, p. 148], their exact relationship is

a difference that depends only on T, K, and p.

Tables 1, 2, 3 and 4 show that this difference can be substantial when both K and p are large relative to T. For moderate sample sizes (e.g., \(T = 500\) or \(T = 1000\)), the discrepancy grows systematically with increasing dimensionality K and autoregressive order p, while it becomes negligible only when T is much larger than Kp. In such high-dimensional, moderate-T settings, which is precisely the context of our application, FPE retains its interpretation as a one-step-ahead prediction error estimate, whereas AIC can deviate markedly. For this reason, FPE provides a more appropriate criterion for prediction-focused order selection in the regimes emphasized in this paper.

This paper adapts the methodology developed in [10, 11] to the context of prediction-based model selection using the log-transformed final prediction error (LFPE) criterion for multivariate AR(p) models. We consider the AR order p as the input to a computer model \(\text {LFPE}(p)\) and aim to efficiently identify the value of p that minimizes this criterion.

Let \(\{p_1, \ldots , p_D\}\) denote a subset of AR orders sampled from the discrete domain \(\{1, 2, \ldots , M\}\) using equally spaced points to ensure the points are well spread across the domain. For each sampled order \(p_i\), we compute \(\text {LFPE}(p_i)\), yielding the response vector:

We assume this vector follows a multivariate normal distribution with constant mean \(\gamma \) and covariance matrix \(\tau ^2 \textbf{R}_D\), where \(\textbf{R}_D\) is the correlation matrix defined using the Matérn 5/2 kernel:

with \(\eta > 0\) as a covariance parameter. The parameters \(\gamma \), \(\tau ^2\), and \(\eta \) are estimated by maximizing the likelihood of \(\textbf{LFPE}\) conditioned on the time series data. While other covariance kernels can be considered (e.g., power exponential), we use the Matérn 5/2 kernel because it produces low-order smooth Gaussian-process realizations, which aligns well with the behavior of the observed LFPE values.

Given an arbitrary input \(p \in \{1, 2, \ldots , M\}\), the kriging predictor for \(\text {LFPE}(p)\) is

with associated prediction variance

where \(\hat{\gamma }\) denotes the maximum-likelihood estimate of the constant mean in the Gaussian-process model for LFPE. The vector \(\textbf{r}_0\) contains the Matérn 5/2 kernel evaluations between the input p and the design points \(\{p_1,\ldots ,p_D\}\) given by

It is important to distinguish between the two types of prediction in this paper: (i) the time-series prediction and its associated error that define LFPE, and (ii) the kriging prediction of LFPE, which provides a fast emulator of the computed LFPE values.

Minimizing only the kriging predictor \(\widehat{\text {LFPE}}(p)\) can yield suboptimal solutions, since predictive uncertainty is ignored. To address this problem, we maximize the Expected Improvement (EI) criterion from [4]:

where \(\text {LFPE}_{\min } = \min \{\text {LFPE}(p_1), \ldots , \text {LFPE}(p_D)\}\), and \(\Phi (\cdot )\) and \(\phi (\cdot )\) denote the standard normal CDF and PDF, respectively. The equivalence between minimizing the objective function and maximizing EI is explained in [4] in general, for any computer model.

This EI criterion balances **exploitation** (evaluating near the current minimum of log(FPE)) and **exploration** (evaluating where uncertainty is high). The resulting Efficient Global Optimization (EGO) algorithm is defined as follows:

1. 1.

   Select D initial equally-spaced AR orders \(\{p_1, \ldots , p_D\}\) from \(\{1, 2, \ldots , M\}\) and compute their corresponding \(\text {LFPE}(p_i)\) values.

2. 2.

   Construct the kriging predictor \(\widehat{\text {LFPE}}(p)\) and variance \(\text {Var}(p)\) over the domain.

3. 3.

   Choose \(p_{D+1}\) to maximize the expected improvement \(\text {EI}(p)\).

4. 4.

   Augment the sample with \(p_{D+1}\) and its associated \(\text {LFPE}(p_{D+1})\) value.

5. 5.

   Repeat steps 2–4 until three consecutive points are selected, as this is indicative of visiting repeatedly (during the exploitation stage) the neighborhood where the global minimum is located.

Since the search space \(\{1, 2, \ldots , M\}\) is finite, EGO is guaranteed to terminate after a finite number of steps. In practice, however, the stopping criterion 5 above is typically achieved well before all values are evaluated, resulting in substantial computational savings. The criterion based on consecutive points is heuristic but reflects convergence of the search to a narrow region and was found to perform well in practice. Stopping after three consecutive points prevents the algorithm from wasting computational time once the acquisition function has essentially flattened out at the bottom of the “bowl.” More details on kriging-based optimization, expected improvement and the EGO algorithm in the general case can be found in [4].

The dominant computational component in the proposed method is fitting the VAR(p) models, which is substantial for large sample size T and large dimension K. In contrast, evaluating LFPE once the VAR(p) model is fitted, and fitting the Gaussian process emulator once the LFPE values are computed, are computationally modest. Consequently, the proposed EGO-based approach achieves computational savings primarily by reducing the number of VAR model fits required for order selection.

To generate realistic high-order VAR(p) models that satisfy stationarity and avoid unrealistically large coefficients, we adopted the following two-step approach. First, we simulated a 30-dimensional stationary vector moving average of order 1, or VMA(1), process,

where \(u_t\) is a white-noise, zero-mean vector with identity covariance matrix, and \(\Theta \) is chosen so that all eigenvalues lie inside the unit circle. Second, we fit a VAR(p) model to the simulated VMA(1) data and used the estimated coefficient matrices \(\left\{ \hat{A}_1,\ldots ,\hat{A}_p\right\} \) as the “true” AR matrices in our simulation experiments. This procedure ensures stationarity because a VAR fitted to a causal VMA(1) process inherits the stability properties of the underlying data-generating mechanism. It avoids hand-tuning AR coefficients, prevents violations of stationarity, and produces realistic high-order VAR structures whose coefficients decay in a manner consistent with empirical applications.

In the simulation study, we consider two representative orders, \(p = 12\) and \(p = 28\). To evaluate the performance of the proposed order selection method, we simulate a vector of \(K = 30\) stationary time series of length \(T = 1600\) and \(T=2500\) from known AR(p) models, resulting in four cases. These values were chosen to reflect moderately high and very high autoregressive orders within the candidate domain \(\{1,\ldots ,40\}\) and moderately large and very large sample sizes, allowing us to assess the performance of the emulated FPE method in scenarios where capturing long-range dependence is challenging and the computational cost of full enumeration is substantial as well as illustrate settings where the AIC and LFPE are roughly the same across the domain and settings where the AIC and LFPE differs as the AR order increases.

For each point of the candidate domain, we used the VAR function from the vars package in R, employing OLS to fit a VAR(p) model. Other R packages (e.g. MTS) that implement MLE could be used, but turned out to be computationally prohibitive despite MLE having a greater degree of accuracy than the more traditional OLS. However, in our implementations, the sample size T is large enough that OLS is sufficiently accurate.

To assess the reliability of our findings, we first used the more traditional exhaustive grid search method i.e. fitting a VAR model to all points of the \(\{1,\ldots ,40\}\) domain and then compared the required computing time and accuracy to our proposed method. For both methods, and for all four cases in each of the two methods, 50 independent simulations were performed.

The plots of AIC versus LFPE over the grid \(\{1,\ldots ,40\}\) for all four cases are given in Fig. 1. As expected, in the larger \(T=2500\) cases, AIC and LFPE both consistently identify the correct order, whearas in the smaller \(T=1600\) cases, LFPE consistently identifies the correct order while AIC does not. As summarized in Table 6, both AIC and LFPE have high prediction interval coverage in the larger \(T=2500\) cases. On the other hand, as shown in Table 5, LFPE has substantially better coverage than AIC does in the smaller \(T=1600\) cases. Computing time on the grid \(\{1,\ldots ,40\}\) was substantially longer for the \(T=2500\) cases than the \(T=1600\) cases.

We then proceeded to employ the proposed method detailed in Sect. 3. To construct the initial design for the Gaussian-process emulator, we select \(D=6\) equally spaced orders spanning the full domain. Specifically, with spacing \(\Delta =8,\) the initial subset is \(\{p_1,\ldots ,p_6\}=\{1,\Delta ,2\Delta ,\ldots ,5\Delta \}=\{1,8,16,24,32,40\}.\) The choice \(D=6\) provides good space-filling coverage while keeping the number of expensive evaluations small. Very small D (e.g. \(D=3\)) leads to inaccurate initial emulator impacting negatively the subsequent EGO algorithm, while a very large D (e.g. \(D=15\)) may waste VAR model fits as the EGO algorithm may stop earlier. For each element of the subset, we use OLS to fit a VAR(p) model, and compute the corresponding LFPE values.

We then use the DiceKriging package in R, employing the default Matérn 5/2 covariance function, to construct a kriging emulator model that emulates the LFPE values over the entire input domain. The Matérn 5/2 kernel provides a practical balance between smoothness and numerical stability, and it reflects the low-order smoothness consistent with the full-grid values observed in both the simulations and the real application. Alternative kernels such as Gaussian or exponential may over-smooth or under-smooth the LFPE values, respectively.

To identify the optimal order, the DiceOptim package is used to perform EGO, iterating until three consecutive points are selected. Further details on the DiceKriging and DiceOptim methodology can be found in [7].

Figures 2, 3, 4 and 5 illustrate a representative simulation run for each of the four cases cases, respectively. The left panel shows the full set of LFPE values, the initially sampled points (blue), the EGO iterations (red), and the final selected order (green). The right panel tracks the predicted LFPE values during the EGO iterations. This example demonstrates the ability of the emulated optimization approach to quickly converge to a region near the true model order using only a few evaluations, significantly reducing the overall computational burden.

The objective function appears relatively smooth over the lag space, which allows EGO to identify the optimal region in just three or four iterations.

To illustrate the proposed method on real-world data, we analyze weekly initial unemployment claims from 1986 to 2019 for \(K=20\) U.S. states obtained from the Federal Reserve Economic Data (FRED) database [3], accessed via the quantmod R package. This resulting in over \(T=1200\) overlapping time points. To address the variance concerns, we applied a log transformation followed by first-order differencing, and we also excluded the post-2020 period to avoid the COVID-related volatility.

To show that this is a setting where AIC tends to overfit the data as opposed to LFPE, we first fit VAR(p) models across the full range \(p \in \{1,\ldots ,40\}\). This exhaustive search required approximately 120 s. As shown in Fig. 6, LFPE selects \(p = 25,\) while AIC continues decreasing and selects \(p = 40,\) suggesting overfitting. This provides a clearer contrast between the two criteria, with FPE appearing more reliable.

To validate our hypothesis about AIC’s overfitting, we conducted an out-of-sample prediction study using a three-step expanding window pseudo-out-of-sample (POOS) forecast for the first three weeks of 2020. This approach provides a robust evaluation of predictive performance across volatile post-holiday data. Over the three steps, the LFPE model (\(p=25\)) achieved an average coverage rate of 80.0%, significantly outperforming the AIC model (\(p=40\)), which achieved an average coverage of 68.3%. While AIC consistently produced slightly tighter average interval lengths (~0.511 vs ~0.537), its poor coverage perfectly demonstrates “false precision.” The \(p=40\) AIC model is overconfident and overfit to the training noise, causing its prediction intervals to miss the true out-of-sample values during volatile periods. LFPE’s properly penalized intervals ultimately provide a safer, more robust forecast.

We then applied the EGO approach using LFPE. According to Fig. 7, it selected \(p = 26\) in 5 iterations and required about 40 s, which is a substantial reduction in computation time while maintaining close agreement with the exhaustive search. This example highlights the ability of the emulated optimization approach to converge quickly to the true model order using only a few evaluations, substantially reducing computational effort.

All computations were conducted in R on a Windows 11 PC equipped with 32 GB RAM and an Intel Core Ultra 7 155 H CPU running at 3.80 GHz.

This paper extends the computationally efficient approach in [10, 11] by applying kriging-based emulator modeling and global optimization to prediction-based order selection for multivariate autoregressive models. By evaluating the log-transformed final prediction error (LFPE) at a strategically chosen subset of autoregressive orders, we construct a kriging emulator to approximate the criterion across the entire search space. Efficient global optimization (EGO) is then applied to locate the order minimizing LFPE, using a predefined stopping rule.

Simulation studies and real data analysis show that the proposed method achieves substantial reductions in computation time–often by a factor of three or more–while accurately recovering the true order. In the real data example of weekly initial unemployment claims across 20 U.S. states, the proposed method reduced the computation time from 120 s for the exhaustive search down to approximately 40 s, while identifying an optimal order that maintained robust out-of-sample prediction coverage.

The accuracy and efficiency of the method depend on the choice of initial design points and the stopping criterion. Future work may also consider adaptive sampling strategies or alternative stopping criteria (e.g., requiring four consecutive points) to further balance accuracy and computation time.