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.
89,019 characters · 21 sections · 77 citation commands
\centerline{\bf Abstract} \abst.arg\\ {KEY WORDS:} \key.arg
\pagestyle{plain}
\setcounter{equation}{0} \baselineskip=24pt Machine learning techniques have revolutionized data analysis across various domains, providing powerful tools for forecasting economic variables and financial outcomes gu2020empirical,babii2022machine,bia2024double. However, their application to interval-valued data remains largely unexplored, despite the increasing prevalence of such data in economic and financial spheres. interval-valued data, or more generally symbolic data, offer richer information (e.g., trends and volatility) than point-valued data within the same time period, and avoid noise contained in high-frequency data bock2000symbolic, billard2002symbolic, bock1999analysis, golan2017interval. In the data-rich environment, interval-valued data are prevalent across diverse situations, including interval-valued salary in the job advertisements zhong2023feature, high-low asset returns Gloria2013Constrained, wang2012cipca, a minimax regret portfolio selection problem GIOVE2006253, high-low livestock prices 2016Interval, Zhang2020hybrid, minimum-maximum daily air quality data YANG2019336, and variation of Cholesterol level wang2012linear. These interval-valued data, coupled with advanced machine learning techniques, offers an opportunity to enhance the analysis and forecasting of complex phenomena. As far as we know, there is little work on developing spares regression models for high-dimensional interval-valued data.
Our attempt in this article is to investigate the crucial problem of interval-valued machine learning methods in high-dimensional contexts, for which we immediately face several unique theoretical challenges in contrast to the case of point-valued time series data. First, interval-valued data introduce inherent complexities due to their distinct algebraic properties and operational rules, necessitating novel approaches beyond standard point-valued techniques. Second, different from quadratic loss with $L_1$ penalties for the existing literature for point-valued machine learning methods, we need to develop a proper loss function to simultaneously capture potential information contained in interval-valued data and yield sparse regression model. Third, the diverging dimensional setting in this paper substantially complicates the theoretical derivations, particularly in establishing the consistency and asymptotic normality of our estimators. A key reason for this is that classical large sample theories, such as the law of large numbers and central limit theorem, developed for point-valued data, cannot be directly applied to interval-valued data.
To address these challenges, we propose a sparse linear regression via machine learning tools to select relevant interval-valued features and optimize model parameters by treating intervals as inseparable sets. The estimation procedure is developed within a penalized minimum distance framework based on $D_K$-metrics with adaptive LASSO type penalty, and we introduce a novel interval-based least angle regression (ILARS) algorithm to solve the resulting optimization problem. Leveraging central limit theories based on $D_K$-distance, we establish asymptotic properties of the proposed parameter estimators. Specifically, under suitable regularity conditions, we prove that the penalized minimum distance estimators demonstrate consistency and asymptotic normality in both fixed and diverging dimensional settings, achieving the oracle properties\footnote{The oracle properties indicate that the estimators perform as well as if the true submodel were known a priori 2001Variable}. In addition, we introduce the interval-based framework of nonnegative garrote regression and ridge regression as additional regularized approaches for interval modeling. Furthermore, when interval-valued data reduces to point-valued data, the proposed loss function with a suitable kernel function is equivalent the quadratic loss. Simulation studies demonstrate the favorable finite sample properties of the proposed estimation. Empirical applications to crude oil prices forecasting and S&P 100 index-tracking highlight the merits of the proposed method compared to other competing approaches, including random forest and multilayer perceptron (MLP).
It is worth discussing some key references and outlining our contributions in relation to the most relevant literature. First, in contrast to existing interval regression methods Gloria2013Constrained,2016Interval,Lin2019extreme,LIMANETO20081500,LIMANETO2010333, our proposed approach treats intervals as inseparable entities to capture the potential information contained within the intervals. Previous methods, such as that of Gloria2013Constrained, construct interval models by modeling upper and lower bounds. While this approach facilitates the direct application of traditional point-data modeling techniques, including maximum likelihood and least squares estimation, it only utilizes the boundary information of intervals, thereby neglecting the potentially valuable information contained within the intervals. Although Gloria2013Constrained's method has been widely adopted, this limitation of considering only boundary points is also prevalent in subsequent studies DIAS20171118,buansing2020information. To overcome this methodological shortcoming, our proposed penalized $D_K$-distance measure utilizes the information of boundaries as well as the interior points of an interval. By utilizing the potential information contained in intervals, our approach is expected to achieve superior forecasting accuracy.
Second, the proposed parsimonious model with adaptive LASSO represents the first attempt to simultaneously achieve estimation and consistent variable selection for high-dimensional ITS. While parsimonious regression techniques for point-valued data have been extensively studied and advanced over the past decades, from the seminal least absolute shrinkage and selection operator tibshirani1996regression to various extensions including adaptive LASSO zou2006adaptive, group LASSO Meier2008grouplasso, with continuing developments in recent years bia2024double,gao2024robust,caner2024should. However, analogous methodologies for interval-valued data remain largely unexplored. While some machine learning approaches for interval-valued data exist, such as artificial neural networks YANG2019336, support vector machines utkin2019imprecise, and visualization techniques zhang2022visualization, none of these methods is designed to handle the dual challenges of high dimensionality and sparsity inherent in modern interval-valued datasets. To the best of our knowledge, only zhong2023feature explicitly addresses high-dimensional interval-valued data. Nevertheless, their approach is restricted to cases where only the response variable is interval-valued while predictors remain point-valued, making it unsuitable for our context where both predictors and response are interval-valued.
Third, we establish the consistency and oracle properties of the proposed penalized estimator, regardless of whether the number of predictors is diverging with the sample size. While several methods treating interval variables as inseparable entities exist, their theoretical frameworks are primarily restricted to low dimensions. For instance, han2016vector introduced an interval-valued vector autoregressive moving average (IVARMA) model, where each component can be viewed as an autoregressive conditional interval model with interval-valued exogenous variables (ACIX). To capture nonlinear features of ITS, sun2018threshold developed threshold autoregressive interval models, which was further extended by yang2024LSTIAR to a logistic smooth transition interval autoregressive model. Similar set-based approaches have been widely used in various fields, including stock market wang2016set, foreign exchange market sun2020assessing, and commodity market wu2023pass,he2021forecasting, among many others. However, none of these existing methodologies addresses the challenges posed by high-dimensional interval-valued data, which is the focus of our study.
The rest of this paper is organized as follows. Section 2 introduces the interval-based machine learning regression via adaptive LASSO. Section 3 describes the proposed estimators' asymptotic properties and introduces the interval-based nonnegative garrote and ridge regression. Section 4 develops the associated asymptotic properties of diverging-dimensional interval regression. Section 5 presents the simulation studies to show the finite sample properties of our method. Section 6 shows two empirical applications on the interval-valued crude oil price forecasting and the interval-based index tracking. Section 7 concludes the paper and discusses its future prospects. The appendix presents technical assumptions, implementation algorithms, and mathematical proofs, with additional simulation results in the supplementary materials.
We first review the ACIX models in the existing literature sun2018threshold,he2021forecasting. Suppose $\{Y_t\}$ and $\{X_{j,t}\}$ are stationary ITS\footnote{The interval variable is defined as a measurable map on a probability space $(\Omega,\mathcal{F},P)$, namely $Y: \Omega \rightarrow I_{\mathbb{R}}$, where $I_{\mathbb{R}}$ is the set of all pairs of ordered numbers in $\mathbb{R}$. Specifically, for any $w$ in $\Omega$, the term of $I_{\mathbb{R}}$ takes the form of $Y(w)=[Y_L(w),Y_R(w)]$, where $Y_R<Y_L$ is allowed. See more discussions in Remark (ref)}. Then, an ACIX model\footnote{For any given intervals $Y_1=[Y_{1L},Y_{1R}]$, $Y_2=[Y_{2L},Y_{2R}]$ and scalar $c$, the operation rules of intervals are defined as follows: (1) addition: $Y_1+Y_2=[Y_{1L}+Y_{2L},Y_{1R}+Y_{2R}]$; (2) Hukuhara's difference: $Y_1-Y_2=[Y_{1L}-Y_{2L},Y_{1R}-Y_{2R}]$; (3) scalar multiplication: $c\cdot Y_1=[c\cdot Y_{1L},c\cdot Y_{1R}]$.} can be expressed as
where ${\bf X}_t=(X_{1,t},...,X_{J,t})'$, $\alpha_0,\beta_j$ and ${\boldsymbol \delta}_j=(\delta_{j,1},\dots,\delta_{j,J})'$ are unknown scalar parameters, $I_0$ is the interval unit element $[-\frac{1}{2},\frac{1}{2}]$, and $u_t$ is an interval martingale difference sequence (IMDS), satisfying $\mathbb{E}[u_t|\mathcal{I}_{t-1}]=[0,0]$ with $\mathcal{I}_{t-1}$ being the information set. Let ${\boldsymbol \theta}=(\alpha_0,\beta_0,\beta_1,\dots,\beta_q,{\boldsymbol \delta}_0',\dots,{\boldsymbol \delta}_s')'=(\theta_1,\dots,\theta_{p})'$ and ${\bf Z}_t=([1,1],I_0,Y_{t-1},\dots,Y_{t-q},{\bf X}_t',\dots,{\bf X}_{t-s}')'$, where $p$ is the dimension of all parameters being either fixed or diverging as $T\rightarrow\infty$. More discussions of ACIX model\footnote{ Similar to the point-valued case, we consider the interval regression model (ref) to be correctly specified in conditional mean, if there exists a true parameter ${\boldsymbol \theta}^0 \in \mathbb{R}^p$ such that $\mathbb{E}[Y_t|{\bf Z}_t] = {\bf Z}_t'{\boldsymbol \theta}^0$. Otherwise, we can define the pseudo-true parameter ${\boldsymbol \theta}^*$ as ${\boldsymbol \theta}^*=\mbox{argmin} \mathbb{E}[\lVert Y_t-{\bf Z}_t'{\boldsymbol \theta}\rVert _K^2]$, where ${\bf Z}_t'{\boldsymbol \theta}^*$ represents the optimal linear combination in the sense of minimizing the $D_K$ distance. Moreover, when (ref) contains no exogenous variables, it reduces to the ACI model.} and the definition of intervals can be found in the existing literature, e.g., he2021forecasting, Yang2016, and among others.
The ACIX model serves as a generalization of the popular ARX-type model, commonly employed for point-valued time series analysis. It provides a framework to capture the temporal dependence of interval processes observed in economics and finance, such as volatility or range clustering and level effects. However, in the era of big data, the dimension of ${\bf Z}_t$ is often sufficiently large to encompass the underlying structure of high-dimensional interval-valued data, and there may exist sparsity within the predictors. For instance, when modeling the dynamics of interval-valued stock returns influenced by multiple factors guo2023statistical, such as supply, demand, geopolitical tensions, and technological advancements, or when analyzing interval-valued macroeconomic indicators using various predictors koop2023bayesian, like GDP growth or inflation forecasting ranges subject to uncertainties arising from consumer confidence, trade policies, and industrial productivity. Additionally, this challenge arises in modeling interval-valued exchange rates premanode2013improving, which are impacted by numerous factors, including economic fundamentals, market sentiment, and central bank interventions. In such scenarios, the ACIX model becomes less suitable due to its fixed number of parameters and inability to handle redundant variables.
To select important interval-valued predictors from a large dataset, we propose a penalized minimum distance estimation for the ACIX model. { Without loss of generality, we assume that the response and covariates are standardized.} Our objective is to estimate the unknown regression coefficients by solving the following penalized regression problem based on $D_K$-distance:
where $\lVert\cdot\rVert_K$ is the $D_K$ norm derived from the $D_K$-distance, $\lambda_T$ is the tuning parameter, $w_j$ is a known adaptive weight for $j=1,\dots,p$, and $T$ is the sample size. In practice, we can set the adaptive weight as $\widehat{w}_j=1/|\tilde{\theta}_j|^\gamma$, where $\tilde{{\boldsymbol \theta}}=(\tilde{\theta}_1,...,\tilde{\theta}_p)'$ is the minimum $D_K$-distance estimator\footnote{ The minimum $D_K$-distance estimator $\tilde{{\boldsymbol \theta}}=\mbox{argmin} \ensuremath{\sum_{t=1}^{T}}\lVert Y_t-{\bf Z}_t'{\boldsymbol \theta}\rVert _K^2$, and it has been proven to be a $\sqrt{T}$ consistent estimator with fixed dimension $p$.} han2016vector, and $\gamma$ is a given constant.
The $D_K$-distance is a metric for measuring the distance between two interval-valued variables korner2002variance,han2016vector,sun2018threshold,he2021forecasting. Specifically, the $D_K$-distance and its derived norm $\lVert\cdot\rVert_K$ are defined as:
where $K(u,v)$ is a symmetric positive definite kernel function for $u, v \in S^0=\{1,-1\}$, and $s_{Y_t}(u)$ is the support function of interval $Y_t$ defined on the unit sphere $S^0=\{-1,1\}$ as: \[ s_{Y_t}(u)=\left\{
\right. \]for $u\in S^0$. It is noteworthy that the operation rules and $D_K$ norm facilitate the construction of a complete normed linear space. Furthermore, we can endow this space with an inner product induced by the $D_K$-distance metric, denoted as $\langle\cdot,\cdot\rangle_K$.\footnote{ For example, suppose $X_1$ and $X_2$ are two intervals. The inner product of them is $\langle s_{X_1},s_{X_2} \rangle_K$. For simplicity of notation, we extend the use of $\langle\cdot,\cdot\rangle_K$ to also denote the multiplication of interval matrices, which can be obtained by replacing the pointed-valued multiplication with inner product for intervals.}
Our penalized minimum $D_K$-distance estimation covers several classical special cases of interval models. The $D_K$-distance is, to a certain degree, equivalent to the $d_W$ distance\footnote{ The $d_W$ distance for intervals is defined as $d_W(A,B)=\sqrt{\int_{[0,1]}(A(\omega)-B(\omega))^2dW(\omega)}$ for $A,B\in I_\mathbb{R}$, where $W(\omega)$ is a probability measure on the real Borel space $([0,1],{\bf B}([0,1]))$.} introduced by bertoluzza1995new, with the advantage of being more computationally tractable. The $d_W$ distance measure involves not only distances between extreme points with weights $W(0)$ and $W(1)$, but also distances between interior points in the intervals with weights $W(\omega),0<\omega<1$. It is interesting to see that the $D_K$ metric, as a equivalence of the $d_W$ metric, preserves this property, which is demonstrated through examples in the special cases. In the following, we investigate various special choices of kernel $K(u,v)$ and discuss the corresponding penalized regression. For notational convenience, we denote the kernel $K$ more concisely as: $K(1,1)=a$, $K(1,-1)=K(-1,1)=b$, and $K(-1,-1)=c$. Let $Y_{m,t}=(Y_{L,t}+Y_{R,t})/2$, $Y_{r,t}=Y_{R,t}-Y_{L,t}$, ${\bf Z}_{m,t}=({\bf Z}_{L,t}+{\bf Z}_{R,t})/2$, ${\bf Z}_{r,t}={\bf Z}_{R,t}-{\bf Z}_{L,t}$ be the midpoints and ranges of $Y_t$ and ${\bf Z}_t$. Then, denote $\Delta_{m,t}=Y_{m,t}-{\bf Z}_{m,t}'{\boldsymbol \theta}$, $\Delta_{r,t}=Y_{r,t}-{\bf Z}_{r,t}'{\boldsymbol \theta}$, $\Delta_{L,t}=Y_{L,t}-{\bf Z}_{L,t}'{\boldsymbol \theta}$, and $\Delta_{R,t}=Y_{R,t}-{\bf Z}_{R,t}'{\boldsymbol \theta}$.
{\bf Case 1.} $a=1/4, b=-1/4, c=1/4$. The $D_K$ norm becomes $||Y_t-{\bf Z}_t'{\boldsymbol{\beta}}||_K^2=D_K(Y_t,{\bf Z}_t'{\boldsymbol{\beta}})^2 = \Delta_{m,t}^2$.\footnote{In this case, the $D_K$-distance is equivalent to the $d_W$ distance, where $W(\omega)$ is a distribution such that $W(1/2)=1$ and $0$ otherwise.} Thus, the penalized minimum distance estimation is obtained by
which is the adaptive LASSO estimation for the midpoints of intervals. Under this kernel function, our interval model effectively utilizes only the midpoint information of the intervals. Especially, if there is no penalty in (ref), note that our method degenerates to the midpoints method proposed by billard2000regression.
{\bf Case 2.} $a=1, b=1, c=1$. The $D_K$ norm becomes $||Y_t-{\bf Z}_t'{\boldsymbol{\beta}}||_K^2=D_K(Y_t,{\bf Z}_t'{\boldsymbol{\beta}})^2 = \Delta_{r,t}^2$. This leads to the following optimization problem: \[ \widehat{{\boldsymbol \theta}}_T=\mbox{argmin}_{\boldsymbol \theta} \ensuremath{\sum_{t=1}^{T}} (Y_{r,t}-{\bf Z}_{r,t}'{\boldsymbol \theta})^2 + \lambda_T\sum_{j=1}^p \frac{1}{|\tilde{\theta}_j|^\gamma}\lvert \theta_j\rvert. \] In this case, our method only use the range information of the ITS. It is equivalent to the adaptive LASSO estimation for the ranges of intervals.
{\bf Case 3.} $a,c>0, b=0$. We have $||Y_t-{\bf Z}_t'{\boldsymbol{\beta}}||_K^2=D_K(Y_t,{\bf Z}_t'{\boldsymbol{\beta}})^2=a\Delta_{R,t}^2+c\Delta_{L,t}^2$.\footnote{If $a+c=1$, the choice of such a kernel $K$ is equivalent to the choice of $W(\omega)$ in $d_W$ distance with $W(\omega)$ follows a Bernoulli distribution with $W(0)=c$ and $W(1)=a$.} The penalized estimation (ref) becomes
In this case, the estimator $\widehat{{\boldsymbol \theta}}_T$ is obtained by minimizing the square errors of weighted bounds with a penalization. Moreover, when there is no penalty (i.e., $\lambda_T=0$), our method encompasses several popular special cases. First, note that (ref) is similar to the constrained Minmax method. The Minmax method estimates the lower and upper bounds of the intervals using different parameter vectors billard2002symbolic, thereby ignoring the dependence between the bounds. Then, (ref) can also be seen as the bivariate regression in brito2007modelling with a constraint. With this case, our model is essentially equivalent to using information about the interval's left and right bounds.
{\bf Case 4.} $a=c, |b|<a$. It follows that, $||Y_t-{\bf Z}_t'{\boldsymbol{\beta}}||_K^2 = D_K(Y_t,{\bf Z}_t'{\boldsymbol{\beta}})^2 = \frac{a+b}{2}\Delta_{r,t}^2+2(a-b)\Delta_{m,t}^2$.\footnote{If $a-b=1$ and $b\leq0$, the choice of such a kernel $K$ is equivalent to the choice of $W(\omega)$ in $d_W$ distance with $W(\omega)$ follows a distribution such that $W(0)=a+b$, $W(1/2)=-4b$, and $W(1)=c+b$.} Then, (ref) takes the form of
In this case, the estimator $\widehat{{\boldsymbol \theta}}_T$ is obtained by the penalized square errors of weighted ranges and midpoints. When $\lambda_T=0$, equation (ref) provides an approach similar to the well-known CRM method proposed by LIMANETO20081500, but with additional constraints. Importantly, our method under this kernel function is equivalent to utilizing both midpoint and range information of intervals.
{\bf Case 5.} $a\neq c, b\neq 0$. We have, $||Y_t-{\bf Z}_t'{\boldsymbol{\beta}}||_K^2 = D_K(Y_t,{\bf Z}_t'{\boldsymbol{\beta}})^2 = a\Delta_{R,t}^2 + c\Delta_{L,t}^2 -2b\Delta_{R,t}\Delta_{L,t}$ or equivalently $||Y_t-{\bf Z}_t'{\boldsymbol{\beta}}||_K^2=(a+2b+c)/4\Delta_{r,t}^2 + (a-2b+c)\Delta_{m,t}^2+(a-c)\Delta_{r,t}\Delta_{m,t}$.\footnote{if $a+c-2b=1$ and $b<0$, the $D_K$-distance is equivalent to $d_W$ distance with distributions as $W(0)=a+b$, $W(1/2)=-4b$, $W(1)=c+b$, and 0 otherwise.} In this case, $\widehat{{\boldsymbol \theta}}_T$ is obtained by solving the following optimization problem:
Here, (ref) could capture the information in $\Delta_{R,t}$, $\Delta_{L,t}$, and $\Delta_{R,t}\Delta_{L,t}$. Utilizing the cross product information will enhance estimation efficiency.
In this section, we examine the asymptotic properties of estimation under the condition that the dimension of the predictors in the penalized ACIX model is large but fixed. In the following, the $L_2$ norm of any vector is denoted by $||\cdot||$. Theorems (ref) and (ref) establish the consistency and asymptotic normality, respectively, of the penalized estimators for interval linear regression. The necessary conditions for these theorems are listed in Appendix A.
Theorem (ref) specifically states that under the imposed assumptions, the distance between the estimated parameter vector $\widehat{{\boldsymbol \theta}}_T$ and the true parameter vector ${\boldsymbol \theta}^0$ in the penalized interval regression converges in probability to zero at a rate of $O_p(1/\sqrt{T})$. This theorem provides a theoretical guarantee for the consistency of the estimator $\widehat{{\boldsymbol \theta}}_T$, ensuring that it converges to the true parameter values at a well-defined rate as the sample size $T$ increases.
In the following analysis, we establish two main asymptotic results: the consistency of interval-valued variable selection and the asymptotic normality of non-zero coefficient estimators. For exposition purposes, we assume without loss of generality that the first $k_0$ predictors are the true variables with non-zero coefficients, while the remaining $m_0=p-k_0$ predictors are redundant with zero coefficients. Let ${\boldsymbol \theta}^0=({\boldsymbol \theta}_1^{0'},{\boldsymbol \theta}_2^{0'})'$, where ${\boldsymbol \theta}_1^0$ is a $k_0\times 1$ vector of non-zero coefficients and ${\boldsymbol \theta}^0_2$ is a $m_0\times1$ vector of zero coefficients, i.e., ${\boldsymbol \theta}^0_1\neq\0$ and ${\boldsymbol \theta}^0_2=\0$. Note that this partition is solely for theoretical analysis, as the true non-zero and zero coefficients are unknown in practice. Let $\widehat{{\boldsymbol \theta}}_T=(\widehat{{\boldsymbol \theta}}_{1T}',\widehat{{\boldsymbol \theta}}_{2T}')'$ denote the estimator corresponding to ${\boldsymbol \theta}^0_1$ and ${\boldsymbol \theta}^0_2$, respectively. Further, let ${\bf Z}_{1t} = (Z_{1t},...,Z_{k_0t})'$ represent the vector of the first $k_0$ covariates.
Theorem (ref) states two asymptotic properties of the penalized interval regression. The first property is the consistency in variable selection, meaning that the probability of correctly identifying the zero coefficients approaches 1 as the sample size increases. The other property is that estimated coefficients of active variables follow an asymptotic normal distribution with mean zero and variance, in terms of $C_{11}^{-1} \mathbb{E}[\langle s_{{\bf Z}_{1t}},s_{u_t}\rangle_K \langle s_{u_t},s'_{{\bf Z}_{1t}}\rangle_K] C_{11}^{-1}$.
Following the spirit of breiman1995better, we propose the interval-based nonnegative garrote, a machine learning technique for adaptive feature selection. We also analyze the relationship between penalized interval regression via adaptive LASSO and via nonnegative garrote. The nonnegative garrote for intervals is equivalent to minimizing the following loss function,
where $c_j$ is a constant, $\lambda_T$ is the tuning parameter, $\widehat{\theta}^j_{garrote}= c_j\tilde{\theta}_j$ is the interval-based nonnegative garrote estimator, and $\tilde{\theta}_j$ denotes the minimum $D_K$-distance estimator as defined in Section 2.1. Based on equation (ref) and Theorem (ref), we can derive the following corollaries.
If the conditions in Corollary (ref) are satisfied, i.e., $\gamma$ in estimation (ref) is set to be 1 and $\widehat{w}_j=1/|\tilde{\theta}_j|$, the penalized estimation for interval regression takes the form:
Let $\theta_j=c_j \tilde{\theta}_j$ and $\theta_j/\tilde{\theta}_j\geq 0$ (or $\theta_j \tilde{\theta}_j\geq 0$). Then, (ref) is equivalent to (ref). Corollary (ref) shows the consistency of variable selection of the interval-based nonnegative garrote.
Next, we propose a ridge regression method for interval-valued data, introducing a regularized approach that directly incorporates interval structures. The interval-based ridge regression is given by
Since the optimization (ref) is a quadratic problem, the following corollary presents a closed-form solution for parameter estimation derived through direct differentiation.
In the preceding sections, we have analyzed the asymptotic properties of the interval linear regression model via adaptive LASSO when the dimension of regressors is fixed. Besides, the number of regressors could be diverging, which has been studied for point-valued data, like fan2004nonconcave and huang2008asymptotic. In this section, we focus on a diverging dimension of regressors as the sample size increases, that is, $p$ grows to infinity at some slower rates than the sample size $T$. We first show the consistency of the minimum $D_K$-distance estimator for ACIX model when the dimension of predictors is diverging.
Theorem (ref) is a generalization of estimation consistency of the conventional ACIX model. It can also be seen as a generalization of the ordinary least square (OLS) estimation consistency with diverging dimension for the interval-valued case.
Without loss of generality, we assume that the coefficients of the first $k_T$ variables are non-zero. It follows that the coefficients of the last $m_T=p-k_T$ variables are zeros. Moreover, we still let ${\bf Z}_{1t}$ denote the first $k_T$ covariates and ${\bf \Sigma}_{1T} = \frac{1}{T}\ensuremath{\sum_{t=1}^{T}}\langle s_{{\bf Z}_{1t}},s_{{\bf Z}_{1t}'}\rangle_K$. Before giving the consistency and oracle properties of the penalized estimation with diverging dimension, we first show the following Lemma.
Note that $\rho_{1T}$ appears in the denominators of $h_T$ and $h_T'$, which makes it possible that $h_T'$ may converge to zero faster than $h_T$ if $\rho_{1T}\rightarrow0$. Additionally, if we suppose that there exists a positive constant $\rho_1$ and $\rho_2$ such that $0<\rho_1<\rho_{1T}<\rho_2<\infty$, Theorem (ref) yields that the convergence rate of $h_T=O_p(\sqrt{p/T})$ and $h_T\leq h_T'$. Thus we have $||\widehat{{\boldsymbol \theta}}_T-{\boldsymbol \theta}^0||=O_p(\sqrt{p/T})$. This condition is also common in the existing literature, such as Condition (F) in fan2004nonconcave. Furthermore, if $p$ is finite, Theorem (ref) degenerates to Theorem (ref).
Theorem (ref) (i) indicates that the penalized linear regression for interval-valued data is consistent in variable selection, that is, the estimators of the zero coefficients are exactly zero with high probability when $T$ is large. Moreover, Theorem (ref) (ii) states that the estimators of the nonzero parameters have the asymptotic normal distribution when the number of parameters diverges. Similar results in the point-valued case can be found in existing literature, such as fan2004nonconcave, huang2008asymptotic and huang2008adaptive. While these studies considered the independent and identically distributed point-valued random variables, this paper proves the oracle properties for martingale difference ITS. Furthermore, Theorem (ref) demonstrates that our proposed penalized model can effectively identify the true non-zero coefficients while shrinking irrelevant coefficients to zero, thus achieving model sparsity. This variable selection property is particularly valuable in high-dimensional interval-valued settings, where it not only enhances model interpretability but also potentially improves prediction accuracy by reducing model complexity.
This section investigates the finite sample performance of the proposed penalized estimation for interval regression. The two-stage minimum $D_K$-distance estimators are considered as estimated adaptive weights of the penalized estimation he2021forecasting, han2016vector, and the tuning parameter $\lambda_T$ is selected by a five-fold cross-validation process.
First, we consider the data generating process (DGP) as follows:
where $Y_t$, $X_{j,t}$, and $u_t$ are all interval variables, ${\boldsymbol \theta}=(\alpha_0,\beta_0,\delta_1,\dots,\delta_{p-2})'$ is the given point-valued coefficients. The ordered pairs $X_{j,t}=(X_{L,j,t},X_{R,j,t}), j=1,...,p-2,$ are generated from the bivariate normal distributions with non-zero covariance matrix. To generate the interval innovations $\{u_t\}$, we employ an ACI(1,0) process following sun2018threshold:
where the parameters $(\alpha_0, \beta_0, \beta_1)'$ are estimated by a two-stage minimum $D_K$-distance method; $Y_t=\ln(P_t)-\ln(P_{t-1})$ and $P_t$ are designated as the time series of the daily S&P 500 index data for the period January 2, 2015 to December 31, 2019, with its bounds as the high and low prices in day $t$. From model (ref), we have the estimated interval innovation, namely $\widehat{u}_t=Y_t-(\widehat{\alpha}_0+\widehat{\beta}_0 I_0+\widehat{\beta}_1 Y_{t-1})$. We then generate $\{u_t\}$ via the naive bootstrapping from $\{\widehat{u}_t\}$, with sample size $T$. In addition, following han2016vector and sun2018threshold, a two-stage minimum $D_K$-distance method is also used here with some given preliminary kernels $K$ in the first stage. We take a preliminary kernel in the first step with $K=(a,b,c)=(5,1,1)$.\footnote{ han2012autoregressive proved that the two-stage $D_K$-distance estimator is asymptotically most efficient among all symmetric positive definite kernels satisfying $K(1,1)>0$, $K(1,1)K(-1,-1)>K(1,-1)^2$ and $K(1,-1)=K(-1,1)$. Thus, as sample size $T$ increases to infinity, the choice of kernel $K$ in the first stage has little impact on the optimal kernel derived by the two-stage estimation.}
In our simulation experiments, we consider two DGPs:
{\bf DGP 1:} The dimension of predictors is fixed at $p=10$. Following zou2006adaptive, the initial coefficients are set as ${\boldsymbol \theta}=(\alpha_0,\beta_0,\delta_1,\dots,\delta_8)'=(0,0,3,1.5,0,0,2,0,0,0)'$. We examine the estimators under sample sizes $T=20, 40, 80$.
{\bf DGP 2:} The dimension of predictors diverges as the sample size increases. We set $p=[3T^{1/3}]$, where $[\cdot]$ represents the largest integer not exceeding $3T^{1/3}$. The initial coefficients are ${\boldsymbol \theta}=(0,0,11/4,-23/6,37/12,-13/9,1/3,0,...,0)'$. We consider sample sizes $T=100,200,400,800$.
Each experiment is repeated 1000 times. To evaluate the performance of our estimation method, we employ bias (Bias), standard deviation (SD), and root mean square error (RMSE) for each estimated parameter $\widehat{\theta}_j$, i.e.,
where $N=1000$ is the number of replications, $\theta_j$ is the true parameter value, and $\bar{\theta}_j=\frac{1}{N}\sum_{i=1}^N\widehat{\theta}_j^{(i)}$ is the average of the estimated $\widehat{\theta}_j$ across all replications.
The LARS algorithm has been used to compute the solution path for the LASSO problem in point-valued cases Efron2004Least,Tibshirani2013The. For our interval-valued case, we propose an interval-based LARS algorithm, as detailed in Appendix B, to obtain the estimators. We set the parameter $\gamma$ to 0.5 and 1 in our experiments.
In this section, the interval innovation $(u_{L,t},u_{R,t})$ is generated from a bivariate normal distribution, $(u_{L,t},u_{R,t}) \sim i.i.d. N(0,{\bf \Sigma}^0)$. The DGP for this case can be expressed as:
where $\alpha_0$, $\beta_0$, and $\delta_j$ are given initial scalar parameters, and ${\bf \Sigma}^0$ is a $2\times 2$ positive definite matrix with diagonal elements equal to 1 and all other elements equal to 0.75. The regressors $(X_{L,j,t},X_{R,j,t})$ $(j=1,\ldots,p-2; t=1,\ldots,T)$ are also generated from bivariate normal distributions with non-zero covariance matrices.
In this section, we also consider two types of DGPs: (1) where $p$ is fixed, and (2) where $p$ diverges as the sample size increases. We define these DGPs as follows:
{\bf DGP 3:} $Y_t=[Y_{L,t},Y_{R,t}]$ and $u_t=[u_{L,t},u_{R,t}]$ are generated according to (ref). We set $p=10$, ${\boldsymbol \theta}=(\alpha_0,\beta_0,\delta_1,\dots,\delta_8)'=(0,0,3,1.5,0,0,2,0,0,0)'$, and $T=20, 40, 80$.
{\bf DGP 4:} $Y_t=[Y_{L,t},Y_{R,t}]$ and $u_t=[u_{L,t},u_{R,t}]$ are generated according to (ref). We set $p=[3T^{1/3}]$, ${\boldsymbol \theta}=(0,0,11/4,-23/6,37/12,-13/9,1/3,0,\ldots,0)'$, and $T=100,200,400,800$.
All other parameter settings remain the same as those in Section (ref). We employ the criteria defined in equations (ref) - (ref) to evaluate the performance of our estimators.
Panel A of Table (ref) reports the Bias, SD and RMSE of our method and the minimum $D_K$-distance estimators of ACIX model based on DGP 1, with $\gamma=0.5$. Several observations emerge from this panel. First, the Bias, SD and RMSE of each estimator whose true value is set to be zero are approaching zero as the sample size $T$ increases, consistent with asymptotic efficiency of variable selection in Theorem (ref). For example, RMSE of $\alpha_0$ decreases from $3.1521 \times 10^{-3}$ to $0.4728\times 10^{-3}$ as $T$ increases from 20 to 80. Second, for the nonvanishing coefficients, the evaluation criteria of estimators of $\delta_1,\delta_2$ and $\delta_5$ converge to zero as $T$ increases. These observations indicate the consistency of the estimated nonvanishing parameters and provide a finite sample evidence for the oracle properties. For example, SD of $\delta_1$ decreases from $0.4867\times 10^{-3}$ to $0.1588\times 10^{-3}$ as $T$ increases from 20 to 80. Third, our method yields substantially improved estimates compared to the minimum $D_K$-distance method. As evidenced in Panel A, our approach demonstrates superior performance relative to ACIX, producing Bias, SD, and RMSE values that more closely approximate zero. This underscores the enhanced efficacy of our penalized method over the minimum $D_K$-distance estimation in scenarios characterized by model sparsity. For example, when the sample size $T=80$, the Bias for $\delta_3$ of our method is $0.0111\times 10^{-3}$, which is smaller than $0.0354\times 10^{-3}$ of the minimum $D_K$-distance method. When $u_t$ follows a bivariate normal distribution, similar results can be obtained from Panel B of Table (ref).
Table (ref) shows the evaluation results of DGP 2 and DGP 4 with $K=(5,1,1)$ and $\gamma=1$. In these cases, our proposed penalized minimum distance estimation still outperforms the minimum $D_K$-distance estimation of ACIX model. In Table (ref), most results of Bias from our estimation provide values closer to zero than those of minimum $D_K$-distance estimation. Nearly all results of SD and RMSE from our estimation are smaller than those of the minimum $D_K$-distance estimation. For example, when $T=400$ and $u_t$ is generated by an ACI process, the Bias, SD, and RMSE of $\delta_1$ are 0.0074, 0.2311, and 0.2312 (all $\times10^{-3}$), which are smaller than those of the benchmark estimation: 0.0079, 0.2312, and 0.2537 (all $\times10^{-3}$). Furthermore, Table (ref) shows that our model makes more accurate estimation of zero coefficients. For example, when $T=400$ and $u_t$ is generated by a bivariate normal distribution, the evaluation results of $\delta_6$ are -0.0003, 0.0109, and 0.0109, which are 70.0%, 36.6%, and 36.6% better than those of benchmark estimation. A possible explanation for this is that our estimation process could shrink the estimators of the zero coefficients to zero. The evaluation results of other parameters are listed in the online appendices of this paper.
Accurate crude oil price forecasting is an important yet controversial issue in economic and management research. Numerous studies have demonstrated that crude oil prices are influenced by a myriad of financial and macroeconomic factors, such as supply and demand dynamics, stock market performance, interest rates, exchange rates, monetary policy, and other commodity prices NASER201675, WEI2017141, he2021forecasting, he2010empirical. Importantly, many of these factors are represented in the form of interval-valued data, reflecting the inherent uncertainty and variability in their measurement. However, most existing work has focused on point-valued data, potentially overlooking valuable information contained within the interval-valued representations. Modeling interval-valued data may capture more comprehensive information, thereby enhancing the accuracy of crude oil price predictions. Consequently, we are motivated to identify and select important interval-valued factors from various potential variables via shrinkage methods, with the expectation of improving the forecasting accuracy of interval-valued crude oil prices.
This section describes the data used in our analysis, focusing on the monthly interval-valued West Texas Intermediate (WTI) crude oil futures prices. These prices are constructed from daily closing prices sourced from the New York Mercantile Exchange (NYMEX), a major marketplace for crude oil futures trading. Denote the ITS of crude oil prices as $Y_t=[Y_{L,t},Y_{R,t}]$, where $Y_{L,t}$ and $Y_{R,t}$ are constructed by taking the logarithm of the minimum and maximum daily closing prices within month $t$, respectively. The sample period spans from January 2006 to December 2019. Figure (ref) illustrates the bounds and range of these interval-valued prices over time.
We draw some interesting observations from Figure (ref). First, the ITS of crude oil price captures intra-month variations that monthly closing prices fail to reflect. Additionally, there appears to be a strong correlation between the lower and upper bounds of the interval. Second, the range of oil prices tends to increase as the price level decreases, suggesting that crude oil prices may become more volatile during downward trends. This pattern demonstrates that price level and volatility are two distinct aspects of crude oil price movements, likely exhibiting a negative correlation. Third, two considerable drops are evident in the trend of interval-valued prices. The first occurs from July to December 2008, attributable to the subprime crisis, while the second spans from June 2014 to January 2016, resulting from shale oil shocks.
As a crucial strategic resource on a global scale, crude oil typically exhibits high volatility in its pricing, influenced by a myriad of factors. Table (ref) presents several factors used in our application, including stock market indicators, monetary market variables, and supply and demand metrics. These factors have been widely used in existing literature wang2016forecasting,chai2018forecasting,Yang2016. To assess the stationarity of our data, we applied the Augmented Dickey-Fuller (ADF) test to $Y_{L,t}$, $Y_{R,t}$, and the bounds of explanatory variables. The results indicate that all point-valued series are first-order stationary.
To explore the performance of sparse model for ITS, we first employ the ACIX model han2016vector with all predictors as a benchmark forecasting model. Next, following the spirit of Gloria2013Constrained and SUN2021MA, the center-range method (CRM) and the constrained center-range method (CCRM) are considered as benchmark forecast methods. These two methods first proposed by LIMANETO20081500,LIMANETO2010333 can be expressed as: $y_t^m = \beta_0^m+\beta_1^m x_{1,t}^m + \cdots +\beta_p^m x_{p,t}^m + \epsilon_t^m$ and $y_t^r = \beta_0^r+\beta_1^r x_{1,t}^r + \cdots +\beta_p^r x_{p,t}^r + \epsilon_t^r$, where $\{y_t^m, x_{i,t}^m\}$ are midpoints of interval observations, and $\{y_t^r, x_{i,t}^r\}$ represent the ranges. Both CRM and CCRM estimators can be obtained from two separate point-valued least squares estimation methods. CCRM imposes restrictions on the coefficients of range, i.e., $\beta_i^r\geq 0$ for $1\leq i\leq p$, to ensure the range variables are nonnegative. A bivariate model for lower and upper bounds (BLU) is also employed as a benchmark forecast model. BLU model is based on estimation in the following system: $y_t^l = \beta_0^l+\beta_1^l x_{1,t}^l + \cdots +\beta_p^l x_{p,t}^l + \epsilon_t^l$ and $y_t^u = \beta_o^u+\beta_1^u x_{1,t}^u + \cdots +\beta_p^u x_{p,t}^u + \epsilon_t^u$. {Furthermore, to evaluate the performance of our proposed interval-based adaptive LASSO model relative to alternative machine learning approaches, we utilize random forest and MLP models as benchmarks for interval-valued data, denoted as IRF and IMLP\footnote{The random forest model is constructed with 100 decision trees using the TreeBagger algorithm. The MLP uses a single hidden layer with 10 neurons, trained for maximum 1000 epochs with error goal of $10^{-5}$ and minimum gradient of $10^{-6}$, with a 70:15:15 data split ratio for training, validation and testing.}. They employ a strategy of estimating the lower and upper bounds of intervals separately.}
The evaluation criteria listed in Table (ref) have frequently been used in interval regression models maciel2017evolving, Rodrigues2015, maciel2023adaptive, yang2024LSTIAR to evaluate the forecast accuracy of these models. The criteria in Panel A measure the gap between predicted and actual intervals, while the other criteria in Panel B measure forecast accuracy of special points in predicted intervals. Specifically, $\omega_1$ measures the nonoverlapping area of the forecasting and actual intervals, and $\omega_{D_K}$ measures the $D_K$-distance between $\widehat{Y}_t$ and $Y_t$. $\omega_{NSD1}$ and $\omega_{NSD2}$ evaluate the nonoverlapping area of $\widehat{Y}_t$ and $Y_t$ relative to their union set, where $\omega(\cdot)$ denotes the width and $R(\cdot)$ the range of an interval, see more details in sun2018threshold. Moreover, $\omega_{MDE}$ is about the mean distance error. In addition, $\omega_{rate}$ refers to the non-efficiency rate. For the point-based criteria, all four are the special cases of RMSE.
We employ a rolling window approach to study the out-of-sample performance of our estimation and other interval-based methods. A rolling estimation scheme is adopted for $60$ months with the first estimation sample spanning from January 2006 to December 2010. We conduct the one-step-ahead out-of-sample forecasts from January 2011 to December 2019, including $T_f=108$ forecasting periods. In addition, we also adopt a rolling estimation scheme for $120$ months, namely from January 2006 to December 2015. The last observations are forecasting sample. For each fixed rolling window, we compare our method's performance with that of the benchmark methods.
Table (ref) shows the out-of-sample performance of seven interval forecasting methods: ACIX, CRM, CCRM, BLU, IRF, IMLP, and our proposed method, focusing on interval-based criteria as outlined in Panel A of Table (ref). Among the seven interval-based methods considered, our method consistently ranks highest in forecasting performance across all interval-based criteria. For instance, with the training sample size of 60, our model's $D_K$-distance measure $w_{D_K}$ is 0.0086, outperforming all the other values. This superior performance can be attributed to two main factors. First, our interval model treats the interval oil price sample as an inseparable set and employs information of distances between both boundaries and interior points. Particularly, the correlation between the interval's lower and upper bounds is taken into account. As Figure (ref) shows, because the bounds of interval crude oil price are not independent, modeling the interval's lower and upper bounds separately does not fully utilize the information of the intervals. Second, the penalized interval regression provides a sparse model with fewer explanatory variables, effectively excluding the influence of unrelated and weakly related variables. Compared with our interval-based machine learning method, the underperformance of classic machine learning algorithms, i.e., IRF and IMLP, in predicting interval-valued oil prices can be attributed to two main factors. First, they simply apply regression to the two endpoints of the interval, neglecting the internal information. Second, the relatively small sample size could lead to overfitting in these complex models, particularly if not adequately tuned for the specific dataset.
The point-based criteria are outlined in Table (ref). It is observed that our method consistently outperforms the other six benchmark interval methods across all point-based criteria. To verify these findings, a Diebold-Mariano (DM) test is conducted on the out-of-sample forecast performance of various points within the forecasting intervals (including lower and upper boundaries, midpoint, and radius). The results of the DM tests, denoted by asterisks in Table (ref), provide strong evidence of the superiority of our PLR method. For example, in Panel A, we observe that our sparse method significantly outperforms all benchmark models at the 1% level for nearly all criteria. The only exception is the radius measure for the IRF model, where the difference is significant at the 5% level. Panel B results remain significant, albeit less so than Panel A. This difference likely stems from a smaller test set sample in Panel B, reducing statistical power.
Index tracking is a crucial technique in portfolio management, enabling investors to replicate the performance of a benchmark index while minimizing tracking error. Traditionally, most existing literature on index tracking have primarily relied on closing prices as the input data source corielli2006factor,wu2014nonnegative,strub2018optimal. However, these approaches may fail to capture the full extent of price fluctuations that occur throughout the trading day. Neglecting intraday price movements can potentially lead to suboptimal portfolio construction and increased tracking error, as the closing price alone may not accurately reflect the true dynamics of the underlying assets, which is a critical aspect in index tracking. To address this issue, incorporating interval-valued data can provide a more comprehensive representation of asset price movements, which can better account for volatility and price variations, potentially leading to improved tracking performance and reduced tracking error. Consequently, we are motivated to develop interval-based index tracking methodologies, aiming to construct portfolios that more accurately replicate the benchmark index's performance. To the best of our knowledge, this paper is the first to propose a replicating strategy based on interval-valued stock prices, marking a novel contribution to the field of index tracking.
In this application, we develop and evaluate an interval-based strategy for tracking the S&P 100 index. The interval-valued log return of stocks is constructed as $[r_{l,t},r_{h,t}]$, where $r_{l,t}=\ln\frac{P_{low,t}}{P_{close,t-1}}$ and $r_{h,t}=\ln\frac{P_{high,t}}{P_{close,t-1}}$. This representation of daily interval-valued returns, also adopted by Gloria2013Constrained, captures investors' high and low expectations and provides more comprehensive trading information. As shown in Figure (ref), the S&P 100 index returns demonstrated notable volatility, particularly during the COVID-19 pandemic period.
Our index tracking strategy comprises two main steps. First, we employ the proposed penalized linear interval regression to select stocks from the S&P 100 constituents, with the number of selected stocks controlled by the tuning parameter. Second, we estimate the weights of the chosen stocks using OLS regression on closing prices.\footnote{Our strategy does not enforce full investment (sum of weights needs not equal one) and allows short selling, enabling direct OLS estimation of weights. See shu2020high for details about this approach.} For comparison, we construct a point-based benchmark strategy that follows the same two-step procedure but uses LASSO tibshirani1996regression for stock selection, with closing prices utilized in both steps.
To implement and evaluate these strategies, we use a rolling window approach with a 250-day training period (approximately one trading year) and a 21-day testing period (approximately one trading month). The process involves selecting 10 stocks from the S&P 100 constituents using both interval-based and point-based methods, followed by weight estimation through OLS regression. To ensure robustness, we examine three different training samples beginning from the first trading day of 2017, 2018, and 2019, respectively, comparing both in-sample and out-of-sample performance. All data are sourced from Wind and Yahoo Finance.
To evaluate the performance of our interval-based strategy and the point-based strategy, we employ the following two widely used measures of tracking error wu2014nonnegative,corielli2006factor:
where $T$ is the sample size, $err_t = r_t - \widehat{r}_t$, $\overline{err}$ is the sample mean of $err_t$, and $r_t$ and $\widehat{r}_t$ are the returns from the tracking portfolio and the index, respectively.
Figure (ref) displays the tracking error and mean absolute deviation of both the interval-based and point-based models. To offer more detailed insights, we illustrate the cumulative error in Figure (ref), where the lines represent $S(\tau)$ and $M(\tau)$ as $\tau$ ranges from 1 to $250$ in the training sample and from 1 to 21 in the test sample.\footnote{ To avoid confusion, it is necessary to clarify that the model is estimated from the entire set of 250 training samples, and we only vary $\tau$ during evaluation.}
Several observations can be obtained from Figure (ref). First, the interval-based strategy outperforms the point-based one in-sample, as measured by both types of tracking errors. At the beginning, the blue line occasionally exceeds the red line, which is attributed to the instability resulting from the limited sample size. As the sample size used for calculating cumulative error gradually increases, the in-sample performance of our method notably outperforms that of the point-based method. Second, the out-of-sample tracking errors in 2018 and 2019 of the interval-based method are mostly smaller than those of the point-based method. The outperformance of our proposed interval-based index tracking strategy illustrates that interval-valued data may contribute to improving conventional portfolio strategies. Third, when the out-of-sample is 2020, neither method demonstrates a significant advantage over the other. One possible reason is the extreme volatility in the stock market caused by the COVID-19 pandemic in early 2020, making it challenging for past performance to capture market behavior during such extreme shocks. Overall, the application of the proposed interval variable selection methods in index tracking demonstrates that interval-valued data can improve portfolio strategies. This inspires us to develop deeper research of interval models in financial studies in the future work.
In this paper, we propose a sparse regression for high-dimensional interval-valued data via machine learning tools. It is shown that the proposed method enjoys the oracle properties, i.e., the consistency in variable selection and the asymptotic normality of estimators. We further extend the proposed method to a diverging-dimensional interval case. Additionally, we also propose an interval-based LARS algorithm to solve the solution path of the estimation. Furthermore, simulation studies confirm the asymptotic properties of our method. Empirical applications highlight that our method improves the out-of-sample forecasts of crude oil price and works well in index tracking.
The proposed machine learning technique for interval linear regression is an interval-valued extension of the adaptive LASSO, originally designed for point-valued data. Several important avenues for future research emerge from this work. One potential extension is to adopt more sophisticated and powerful machine learning techniques, such as neural networks YANG2019336 and random forests, to handle interval-valued data based on random set theory, rather than simply applying existing tools to model the interval bounds. Moreover, other dimension reduction techniques for interval-valued data can be proposed, including new principal component analysis based on $D_K$-distance and interval factor models He2022Large,wang2012cipca.
{ \baselineskip=16pt
}
\setcounter{table}{0} \setcounter{equation}{0} \setcounter{remark}{0}
{