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.
79,118 characters · 25 sections · 66 citation commands
Long-term prediction intervals with many covariates
Keywords: Forecasting, Heavy-tailed distribution, Long-range dependence, Electricity prices, Bootstrap, Time-aggregation
Prediction intervals (PI hereafter) help forecasters to access the uncertainty concerning the future values of time series. This has been a central topic of interest in the analysis of time series from both theoretical and algorithmic point of view. The $k$-step ahead prediction of a time-series after observing $X_1, \ldots, X_t$ has been discussed from many different perspectives. While estimation of the conditional mean $\mathbb{E}(X_{n+k}|X_1,\ldots,X_n)$ has received most of the attention since 11; other approaches were also developed significantly in the last few decades, see 20, 29, 18, 48, 31 33, 8, 17, 35. Nonetheless in the eye of an applied scientist, evaluation of prediction intervals ( PIs henceforth) is challenging c93,ct03, since a higher empirical coverage probability often comes with the cost of larger width, thus less precision. In this paper, we focus on the construction and evaluation of PIs for a single univariate series using information from many predictors under a possibly high-dimensional linear model with a general error process. Moreover, the target for which we construct the PIs is not a single future value of the univariate series, but a sum of its values over the entire forecasting horizon. This setup has practical applications in telecommunication and data services, energy, and finance. The characteristic feature of these sectors is that data about sales, prices, or returns are collected with hourly frequency while the management is interested in predictions of aggregated volumes over horizon spanning several weeks or months. In the volatility forecasting of stock returns using popular ARCH/GARCH models, often aggregated (over time) mean squared error (AMSE) is an evaluation metric (See starica2003garch,subbarao2008,sayar2020 etc.). When it comes to predicting energy or electricity consumption, policy makers want to predict future usages for entire one month or several months ahead. The aggregated predictions help them in setting the price of such commodities. See some of the applications mentioned in zhou10. Time-aggregation is also useful to arrive at strategic decisions made by trust funds, pension management and insurance companies, portfolio management of specific derivatives kw13 and assets bky16 among others.
Towards including external regressors, lfn15 found that the inclusion of disaggregated wind speed and temperature (measured at more than 70 weather stations across Germany) leads to improvements of forecasts for EPEX SPOT hourly day-ahead electricity prices.\footnote{In electricity price forecasting (EPF), one has to control for weather conditions, local economy, and environmental policy kr05,hrz12. Additionally, EPF is challenging due to complex seasonality (daily, weekly, and yearly), heteroscedasticity, heavy-tails, and sudden price spikes w14.}. We adapt their framework to allow for long-run dependence which is important for long-horizon and medium-horizon forecasting. Similar applications can be found in in power portfolio risk management, derivatives pricing, medium-term and long-term contract evaluation, and maintenance scheduling. To challenge our forecasts, we provide out-of-sample comparison with several forecasting methods, which include Bayes PIs of mw16 and bootstrap PIs obtained from methods such as exponential smoothing, neural networks, and regression with auto-correlated errors implemented in the R-package forecast hk08.
An interesting advantage of predicting a sum of future values instead of just a single $k$-step ahead prediction is it allows for some weak law of large numbers to kick in if the horizon length of prediction is large. Towards the methodological development and corresponding theoretical advancement, we explore the properties of the quantile-based PIs of zhou10 in a high-dimensional regression setting. First we identify that zhou10 uses linear processes in their modelling of innovation processes. This limits the scope of applicability. Using the idea of predictive density and the functional dependence framework as proposed in a seminal paper by wu05, we extend the quantile consistency results to a large class of nonlinear processes which includes thresholded/transition autoregressive ltd03 processes and some nonlinear version of GARCH processes.
On the regression front, the quantile consistency results obtained so far were limited to only a low-dimension regression setting. This builds a strong motivation to explore a high-dimensional regression setting as this would allow us to incorporate several predictors that better explain the variability in the data. In particular, we put a special focus on the scenario where the number of predictors grows much faster than the sample size ($\log p=o(n)$). Using the popular LASSO estimator, we were able to provide theoretical guarantees for these PIs. One interesting contribution of this paper is how our consistency rates explicitly depict the price of having short or long-range dependence and lighter or heavier tails. We also show sharp consistency results for stochastic design in the high-dimensional regime which can be seen as a contribution that is important on its own. Thus on one hand, this paper advances the theory of future-aggregated prediction by enlarging the scope of time-series; at the same time, it allows for many covariates in a high-dimensional regression setting. In particular, this paper can be seen to provide some theoretical justification for how a simple LASSO-based prediction routine can sustain such general time-series dependence.
Besides the theoretical contributions, we tackle the practical validity of PIs in the scenario of long-horizon forecasting with only a relatively short sample. One new discovery was that the original quantile-based PIs fail when the length of time-horizon to be aggregated for forecasting grows compared to the length of the past data. In light of this shortcoming, we employ a bootstrap-assisted step to improve the out-of-sample coverage probability. Apart from a conjectural viewpoint justifying our bootstrap procedure, these simulation results support our choice of PIs for the final out-of-sample experiment using real-world data on hourly electricity prices.
The rest of the paper is organized as follows: Section (ref) shows the construction of quantile-based PIs and also details a novel bootstrap adjustment to improve coverage. Section (ref) states the quantile consistency results in an exhaustive number of cases: the error/innovation process can be linear or non-linear, short-range or long-range and can have light or heavy tail. Building on the consistency for just the error process, we use Section (ref) to consider consistency results for LASSO-fitted residuals. Apart from the traditional fixed design, we also provide discussion on stochastic design in this section. Moving on, the simulations are shown in details in Section (ref) where we compare OLS, LAD and LASSO for both low-dimension and high-dimension regression setup. Section (ref) is used to analyze EPEX spot electricity using several methods including ours. For contrasting our method with the competing ones in literature we use a Pseudo-out-of-sample (POOS hereafter) approach. Finally Section (ref) provides concluding remarks. We defer some theoretical results and all proofs of theorems to the appendix.
We now introduce some notation. For a random vector $Y$, write $Y \in \mathcal{L}_p$, for $p > 0$, if $\|Y \|_p := E(|Y |^p ) ^{1/p} < \infty$. For the $\mathcal{L}_2$ norm write $\|\cdot \| =\| \cdot \|_ 2$. Throughout the text, $c_p$ denotes a constant that depends only on $p$ and $c$ denotes a universal constants. These might take different values in different lines, unless otherwise specified. Then, $x^+=\max(x,0)$ and $x^{-}=-\min(x,0)$. For two positive sequences $a_n$ and $b_n$, if $a_n/b_n \to 0$, write $a_n =O(b_n)$. Write $a_n \lesssim b_n$ if $a_n \leq c b_n$, for some $c<\infty$. For a random variable sequence $X_n$ and a possibly random or non-random sequence $Y_n$, we say $X_n=o_{\mathbb{P}}(Y_n)$ if $X_n/Y_n \to 0$ in probability. Respectively if $X_n/Y_n$ is a tight random variable we say $X_n=O_{\mathbb{P}}(Y_n)$. The $d$-variate normal distribution with mean $\mu$ and covariance matrix $\Sigma$ is denoted by $N(\mu, \Sigma)$. Denote by $I_d$ the $d \times d$ identity matrix. For a matrix $A= (a_{ij})$, we define its element-wise $\ell-\infty$ norm as $|A|_{\infty}= max_{ij} |a_{ij}|$ and its $\ell_1$ norm as $|A|_1= max_{j} \sum_i |a_{ij}|$.
Suppose an univariate target time series $y_i,i=1,\ldots,n$ follows a regression model
where $e_i$ stands for the mean-zero temporally-dependent error process. Assuming model (ref), we wish to construct PI for $y_{n+1}+\ldots+y_{n+m}$ after observing $(y_i,\bm x_i);i=1,\ldots,n$. We first discuss the scenario $\bm \beta=0$, i.e., $y_i=e_i$ is a zero-mean noise process. This serves as a primer to the high-dimensional regression problem and introduces a basic quantile-based method of constructing the prediction interval.
We will appropriately define our notion of short-range and long-range and light-tailed and heavy-tailed distribution later in Section (ref). Loosely speaking, short-range dependence stands for faster decay of dependence as the lag increases where long-range dependent process kills dependence slowly. Heavy-tailed refers to the scenario where the error process does not possess even the finite second moment. Depending on the strength of dependence and the tail behavior, zhou10 proposed two different types of PIs for $m$-step ahead aggregated response $e_{n+1}+\ldots+e_{n+m}$. We present them here, with some suggestive modifications especially for the first approach where we estimate the long-run variance differently.
If the process $e_t$ shows short-range dependence and light-tailed behavior, then in the light of a quenched central limit theorem, zhou10 proposed the following PI for $\frac{1}{m}(e_{n+1}+\cdots+e_{n+m})$,
where $\sigma$ is the long-run standard deviation (sd) of $e_t$. However, since $\sigma$ is unknown, it must be estimated. One can estimate the long-run variance $\sigma^2$ of the $e_i$ process using the sub-sampling block estimator dfsvw13
with block length $l$ and number of blocks $\kappa=\lceil n/l\rceil$ in order to obtain $100(1-\alpha)\%$ asymptotic PI $$[L,U]=\pm \hat \sigma Q^t_{\kappa-1}(\alpha/2) \sqrt{m},$$ where $Q^t_{\kappa-1}$ is student-t quantile with $\kappa-1$ degrees of freedom.
A substantially more general method that can account for long-range dependence or heavy-tailed behavior of the error process uses the following quantile
The PI in this case is
This approach enjoys reasonable coverage for moderate rate of growth for $m$ compared to sample size $n$.
Now assume the scenario where $\bm \beta$ is possibly non-zero in ((ref)). For $p<n$, with dependent but linear error process, after estimating $\bm{\bm \beta}$, zhou10 constructed PI as
where $\hat{e}_i=y_i-\bm{x}^{\mkern-1.5mu\mathsf{T}}_i\hat{\bm{\bm \beta}}$ are the regression residuals. Regarding the choice of estimator for $\bm{\bm \beta}$, if the error process shows light tailed behavior and short-range dependence, we typically use OLS estimator $ \bm{\hat{\bm \beta}}=\mathop{\rm argmin} \sum_i (y_i-\bm{x}^{\mkern-1.5mu\mathsf{T}}_i \bm{\bm \beta})^2 $. For heavy-tailed or long-range dependent errors huber, it is better to use robust regression with general distance $\rho(\cdot)$ and $\bm{\hat{\bm \beta}}= \mathop{\rm argmin} \sum_i \rho (y_i-\bm{x}^{\mkern-1.5mu\mathsf{T}}_i \bm{\bm \beta}).$ Examples of distance include the $\mathcal{L}^q$ regression for $1 \leq q \leq 2$. In this paper our focus is on a LASSO-type least square estimation (cf. (ref)) in presence of possibly dependent and nonlinear errors where the number of covariates $p$ is much larger than the sample size $n$.
Consider the high dimensional regression situation i.e. $p \gg n$. Here we use the popular LASSO estimator
with the penalty coefficient $\lambda$. Then we get PIs ((ref)) with $\bm{\hat{\bm \beta}}$ replaced by the LASSO estimator. One of the key contributions of our paper lies in the fact that the intuitive and computationally fast lasso estimator adapts well for the time-aggregated response scenario but still provides theoretical consistency guarantee as we show in Section (ref).
Note that, the PI in ((ref)) requires the future values of covariates $\bm x_i$ namely $\bm{x}_{n+1},\ldots, \bm{x}_{n+m}$. zhou10 discussed the scenario where these predictors are trigonometric and thus perfectly predictable. For practical implementation in econometrics, it is more appealing to consider stochastic design so that one can allow for the covariates $\bm x_i$ to be observed random variables as well. This naturally raises the question of how to predict the future values of the covariates. For the theoretical part, we shed some light on consistency results for the stochastic design in Section (ref). However, in the implementation in Section (ref), we choose to simply use a year-over-year mean forecast for the weather covariates.
In our actual implementation for real data, we propose a bootstrap adjusted version of the algorithm in Section (ref). In the stationary bootstrap, one randomly draws a sequence of starting point uniformly from $\{1,2,\ldots,n\}$ and a sequence of Geometric random variable $L_1,L_2,\cdots$. Depending on the values of $L_i$ and the starting value, we draw a consecutive $L_i$-length block and then finally concatenate these blocks together to arrive at a replicated series of length $n$. Then finally we look at the final $m$ for each such series, whereas it is understandable that it was not particularly important to only look at the last $m$ since the starting point is uniform. We conjecture that this ensures consistency of the quantiles despite providing some more dispersion to the average. This extra dispersion arises as this method concatenates blocks containing same elements with a nontrivial probability and thus adding those covariances in the dispersion of the overall average. We postpone a rigorous theoretical justification of this innovative bootstrap technique to a future work since this paper focuses more on the exploration of the performance of LASSO fitted residuals. For the convenience of the readers we summarize the algorithm down below:
Our theoretical results are concerned with the consistency of the usual quantiles from the original series. But we conjecture that the stationary bootstrap technique to obtain the replicated series retains the asymptotic dependence structure of the original series and thus the quantiles of the $m$-length average from the original series and that from the final $m$ of the replicated series are close to each other. Additionally, using the Gaussian kernel density to obtain the kernel quantile estimator (See sm90) further improves the performance in prediction. These adjustments are supported by the empirical evidence given in chudyold in a univariate setup.
Next, we exhibit some theoretical consistency results addressing various different dependency cases, tail-decay and high dimensional regression in the next two sections.
Before moving on to a more general discussion with a large number of covariates compared to sample size, we use this section to discuss a primer without covariates i.e. our response here is just a mean-zero error process. This section also elaborately describes the model specifications on the error process i.e. whether the process is short/long-range dependent or has light/heavy tails. Under short-range dependence, if window size $m$ is long enough then the dependence of $y_{n+1}+\ldots+y_{n+m}$ on $y_1,y_2,\ldots, y_n$ diminishes and the conditional distribution of $(y_{n+1}+\ldots+y_{n+m})/\sqrt{m}$ given $y_1,\ldots,y_n$ is almost similar to the unconditional distribution and thus one can obtain a simple central limit theorem to quantify the uncertainty surrounding the prediction. Consider $(e_i)$ a mean-zero stationary process and let $S_m=e_1+\ldots+e_m$. Wu04 proved that for $q>5/2$, the condition
gives the a.s. convergence
where $\|\cdot \|_2$ denotes the $L_2$-norm, $\Delta$ denotes the Levy distance, $\mathcal{F}_i$ is the $\sigma-$field $\sigma(\ldots, e_{i-1},e_i)$, $m \to \infty$ and $\sigma^2= \lim_{m \to \infty}\|S_m\|_2^2/m$ is the long-run variance. We start by collecting an asymptotic normality result from chudyold. When the error process has possibly long-range dependence such asymptotic normality fails. Keeping the linear structure intact, we provide some empirical consistency result for the quantile-based methods from Section (ref). These results depict all possible scenarios of heaviness of tail and range of dependence. Finally we extend these results to the more practically applicable non-linear case under the functional dependence framework. Throughout this section, the observed process $e_i$ are either direct linear sum of independent and identically distributed (i.i.d. henceforth) innovations that are unobserved or we assume $e_i$ to be a general non-linear function of these. For the latter, we define some tractable moment-based coupling measure which we call functional dependence measure.
Assuming linearity of the mean-zero noise process $e_i$ in the following manner
for some $p>2$, it is easy to derive the following conditions on $a_i$ to ensure the convergence in ((ref)). The proof can be found in the appendix of chudyold.
Next, as outlined above we proceed for a quantile consistency result for linear processes to validate our proposed methods. We start by formally defining short-range and long-range dependence for linear processes with the representation in ((ref)). For short-range dependence we assume
$$\text{(SRDL): }\sum |a_i|<\infty.$$
For long range dependence, we first revisit the definition of slowly varying function (s.v.f.). A function $g(\cdot)$ is called a s.v.f. if for all $a>0$, $\lim_{x \to \infty}g(ax)/g(x)=1$. Assume
$$\text{(LRDL($\gamma$))}: l^*(i)=a_i i^{-\gamma} \text{ is a slowly varying function (s.v.f.)} $$
for some $q<\gamma<1$ where $1/q= \sup \{t:\mathbb{E}(|\epsilon_j|^t)<\infty\}$. Note that, the definition of $(\text{LRDL}(\cdot))$ also takes into account the possible heavy-tailed distribution of $\epsilon_j$. We also assume $\epsilon_i$ admits a density $f_{\epsilon}$ and
$$\text{(DENL):} \sup_{x \in \mathbb{R}}(f_{\epsilon}(x)+|f'_{\epsilon}(x)|)<\infty.$$
For a fixed $0<u<1$, let $\hat{Q}(u)$ and $\tilde{Q}(u)$ denote the $u$-th sample quantile and the actual quantile of $\tilde{S}_i$, $i=m,\ldots,n,$ respectively, where
and
where Case 1-2 denotes SRDL holds and $\mathbb{E}(\epsilon_j^2)<\infty$ or $\mathbb{E}(\epsilon_j^2)=\infty$ respectively whereas Case 3-4 stands for $\text{LRDL}(\gamma)$ holds and $\mathbb{E}(\epsilon_j^2)<\infty$ or $\mathbb{E}(\epsilon_j^2)=\infty$ respectively. We collect the following theorem from zhou10 for the rates of convergence of quantiles depending on the nature of the error process in terms of tail behaviour and dependence:
A technical fact about these results is that their proofs are heavily dependent on the linear structure. We believe that it is an important task to extend these to the more general non-linear scenarios but it comes at the cost of some more abstraction in how these nonlinear processes are posited. We use a very general framework of functional dependence form wu05 to first describe the class of non-linear process and then provide analogs of central limit theorem and quantile consistency results.
Economic and financial time series are often subject to structural changes and thus the linear models often fail to capture the diverse range of data generating process. Comparatively possible non-linear class of time-series models is significantly richer in its scope. Useful nonlinear time-series models, among many, include the regime-switching autoregressive processes, which assume that the series change their dynamics when passing from one regime to another and neural-network models. We provide extension of the asymptotic normality result (ref) to nonlinear error process. We lay down the definition of non-linearity through the following structure: Let $e_i$ be a stationary process that admits the following representation
where $H$ is such that $e_i$ is a well-defined random variable, $\epsilon_i, \epsilon_{i-1}, \ldots$ are i.i.d. innovations and $\mathcal{F}_i$ denotes the $\sigma$-field generated by $(\epsilon_{i}, \epsilon_{i-1}, \ldots)$. One can see that it is a vast generalization from the linear structure in $H$. In order to derive a result similar to (ref) but in the non-linear regime, we define the following functional dependence measure for $e_i$ in (ref), by which we follow wu05's framework to formulate dependence through coupling:
where $\mathcal{F}_{i,k}$ is the coupled version of $\mathcal{F}_i$ with $\epsilon_k$ in $\mathcal{F}_i$ replaced by an i.i.d. copy $\epsilon_k'$, $\mathcal{F}_{i,k}= (\epsilon_i, \epsilon_{i-1}, \ldots, \epsilon_k', \epsilon_{k-1}, \ldots )$ and $e_{i,(i-j)}= H(\mathcal{F}_{i,(i-j)})$. Clearly, $\mathcal{F}_{i,k}= \mathcal{F}_i$ is $k>i$. As wu05 suggests, $\|H(\mathcal{F}_i)- H(\mathcal{F}_{i, (i-j)}) \| _p$ measures the dependence of $e_i$ on $\epsilon_{i-j}$. This dependence measure can be seen as an input-output system and a natural analogue of the linear coefficient at a certain lag. It facilitates mild moment conditions on the dependence of the process which are easily verifiable compared to the more popular strong mixing conditions. Define the cumulative dependence measure
which can be thought as cumulative dependence of $(e_j)_{j \geq k}$ on $\epsilon_k$.
Apart from the asymptotic normality result from Theorem (ref), we will next show that the estimated quantiles from the moving blocks are consistent. However, the same notion of non-linear dependence falls short of establishing quantile consistency due to a technical reason and thus we resort to first define some form of predictive dependence.
For a general nonlinear process with possibly heavy tails and long-range dependence, the central limit theorem fails. We present one of the main results in this paper by showing the empirical quantile consistency that validates PIs of the form (ref). Towards that, we need to control the latent dependence of $e_i$ on $\epsilon_{i-j}$ keeping in mind a representation like ((ref)). Therefore, we introduce the predictive density-based dependence measure. Let $\mathcal{F}'_k = (\epsilon_k, \ldots , \epsilon_{1}, \epsilon_0', \epsilon_{-1}, \ldots ,),$ be the coupled shift process derived from $\mathcal{F}_k$ by substitution of $\epsilon_0$ by its i.i.d. copy $\epsilon_0'$. Let $F_1(u|\mathcal{F}_k) = P\{H(\mathcal{F}_{k+1}) \leq u|\mathcal{F}_k\}$ be the one-step ahead predictive or conditional distribution function and $f_1(u|\mathcal{F}_k) = d F_1(u|\mathcal{F}_k)/ du,$ be the corresponding conditional density. We define the predictive dependence measure
The quantity $\psi_{k,q}$ ((ref)) measures the contribution (read change in $q$-th norm) of $\epsilon_0$, the innovation at step 0, on the conditional or predictive density at step $k$. We shall make the following assumptions:
The (DEN) implies that the marginal density $f(u) = \mathbb{E}{f_1(u|\mathcal{F}_0)} \leq c_0.$ Recall the sufficient conditions for the linear cases in zhou10 were based on the coefficients of the linear process. Here, the conditions for both short-range and long-range dependent errors here were transferred onto predictive dependence measure. We assume:
where $1/q= \sup \{t:\mathbb{E}(|\epsilon_j|^t)<\infty\}$. For a fixed $0<u<1$, let $\hat{Q}(u)$ and $\tilde{Q}(u)$ denote the $u$-th sample quantile and actual quantile of $\tilde{S}_i$; $i=m,\ldots,n,$ where
and
Here, Case 1-2 denotes $SRD$ holds and $\mathbb{E}(\epsilon_j^2)<\infty$ or $\mathbb{E}(\epsilon_j^2)=\infty$ respectively whereas Case 3-4 stands for $LRD(\gamma)$ holds and $\mathbb{E}(\epsilon_j^2)<\infty$ or $\mathbb{E}(\epsilon_j^2)=\infty$ respectively.
Then we have following rates of convergence of quantiles depending on the nature of the error process in terms of tail behaviour and dependence:
Note that the results presented above talks about the scenario without any predictor. For the low dimensional regression case, the proofs follow exactly similar to how zhou10 proceeded by showing $$\sup_{i \leq n}|\hat{e}_i-e_i|=O_{\mathbb{P}}(\Pi(n)) $$for some suitably chosen small $\Pi(n)$. We skip the details here and directly move on to the scenario where the number of predictors far outgrows the number of past time-points available.
In this section, we focus on results concerning theoretical guarantees of the prediction intervals for high-dimensional regression where the error process is possibly nonlinear and shows temporal dependence. Note that all the results from this section can be easily extended to situations where the error process has exponentially decaying tails such as sub-exponential and sub-gaussian, but we restrict ourselves to conditions of only finitely many moments which is a significantly weaker assumption than that is present in the literature.
The results for this subsection heavily depend on optimal concentration inequalities for $S_{n,b}=\sum_{i=1}^nb_ie_i$ that have not been established for possibly nonlinear dependence before this. Since this is just a proof technique tool we postpone the concentration results here. However, as one of our final aim is to point out the price of dependence in a wide range of general scenarios, it is important to define dependence-adjusted norm to be able to state the concentration and consistency results. Recall ((ref)) and assume the short-range dependence, $\Theta_{0,q}<\infty$ holds with $q$ being less or more than 2 depending on the tail-behavior of the error process. Further, we define dependence adjusted norm, for $\alpha>0$,
It is easy to note that, finiteness of the dependence-adjusted measure is a stronger ask than the finiteness of $\Theta_{0,q}$. Next, we show that for short-range dependent nonlinear error processes the error bounds obtained in Theorem (ref) remain intact under a proper choice of the sparsity condition.
For the model in ((ref)), we first assume that the $\bm x_i$'s are fixed and the future $\bm x_i$'s are known. Under this setting, we next show the quantile consistency for the nonlinear process. Note that a very similar result can be shown for the linear process with the error process admitting a simpler representation ((ref)), however we skip writing that as a separate theorem here to avoid repetitiveness.
One can note the $\sqrt{\log p/n}$ term we have in our definitions for $r=\lambda/2$. This allows us to capture the ultra-high dimensional scenario where $\log p =o(n)$, the usual benchmark in the high-dimensional literature. The additional terms involving $\|e_{.}\|_{.,\alpha}$ are due to the dependence present in the error process. The sparsity condition for the light tail case, i.e. $q \geq 2$, in ((ref)) in the view of the choice of $r$ in ((ref)) can be written as: For $q\geq2$
for the choice of $m=o(n^{1/3})$. Thus we can allow $p$ to grow a bit faster than the usual ultra-high dimensional benchmark $e^{O(n)}$ rate. This is an interesting result and our conjecture is that this added advantage is due to considering future aggregation instead of just $k$-step ahead forecast for a fixed $k$. In other words, if prediction horizon $m$ is allowed to grow to $\infty$, the $m$-length average of residuals can automatically provide some concentration. Thus it can allow for scenarios where estimation of $\bm \beta$ is not very precise. We believe that this is an interesting exploration of the relaxation of the sparsity condition compared to the usual LASSO literature.
For the special case of the linear process (See (ref)), the conditions in ((ref)) will remain identical and one can provide some more specifications in the definition of $r$ in ((ref)) using the linear coefficients $a_i$ from ((ref)) in the view of the Nagaev-type concentration inequalities derived in Result (ref). Moreover, for the linear process, one can also state the corresponding results for the (LRDL) case, however since the condition on sparsity is unaffected by this nature of dependence, we do not state them separately here.
For the high-dimensional regression scenario, it is very natural to ask whether one could relax the fixed predictor settings to random design as it is practically impossible to find many perfectly predictable covariates. We show that under a mild condition on the possibly stochastic covariate process, the prediction intervals based on LASSO-optimized residuals also work in a stochastic design set-up. This is a particularly interesting extension from zhou10 since in the high-dimensional regime it is impractical to assume an exponentially growing number of covariates to be perfectly predictable for its future values. Moreover, apart from allowing a random design, we also allow the covariates to be temporally dependent. For this subsection, we restrict ourselves to only nonlinear processes for the sake of clarity. We assume
where $\epsilon_i^x$ are i.i.d., $\textbf{G}_x$ is a measurable function and $\bm x_t$ is $p \times 1$ high-dimensional mean-zero stationary process. Let $\{\epsilon_{i}^{x'}\}$ be an i.i.d. copy of $\{\epsilon_i^x\}$. Due to the high dimension of $\bm x_t$ and the dependence as a stochastic time-series, we opt for a dependence-adjusted uniform (over co-ordinates) functional dependence measure. For a fixed co-ordinate $1 \leq j \leq p$ and $q \geq 1$, define
Use this to define the dependence-adjusted norm as
The quantity $\Phi_{q,\alpha}$ provides a concise and natural measure of dependence which can effectively account for high dimensionality and temporal dependence. Let also the error process $(e_i)$ admit the following representation
where $\epsilon_i^x$ are i.i.d. and $G_e$ is a measurable function. One can then define the cumulative dependence using this functional dependence measure. However, thanks to the structure of the regression problem we can skip defining that for the $e_i$ process. Instead, for the quantile consistency in the case of stochastic design, we will need a notion of functional dependence on the cross-product process $x_{.j}e_.$ as follows
where $e_{k,l}^*=G_e(\epsilon_l^e,\epsilon_{i-1}^e,\ldots,\epsilon_{i-k}^e{'},\ldots)$ and $\{\epsilon_{i}^e\}'$ is an i.i.d. copy of $\{\epsilon_i^e\}$. Using ((ref)), we define the dependence adjusted norm as ((ref)) for the $x_{.j}e_.$ process uniformly over $1 \le j \le p$ as follows:
With this background on the stochastic covariate and the error process we now state the quantile consistency result below.
The prediction intervals based along the line of ((ref)) would need the future values of $\bm{x}_i$ which we do not observe. One possible solution to this is to fit a vector-autoregressive (VAR) model with appropriate lags and then estimate the $k$-step ahead predictions for $ 1 \leq k \leq m$ using the estimated matrix coefficients of the VAR process. It is also possible to lay down assumptions on the $\bm{x}_i$ and $e_i$ process and handle this in a much more rigorous way. But since our focus is on the relatively easier but asymptotically valid methods of estimating the prediction intervals based on quantiles, we omit that discussion. Instead, for practical implementation, (cf. our data analysis from Section (ref) where, for covariates, we used wind and temperature data that are stochastic in nature) we just use the past year-over-year mean as to substitute for $\bm{x}_i^{\mkern-1.5mu\mathsf{T}}\bm \beta $ for $i =n+1 ,\ldots,n+m$. Since $m$ is also growing, it is natural that the values of aggregated $\bm{x}_i^{\mkern-1.5mu\mathsf{T}}\bm \beta$ would concentrate around aggregated past year-over-year mean. We plan to discuss the issue of simultaneously estimating the future $\bm{x}_i$'s in a future work.
In this section, we compare the predictive performance of methods discussed above using OLS, LAD and LASSO estimator in both low-dimensional and high-dimensional setup. In the low-dimension setup the comparison between CLT based methods (Quenched CLT as described in Section (ref)) and QTL- Quantile based methods will be evaluated.
The focus here is on evaluation of PIs discussed in the previous section based on their coverage probability. We start by generating the error process $(e_t)$ as:
with $\epsilon_i$ i.i.d. from an $\alpha^*$-stable distribution. The heavy-tails index $\alpha^*=1.5$, autocovariance decay parameter $\gamma=-0.8$, speed-of-transition parameter $\delta=0.05$, autoregressive coefficients $\phi_1=0.6$ and $\phi_2=-0.3$, the noise standard deviation $\sigma=54.1$, and threshold $T=0$, were all selected based on the autoregressive models fitted to the electricity prices used later in the empirical part. The logistic transition function is given by $G(e_{i-1};\delta,T)=(1+\exp(-\delta(e_{i-1}-T)))^{-1}$. These three specifications represent
respectively. Eventually, we add a large number of exogenous covariates to the error process, obtaining $y_i= \bm{x}^{\mkern-1.5mu\mathsf{T}}_i\bm \beta +e_i, i=1,\ldots,n+m.$ We compute our PIs based on $(y_1,\bm x_1)\ldots,(y_n,\bm x_n)$ and evaluate them on $\bar{y}_{+1:m}=1/m\sum_{i=1}^my_{n+i}$. Note that we predict the averages instead of sums. This is motivated by easier comparison of predictive performance across different forecast horizons, and also turns out as more appropriate in the following empirical part.
Regarding covariates, we consider two scenarios (i) $p<n$ and (ii) $p>n$. In scenario (i) (See Table (ref)), we compare PIs based on OLS, LAD, and LASSO estimators. We set $n=8736$ ($\approx$ 1 year of hourly data), $m=168,336,504,672$ (1,2,3,4 weeks of hourly data), and $p=319$ (151 weather variables and 168 periodic variables), similarly to our empirical application (except that the horizon there spans up to 17 weeks) described in Section (ref). In scenario (ii) (See Table (ref)), we only use the LASSO as the other two estimators are not uniquely identified. We set\footnote{We reduce $n$ for computational convenience so that the $p>n$ does not have to be very large.} $n=336$ (2 weeks of hourly data), $m=24,48,72,96$ (1,2,3,4 days of hourly data) and $p=487$ (151 weather variables and 336 periodic variables).
The elements of $\bm \beta\in\mathbb{R}^p$ are i.i.d. from the uniform distribution\footnote{Previous version of this paper also contained results based on Cauchy distribution. Leading to the same general conclusions, we have omitted them for the sake of brevity.} $U[-1,1]$. Moreover, as properties of the LASSO estimator depend on the sparsity of $\bm \beta$, we assume $s=(1-\|\bm \beta\|_0/p)=50\%$. Supplementary online material contains results for sparsity $s=90\%$ and $20\%$, i.e., for high and low sparsity set-up. Throughout the experiment, we keep the (sparse) $\bm \beta$ fixed for all $1000$ repetitions. We compute PI for nominal coverage $(1-\alpha) = 60\%, 80\%, 90\%, 95\%$ and compare the quenched CLT method and QTL method proposed in Section (ref), based on their coverage probabilities (CP hereafter with the plural being CPs) \[ (\widehat{1-\alpha}) = \frac{1}{1000} \sum_{j=1}^{1000} \mathbb{I}\left([L,U]_{j,\hat{\bm \beta}}\ni\bar{y}_{j,+1:m}\right), \] and based on the average Winkler loss winkler72 \[ \mathcal{L} = |[L,U]_{j,\hat{\bm \beta}}| + \frac{2}{\alpha} \inf_{z\in [L,U]_{j,\hat{\bm \beta}}} |\bar{y}_{j,+1:m}-z|, \quad j=1,\ldots, 1000, \] where $\mathbb{I}$ for the $j$-th trial is 1 when $\bar{y}_{j,+1:m}$ is covered by the interval $[L,U]_{j,\hat{\bm \beta}}$ and 0 otherwise. The Winkler loss is a commonly used quantile-type loss, which penalizes the width of the PI and the size of misses, thus being appropriate for comparing PIs based on known quantiles askanazi18. Notably, the higher is the nominal coverage $1-\alpha$ the more weight is assigned to the misses via the factor $2/\alpha$.
Two general conclusions may be drawn from the experiment, independent of whether the dimensionality of covariates is high or low:
Further results require the distinction of the low-dimensional and high-dimensional scenarios. In the low dimensional scenario (see Table (ref)), in general, the following holds:
For the high-dimensional scenario, some of the previous statements do not hold. Similar to the design in chudyold, the horizon/sample ratios become very high, i.e., $m/n>1/4$. Moreover, we must face the curse of dimensionality for which we use LASSO\footnote{Details concerning the selection of tuning parameter $\lambda$ for LASSO can be found in Appendix C}. Still, this set-up has a largely negative impact on QTL's performance (see Table (ref)), and therefore we exploit a data-driven adjustment based on replication of the residual $\hat{e}_i=y_i-\hat{y}_i$ using stationary bootstrap. We denote the adjusted QTL by ADJ. Details concerning implementation are in Section (ref), where ADJ is used under similar a set-up as here. The simulations provide us with the following results:
Additionally, we may wonder what the impact of the sparsity and the generating distribution of $\bm{\beta}$ is. In general, the CPs are slightly higher when $\bm{\beta}$ is very sparse. In turn, $\bm{\beta}$ drawn from the Cauchy distribution leads to slightly smaller CPss. Still, both alternative set-ups lead to the same general conclusions with differences between CPss of identical methods (for identical series and horizons) within the range of two percentage points.
Next we compare and contrast forecasts obtained by our methods with the existing ones through POOS in a real-life data. Following is a list of competing methods we will explore:
We first give a quick overview of the last three methods above so that they are comparable to the ADJ as specified in Section (ref).
Robust Bayes PIs RBS:\\ For this sophisticated univariate approach, we focus on intuition and refer to the supplementary Appendix of mw16 for more details about the implementation. The robust Bayes PIs are specifically designed for long-horizon predictions, e.g., when $m/n\approx1/2$. First, the high-frequency noise is extracted out from $y_t$ using low-frequency cosine transformation. Projecting $\bar{y}_{+1:m}$ on the space spanned by the first $q$ frequencies is the key to obtaining the conditional distribution of $\bar{y}_{+1:m}$. In order to expand the class of processes for which this method can be used while keeping track of parameter uncertainty, mw16 employed a Bayesian approach. In addition, the resulting PIs are further enhanced to attain the frequentist coverage using the least favorable distribution. This requires advanced algorithmic search for quantiles of non-standard distributions, which is its main drawback in terms of implementation. On the other hand, their supporting online materials provide some pre-computed inputs which make the computation faster.
Bootstrap PIs for ARX, ETS and NAR:
We forecast $\bar{y}_{+1:m}=1/m\sum_{t=1}^m y_{n+t}$, i.e. the average of $m$ future hourly day-ahead spot electricity prices for Germany and Austria - the largest market in the European Power Exchange (EPEX SPOT). One of the reasons why we decided to forecast future averages was that the Bayes approach of mw16 is designed specifically for the means. Since all other methods are flexible, we used the means as a common basis for the comparison.
The prices arise from day-ahead hourly auctions where traders trade for specific hours of the next day. With the market operating 24 hours a day, we have $11 640$ observations between 01/01/2013 00:00:00 UTC\footnote{Coordinated Universal Time.} and 04/30/2014 23:00:00 UTC. We split the data into a training period spanning from 01/01/2013 00:00:00 UTC to 12/31/2013 23:00:00 UTC and an evaluation period spanning from 01/01/2014 00:00:00 UTC to 04/30/2014 23:00:00 UTC (see Figure (ref)A). The forecasting horizon is $m=1,\ldots,17$ weeks ($168,\ldots,2856$ hours).
Inspection of the periodogram for the prices in Figure (ref)C reveals peaks at periods 1 week, 1 day and $1/2$ day. The mixed seasonality is difficult to model by SARIMA or ETS models which are suitable for monthly and quarterly data or by dummy variables. Instead, we use sums of sinusoids with seasonal Fourier frequencies at $\omega_k=2\pi k/168$, $k=1,2,\ldots,\frac{168}{2}$ corresponding to periods 1 week, $1/2$ week, $\ldots$ , 2 hours bmrt07,wm08,cf05. The coefficients of linear combination $\bm \beta_k^{(s)}, \bm \beta_k^{(c)}$ can be estimated by least squares. In addition, we use 2 dummy variables as indicators for all weekends.
As mentioned in Section (ref), the local weather variables are also used as covariates. The weather conditions implicitly capture seasonal patterns longer than a week, which is very important when forecasting long horizons. Local weather is represented by 151 hourly wind speed, and temperature series is observed throughout 5 years (2009-2013), i.e., including the training period but not the evaluation period (see above). In order to approximate some missing in-sample data and unobserved values for the evaluation period, we take hourly-specific-averages\footnote{See alternative approximation of future values by bootstrap hf10} of each weather series over these 5 years. In total, we have 168 trigonometric covariates, 151 weather covariates and 2 dummies which gives a full set of 321 covariates.
Before we compare the ADJ with the other competitors, we will address a few issues that are usual with analysis of any real-life datasets. In Figure (ref)B, we see a drop in the price level during December 2013. The forecasts based on the whole training period would therefore suffer from bias. By contrast, using only the post-break December data would mean a loss of potentially valuable information. An optimal trade-off in such situations can be achieved by down-weighing older observations ppp13, also called exponentially weighted regression t10. In order to achieve better forecasting performance, we use the exponentially weighted regression with standardized exponential weights $v_{n-t+1}= \delta^{t-1}((1-\delta))/(1-\delta^t )$, $t=1,\ldots,n$ and with $\delta=0.8$. This applies to ADJ and NAR methods. The ETS and ARX models provide exponential down-weighing implicitly, but with optimally selected weights. mw16 showed that the RBS is robust to structural changes. We would also like to see if there are actual benefits from using disaggregated weather data instead of weather data aggregated across the weather stations. Therefore, we compute the ADJ PIs using no regressors as in Figure (ref)A, using only deterministic regressors as in Figure (ref)B, using deterministic regressors and aggregated weather variables defined as $\bar{w}_t=\sum_{k=1}^{73}w_{k,t}$, $\bar{\tau}_t=\sum_{l=1}^{78}\tau_{l,t}$ as in Figure (ref)C and finally, using all 321 covariates as in Figure (ref)D. As we can see, there is only very little difference between the first three plots, which means that using only deterministic regressors with or without the aggregated weather data does not prevent the bias at the end of the evaluation period. On the other hand, if we use the disentangled local weather data, significant improvement is achieved.
Finally, we get to the comparison with the alternative PIs denoted as RBS, ETS, NAR, and ARX. All these PIs are given in Figure (ref). Of the four methods, only RBS gives sensible PIs. RBS works consistently well over the whole $17$-weeks-long evaluation period (Figure (ref)A). However, when compared to the ADJ, the p.i's seem too conservative. Hence the ADJ provides more precision on top of decent coverage. Prediction intervals by ETS get too conservative as the horizon grows and do not provide a valid alternative to ADJ. The NAR is even more biased than the ADJ without covariates, especially for large $m$. Not so surprisingly, the ARX perform worst of all methods, presumably because the exponential down-weighing implied by the simple autoregression is too mild. Besides, the narrow PIs are the result of ignoring the parameter (among other types of) uncertainty.
We constructed quantile-based prediction intervals in a regression framework. From a theoretical perspective, we have extended the results of zhou10 to high-dimensional set-up and also to the case of the nonlinear error process. We showed the consistency of fitted residuals for the ultra-high dimensional case $\log p=o(n)$. Under some mild conditions on the fixed or stochastic covariate process, we were able to establish quantile consistency for the normalized average of fitted residuals and thus provide a significant extension to theoretical validity of the non-parametric and simple quantile-based prediction intervals. The quantile method has been additionally adjusted for short sample and long horizon and was successfully applied to predict spot electricity prices for Germany and Austria using a large set of local weather time series. The results have shown the superiority of the adjusted method over selected conventional methods and approaches such as exponential smoothing, neural networks as well as the recently proposed low-frequency approach of mw16.
Regarding future work, some interesting extensions can include multivariate target series and subsequent construction of simultaneous prediction intervals. Applications of such simultaneous intervals could include the prediction of spot electricity prices for each hour simultaneously in the spirit of rbd15.
We are thankful to Stefan Feuerriegel for providing data from their paper lfn15 for the analysis in a direct communication.
We are thankful to two anonymous referees and the editor for their helpful comments and corrections that helped improve the presentation of this paper significantly. The first and third authors are partially funded by NSF-DMS 2124222 NSF-DMS 1405410 respectively.