Generalized estimation equations method for fixed effects panel interval-valued data models

1.   Introduction

In the real world, there are a large number of interval-valued data. For example, in the stock market, the K-line chart, as an important decision-making tool, includes the opening price, closing price, maximum, and minimum of the day. In temperature prediction, the minimum and maximum temperatures can be regarded as an interval. On the other hand, due to the complexity of the system, the data obtained is an interval rather than a real point. For example, experts often predict the economic growth rate of the next year in the form of an interval, such as 4–5%. As far as we know, interval-valued data can provide a tool for dealing with large datasets, which were proposed by Diday ^[1]. It not only contains richer information than point data but also has better interpretability. So, it is meaningful to analyze interval-valued data. Billard and Diday ^[2] proposed the first algorithm to fit the interval-valued linear regression, which fits the linear regression model to the midpoint of the interval values and obtains the parameters by minimizing the midpoint error. Billard and Diday ^[3] proposed the min-max method, which defined two models, each corresponding to a response bound. The lower bounds of the response depend on the lower bounds of the regressive variables, and the upper bounds of the response depend on the upper bounds of the regressive variables. Lima Neto et al. ^[4] introduced the center and range Method (CRM) and proposed two linear models: one for the midpoint and another for the range. Kong and Gao ^[5] proposed a regularized interval MM estimate (RIMME) for interval-valued regression, which can achieve a good balance between the prediction accuracy and mathematical coherence of the predicted intervals. Xu and Qin ^[6] proposed a bivariate Bayesian regression model based on CRM with known or unknown covariance matrices. Beyaztas et al. ^[7] introduced the functional forms of some well-known regression models that take interval-valued data. Interval analysis technologies can also be found in other literature ^{[8,9,10,11,12]}. To sum up, it is not difficult to find that the above research on interval-valued data is mainly based on cross-sectional data.

Panel data or longitudinal data is obtained by multiple observations of individuals at different times. Statistical models that combine cross-section and time series real value data are becoming more and more popular in economic research, and panel data sets have certain advantages over traditional pure cross-section or pure time series data sets. Galvao Jr. ^[13] studied a quantile regression dynamic panel model with fixed effects. Aristodemou ^[14] studied semiparametric identification in linear index discrete response panel data models with fixed effects. Hamiye Beyaztas and Bandyopadhyay ^[15] proposed a new, weighted likelihood-based robust estimation procedure for linear panel data models with fixed and random effects. Bai et al. ^[16] studied generalized least squares (GLS) estimation for linear panel data models. Baltagi and Li ^[17] focused on prediction in spatial models based on panel data. Elhorst ^[18] studied the specification and estimation of spatial panel data models. Zhang et al. ^[19] studied a penalized quantile regression for spatial panel models with fixed effects. In particular, in recent years, a few scholars have begun to study interval panel data models. For example, Ji et al. ^[20] introduced fixed effects panel data regression models for interval-valued data, and presented three kinds of panel interval-valued data regression models and the estimation of their parameters. Zhang et al. ^[21] proposed a robust estimation method based on the iterative weighted least squares technique to reduce the impact of outliers on models. Considering the correlation between the center and range or the upper and lower bounds of intervals, Ji et al. ^[22] proposed the bivariate maximum likelihood (BML) method for estimating unknown parameters of the model. Usually, the observation data of different individuals are independent of each other, but there is correlation between the observation data of the same individual at different times. As far as we know, current research on panel interval-valued data has not considered the correlation within a group of the data.

Generalized estimation equations (GEEs) proposed by Liang and Zeger ^[23] are commonly used to analyze the longitudinal data with correlation within individuals. Crowder ^[24] studied the (weak) consistency and inconsistency of the solutions of general estimating equations. Balan and Schiopu-Kratina ^[25] studied the consistency and asymptotic normality of the GEEs estimation, with the covariate dimension fixed. Wang ^[26] studied the consistency and asymptotic normality of GEEs estimates when the sample size of the covariate dimension tends to infinity. Wang and Carey ^[27] showed that the discrepancy between the working correlation structure and the true correlation structure and the estimation method of the correlation coefficient affect the asymptotic relative efficiency. Pan ^[28] proposed a new model selection criterion that minimizes the expected predictive bias of estimating equations. Wang ^[29] presented a systematic review on GEEs, covering foundational concepts as well as several recent developments due to practical challenges in real applications. Seaman and Copas^[30] described doubly robust (DR) GEEs, and illustrated their use on simulated data. To sum up, previous studies have only applied GEEs to point value data. Therefore, we intend to apply it to panel interval-valued data in this paper.

In this paper, unlike the independent hypothesis of other studies, we consider the correlation within a group of the panel interval-valued data to obtain better parameter estimates. GEEs are commonly used to analyze the longitudinal data with correlation within individuals. We intend to apply it to panel interval-valued data in this paper. The results of Monte Carlo simulations reveal that our proposed method outperforms the LSDV method based on independent assumptions in terms of estimation accuracy and robustness. When the working correlation matrix is assumed correctly, the GEEs method performs better than the LSDV method. When the working correlation matrix is assumed incorrectly, the performance of the two methods is similar. In addition, the higher the degree of correlation within the group, the better the performance of the GEEs method. The empirical results also demonstrate the superiority of our proposed method.

The rest of this paper are organized as follows: Current regression models for fixed effects panel interval-valued data and the LSDV method are provided in Section 2. The regression method based on the generalized estimation equations is presented in Section 3. In Section 4, some Monte Carlo simulations and real data analysis show that, in contrast with the least-squares dummy-variable (LSDV) method, the proposed GEEs method has advantages in forecasting performance. Finally, Section 5 gives the conclusion and discussion.

2.   Fixed effects panel interval-valued data models

In this section, the LSDV method and three existing kinds of fixed effects panel interval-valued data models are introduced: the min-max (Min-Max or MM) model, the center and range (CRM) model, the center (CM) model. These three kinds models hold their own advantages and disadvantages.

For the panel interval-valued dataset $Z = \{(X_{it}, y_{it})\, |\, i = 1, 2, ..., N;t = 1, 2, ..., T\}$, let $X_{it} = (x_{it1}, x_{it2}, ..., x_{itp})$, where $x_{itk} = [x_{itk}^l, x_{itk}^u], y_{it} = [y_{it}^l, y_{it}^u], x_{itk}^c = \frac{x_{itk}^l+x_{itk}^u}{2}, x_{itk}^r = \frac{x_{itk}^u-x_{itk}^l}{2}, y_{it}^c = \frac{y_{it}^l+y_{it}^u}{2}, y_{it}^r = \frac{y_{it}^u-y_{it}^l}{2}, k = 1, 2, ..., p$. $x_{itk}^l$($x_{itk}^u, x_{itk}^c$ or $x_{itk}^r$) are considered as explanatory variables and $y_{it}^l$($y_{it}^u, y_{it}^c$ or $y_{it}^r$) are considered as response variables. In addition, let $u_i^l(u_i^u)$ be the individual effect and $u_i^c = \frac{u_i^u+u_i^l}{2}, u_i^r = \frac{u_i^u-u_i^l}{2}$.

2.1.   The min-max model

The min-max model of fixed effects panel interval-valued data regression consists of two different models. Two regression models are established for the upper and lower bounds of the response variable by using the upper and lower bounds of the explanatory variables. The specific form is given as follows:

where the errors $\epsilon_{it}^l\sim\ N(0, \sigma^2)$, $\epsilon_{it}^u\sim\ N(0, \sigma^2), i = 1, 2, ..., N; t = 1, 2..., T.$

Based on the lower bounds of $y_{it}$ and $x_{itj}\, (i = 1, 2..., N; t = 1, 2, ..., T; j = 1, 2, ..., p)$, the parameters $u_i^l$ and $\beta^l = (\beta_1^l, \beta_2^l, ..., \beta_p^l)^T$ are estimated by the least-squares dummy-variable (LSDV) method:

where $A^l = \sum_{i = 1}^{N}\sum_{t = 1}^{T}(X_{it}^l-\overline{X}_i^l)(X_{it}^l-\overline{X}_i^l)^T, B^l = \sum_{i = 1}^{N}\sum_{t = 1}^{T}(X_{it}^l-\overline{X}_i^l)(y_{it}^l-\overline{y}_i^l)^T, X_{it}^l = (x_{it1}^l, x_{it2}^l, ..., x_{itp}^l)^T, \overline{X}_i^l = \frac{1}{T}\sum_{t = 1}^{T}X_{it}^l, \overline{y}_i^l = \frac{1}{T} \sum_{t = 1}^{T} y_{it}^l, \overline{x}_{i.j}^l = \frac{1}{T} \sum_{t = 1}^{T}x_{itj}^l.$

The estimates of $\beta^u$ and $u_i^u$ have the following similar conclusion:

where $A^u = \sum_{i = 1}^{N}\sum_{t = 1}^{T}(X_{it}^u-\overline{X}_i^u)(X_{it}^u-\overline{X}_i^u)^T, B^u = \sum_{i = 1}^{N}\sum_{t = 1}^{T}(X_{it}^u-\overline{X}_i^u)(y_{it}^u-\overline{y}_i^u)^T, X_{it}^u = (x_{it1}^u, x_{it2}^u, ..., x_{itp}^u)^T, \overline{X}_i^u = \frac{1}{T}\sum_{t = 1}^{T}X_{it}^u, \overline{y}_i^u = \frac{1}{T} \sum_{t = 1}^{T} y_{it}^u, \overline{x}_{i.j}^u = \frac{1}{T} \sum_{t = 1}^{T}x_{itj}^u.$

For a new observation $X_{it} = (x_{it1}, x_{it2}, ..., x_{itp})^T$, with $x_{itj} = [x_{itj}^l, x_{itj}^u]$, the prediction of $\hat{y}_{it} = [\hat{y}_{it}^l, \hat{y}_{it}^u]$ is given by:

Sometimes the min-max model does not guarantee the mathematical coherence of the predicted interval bounds, and then the response variables are defined as follows:

2.2.   The center and range model

The center and range model of fixed effects panel interval-valued data regression is as follows:

where the errors $\epsilon_{it}^c\sim\ N(0, \sigma^2)$, $\epsilon_{it}^r\sim\ N(0, \sigma^2), i = 1, 2, ..., N; t = 1, 2..., T.$ The estimates of $\beta^c$ and $u_i^c$ given by the LSDV method are:

where $A^c = \sum_{i = 1}^{N}\sum_{t = 1}^{T}(X_{it}^c-\overline{X}_i^c)(X_{it}^c-\overline{X}_i^c)^T, B^c = \sum_{i = 1}^{N}\sum_{t = 1}^{T}(X_{it}^c-\overline{X}_i^c)(y_{it}^c-\overline{y}_i^c)^T, X_{it}^c = (x_{it1}^c, x_{it2}^c, ..., x_{itp}^c)^T, \overline{X}_i^c = \frac{1}{T}\sum_{t = 1}^{T}X_{it}^c, \overline{y}_i^c = \frac{1}{T} \sum_{t = 1}^{T} y_{it}^c, \overline{x}_{i.j}^c = \frac{1}{T} \sum_{t = 1}^{T}x_{itj}^c.$

Similarly, we can obtain the estimates of $\beta^r$ and $u_i^r$ :

where $A^r = \sum_{i = 1}^{N}\sum_{t = 1}^{T}(X_{it}^r-\overline{X}_i^r)(X_{it}^r-\overline{X}_i^r)^T, B^r = \sum_{i = 1}^{N}\sum_{t = 1}^{T}(X_{it}^r-\overline{X}_i^r)(y_{it}^r-\overline{y}_i^r)^T, X_{it}^r = (x_{it1}^r, x_{it2}^r, ..., x_{itp}^r)^T, \overline{X}_i^r = \frac{1}{T}\sum_{t = 1}^{T}X_{it}^r, \overline{y}_i^r = \frac{1}{T} \sum_{t = 1}^{T} y_{it}^r, \overline{x}_{i.j}^r = \frac{1}{T} \sum_{t = 1}^{T}x_{itj}^r.$

For a new observation $X_{it} = (x_{it1}, x_{it2}, ..., x_{itp})^T$, with $x_{itj} = [x_{itj}^l, x_{itj}^u]$, the predictions of $\hat{y}_{it}^c$ and $\hat{y}_{it}^r$ are given by:

Then, the predictions of $\hat{y}_{it}^l$ and $\hat{y}_{it}^u$ are as follows:

2.3.   The center model

Based on the center value of the panel interval-valued data, the liner regression model with fixed individual effects is as follows:

where the errors $\epsilon_{it}^c\sim\ N(0, \sigma^2), i = 1, 2, ..., N; t = 1, 2, ..., T$.

As described earlier, the estimates of $\beta^c$ and $u_i^c$ are given by the LSDV method:

For a new observation $X_{it} = (x_{it1}, x_{it2}, ..., x_{itp})^T$, with $x_{itj} = [x_{itj}^l, x_{itj}^u]$, the prediction of $\hat{y}_{it} = [\hat{y}_{it}^l, \hat{y}_{it}^u]$ is given by:

3.   Fixed effects panel interval-valued data models based on the GEEs method

This section introduces generalized estimation equations(GEEs) for fixed effects panel interval-valued regression models that consider the data correlation within individuals.

For panel interval-valued data set $Z = \{(X_{it}, y_{it})\, |\, i = 1, 2, ..., N;t = 1, 2, ..., T\}$, $X_{it} = (x_{it1}, x_{it2}, ..., x_{itp})^T$ is assumed as $p \times 1$ interval-valued explanatory variables with $x_{itj} = [x_{itj}^l, x_{itj}^u]$, $y_{it} = [y_{it}^l, y_{it}^u]$ as the observed interval-valued dependent variables, $i = 1, 2, ..., N; t = 1, 2, ..., T; j = 1, 2, ..., p$. $u_i$ is an unobserved effect for individual $i$ and is time invariant. Let $Y_i^a = (y_{i1}^a, y_{i2}^a, ..., y_{iT}^a)^T, X_i^a = (X_{i1}^a, X_{i2}^a, ..., X_{iT}^a)^T$, where $a = l, u, c,$ or $r$. Let $Y^a = ((Y_1^a)^T, (Y_2^a)^T, ..., (Y_N^a)^T)^T, X^a = ((X_1^a)^T, (X_2^a)^T, ..., (X_N^a)^T)^T, u^a = (u_1^a, u_2^a, ..., u_N^a)^T$.

3.1.   The GEEs method of the min-max model

The matrix form of the min-max model of fixed effects panel interval-valued data regression is as follows:

To save space, we only present the estimation method of the model (3.1). Considering the correlation between the observations of the same individual at different times, and the independence between the observations of different individuals, we suppose that $\epsilon^l\sim N(0, \Sigma^l), cov(\epsilon_{ik}^l, \epsilon_{jl}^l) = 0 (i\neq j)$, where $\epsilon^l = (\epsilon_1^l, \epsilon_2^l, ..., \epsilon_N^l)^T, \epsilon_i^l = (\epsilon_{i1}^l, \epsilon_{i2}^l, ..., \epsilon_{iT}^l)$,

$\Sigma^l = \left(\begin{array}{cccc} \Sigma_1^l & 0 & \cdots & 0 \\ 0 & \Sigma_2^l & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \Sigma_N^l \end{array}\right),$ $\Sigma_i^l = \left(\begin{array}{cccc} var(\epsilon_{i1}^l) & cov(\epsilon_{i1}^l, \epsilon_{i2}^l) & \cdots & cov(\epsilon_{i1}^l, \epsilon_{iT}^l) \\ cov(\epsilon_{i2}^l, \epsilon_{i1}^l) & var(\epsilon_{i2}^l) & \cdots & cov(\epsilon_{i2}^l, \epsilon_{iT}^l) \\ \vdots & \vdots & \ddots & \vdots \\ cov(\epsilon_{iT}^l, \epsilon_{i1}^l) & cov(\epsilon_{iT}^l, \epsilon_{i2}^l) & \cdots & var(\epsilon_{iT}^l) \end{array}\right),$

and similarly, $\epsilon^u\sim N(0, \Sigma^u)$.

There is a problem for us to find the estimation of $\beta^l$ in the model (3.1) because of the fixed effects $u^l$. Specifically, simultaneously estimating $u^l$ and $\beta^l$ directly may incur "curse of dimensionality" since $u^l$ is high dimensional when $N$ becomes too large. First, similar to Section 2.1, we use the projection techniques to eliminate fixed effects in the model (3.1). Specifically, let $Q = I_{NT}-P$, $P = I_N\otimes(J_T/T), J_T = \tau_T \tau_T'$, where $I_N$ is an $N \times N$ identity matrix and $\tau_T$ is an $N \times 1$ vector of ones. Multiplying by $Q$ on the left for model (3.1), we have

where $\widetilde{Y}^l = QY^l, \widetilde{X}^l = QX^l, \widetilde{\epsilon}^l = Q\epsilon^l$, $Q(I_N\otimes \tau_T) = 0$, and $\widetilde{\epsilon}^l\sim N(0, Q\Sigma^l Q^T), Q = diag(Q_1, Q_2, ..., Q_N), Q_i = I_T-J_T/T$,

Then, the generalized estimating equations (GEEs^[23]) can be defined as follows:

where $\widetilde{Y}^l_i = Q_iY_i^l = (\widetilde{y}^l_{i1}, \dots, \widetilde{y}^l_{iT})^T$, $\widetilde{X}^{l}_i = Q_iX_i^l$, $V_i^l$ is the covariance matrix of $\widetilde{Y}^l_i$. According to Liang and Zeger (^[23]), we simplify the covariance of the $i$th subject $V_i$ by taking $V_i^l = (A_i^l)^{1/2}R_i(\alpha)(A_i^l)^{1/2}$, where $A_i^l = {\rm diag}\{var(\widetilde{Y}^l_{i1}), \dots, var(\widetilde{Y}^l_{iT})\}$, and $R_i(\alpha)$ is a common working correlation with nuisance parameter $\alpha$. In general, $V_i^l$ is unknown, and we often use the empirical estimators $\hat{V}_i^l$ for the $V_i^l$. Specifically,

where $\hat{\alpha}$ can be obtained by following the method proposed by Liang and Zeger (^[23]) in Section 3.4.1. $A_i^l(\hat{\beta^l})$ is the empirical estimator of $A_i$ based on the current estimator of $\beta^l$.

Thus, we implement the iterative algorithm to solve Eq (3.4). Given the current estimate $\hat{\alpha}$, we use the following iterative procedure for $\beta^l$:

where $\hat{V}_{i}^l = (A_i^l)^{1/2}(\hat{\beta^l}^{(k)})R_i(\hat{\alpha})(A_i^l)^{1/2}(\hat{\beta^l}^{(k)})$.

Then, the estimates of the fixed effects $u^l$ are given as follows:

where $\overline{Y}^l = PY^l$, $\overline{X}^l = PX^l$.

Using similar methods, we can obtain the estimates of $\hat{\beta}^u$ by iterating the following ralationship to convergence:

where $\hat{V}_{i}^u = (A_i^u)^{1/2}(\hat{\beta^u}^{(k)})R_i(\hat{\alpha})(A_i^u)^{1/2}(\hat{\beta^u}^{(k)})$, and the estimates of $\hat{u}^u$ are given by

where $\overline{Y}^u = PY^u$, $\overline{X}^u = PX^u$.

For a new observation $X = [X^l, X^u]$, the prediction $\hat{Y} = [\hat{Y}^l, \hat{Y}^u]$ is given by

3.2.   The GEEs method of the center and range model

The matrix form of the center and range model of fixed effects panel interval-valued data regression is as follows:

where $\epsilon^c\sim N(0, \Sigma^c), \epsilon^r\sim N(0, \Sigma^r), \Sigma^c = \frac{\Sigma^l+\Sigma^u}{4} = \Sigma^r$.

Similar to Section 3.1, we only present the estimation method of model (3.9). To eliminate the fixed effects, let $\widetilde{Y}^c = QY^c, \widetilde{X}^c = QX^c, \widetilde{\epsilon}^c = Q\epsilon^c$, and $Q(I_N\otimes \tau_T) = 0$. Multiplying (3.9) by $Q$ on the left gives the following intra-group regression model equation:

where $\widetilde{\epsilon}^c\sim N(0, Q\Sigma^c Q^T).$ Then, we can obtain the estimates of $\hat{\beta}^c$ by iterating the following ralationship to convergence:

where $\hat{V}_{i}^c = (A_i^c)^{1/2}(\hat{\beta^c}^{(k)})R_i(\hat{\alpha})(A_i^c)^{1/2}(\hat{\beta^c}^{(k)})$, and the estimates of $\hat{u}^c$ are given as follows:

where $\overline{Y}^c = PY^c$, $\overline{X}^c = PX^c$.

Similarly, estimates of $\beta^r$ can be obtained as follows:

where $\hat{V}_{i}^r = (A_i^r)^{1/2}(\hat{\beta^r}^{(k)})R_i(\hat{\alpha})(A_i^r)^{1/2}(\hat{\beta^r}^{(k)})$, and the estimates of $\hat{u}^r$ are given by

where $\overline{Y}^r = PY^r$, $\overline{X}^r = PX^r$.

For a new observation $X = [X^l, X^u]$, the prediction $\hat{Y} = [\hat{Y}^l, \hat{Y}^u]$ is given as follows:

where

3.3.   The GEEs method of the center model

The matrix form of the center model of fixed effects panel interval-valued data regression is as follows:

where $\epsilon^c\sim N(0, \Sigma^c), \Sigma^c = \frac{\Sigma^l+\Sigma^u}{4}$.

As described earlier, the estimates of $\beta^c$ and $u^c$ given by GEEs are as follows:

where $\hat{V}_{i}^c = (A_i^c)^{1/2}(\hat{\beta^c}^{(k)})R_i(\hat{\alpha})(A_i^c)^{1/2}(\hat{\beta^c}^{(k)})$, and the estimates of $\hat{u}^c$ are obtained as follows:

where $\overline{Y}^c = PY^c$, $\overline{X}^c = PX^c$.

For a new observation $X = [X^l, X^u]$, the prediction $\hat{Y} = [\hat{Y}^l, \hat{Y}^u]$ is given by

3.4.   Computation algorithm

3.4.1.   The working correlation matrix

The following introduces the structures of several commonly used working correlation matrices ^[23,31,32]:

● Independence: For the independence working correlation matrix, it is assumed that all the response variables are mutually independent. $\alpha = 0, R_i(\alpha) = I_T.$

● Exchangeable: Suppose that the correlation between any two observed values of the same individual are the same. The form of the correlation matrix is as follows:

where $\alpha$ is the unknown correlation coefficient,

where $r_{it} = \widetilde{y}_{it}-\widetilde{X}_{it}\beta$.

● AR(1): Suppose that the correlation depends on the time. The form of the correlation matrix is as follows:

where

3.4.2.   Algorithm

In this subsection, we present a feasible algorithm as follows. Consider the general form of the model: $Y = X\beta+(I_N\otimes \tau_T)u+\epsilon$.

Step 1: Use the matrix $Q$ to eliminate fixed effects $u$. Then the initial value $\hat{\beta}^{(0)}$ is obtained from the independent assumption.

Step 2: The estimates $\hat{\alpha}^{(k)}$ and $A_i^{1/2}(\hat{\beta}^{(k)})$ are obtained based on the assumed structure of the correlation matrix as in Section 3.4.1. Then $\hat{V}_{i}^{(k)} = A_i^{1/2}(\hat{\beta}^{(k)})R_i(\hat{\alpha})A_i^{1/2}(\hat{\beta}^{(k)})$ is obtained. Update $\hat{\beta}^{(k+1)} = \hat{\beta}^{(k)}+\left(\sum_{i = 1}^{N} (\widetilde{X}_{i})^{T} (\hat{V_{i}}^{(k)})^{-1} \widetilde{X}_{i}\right)^{-1}\left(\sum_{i = 1}^{N} (\widetilde{X}_{i})^{T} (\hat{V_{i}}^{(k)})^{-1}({\widetilde{Y}_i}-\widetilde{X}_i\hat{\beta}^{(k)})\right)$.

Step 3: Iterate Step 2 until convergence, and denote the final estimators as the GEEs estimators. Then $\hat{u} = \left((I_N\otimes \tau_T)^{T}(I_N\otimes \tau_T)\right)^{-1}(I_N\otimes \tau_T)^{T}(\overline{Y}-\overline{X}\hat{\beta})$.

Here we can also utilize the existing algorithm to obtain estimates of $\beta$, which can be obtained by using the geepack function package in the statistical software R language.

4.   Numerical illustrations

This section compares the proposed GEEs method with the LSDV method proposed by Ji et al.^[20] through computational experiments demonstration. In computational experiments, two types of panel interval-valued datasets are considered. The first is a dataset of synthetic panel interval-valued data, while the second is the air pollution dataset of the 26 cities in China's Yangtze River Delta.

4.1.   Monte Carlo simulations

Without loss of generality, we provide two configuration schemes for synthetic panel interval-valued datasets: one from the min-max model and another derived from the CRM model. In the simulation framework, 20 or 50 individuals are considered, and each individual consists of three attributes of explanatory variables over 14 continuous periods, that is, $N = 20$ or $50$ and $T = 14$. Let the first 10 periods of data be the train set, and the rest be the test set. Considering the actual problem scenarios, let the correlation matrix be the first-order autoregressive (AR(1)) correlation setting $R_1(\alpha)$ or exchangeable correlation setting $R_2(\alpha)$.

4.1.1.   Dataset configuration based on the CRM model

In this section, we consider the following panel interval-valued data-generating processes.

where $\beta^c = (\beta_1^c, \beta_2^c, \beta_3^c) = (1, 3, 2), \, \beta^r = (\beta_1^r, \beta_2^r, \beta_3^r) = (1, 2, 1)$; $x_{ijt}^c \sim U[5, 10]$, $x_{ijt}^r \sim U[1, 3]$; fixed effects $u_i^c$ and $u_i^r \sim U[0, 0.5]$; $\epsilon_i^c = (\epsilon_{i1}^c, \epsilon_{i2}^c, ..., \epsilon_{iT}^c) \sim N(0, \sigma_c^2 R_i(\alpha))$, $\epsilon_i^r = (\epsilon_{i1}^r, \epsilon_{i2}^r, ..., \epsilon_{iT}^r) \sim N(0, \sigma_r^2 R_i(\alpha))$, where $\sigma_c^2 = \sigma_r^2 = 4, i = 1, 2$; and let the correlation coefficient $\alpha = 0.9, 0.5, 0.3$ to generate data with different degrees of correlation.

4.1.2.   Dataset configuration based on the min-max model

In this section, we consider the following panel interval-valued data-generating processes.

where $\beta^u = (\beta_1^u, \beta_2^u, \beta_3^u) = (3, 2, 2), \beta^l = (\beta_1^l, \beta_2^l, \beta_3^l) = (2, 1, 1)$; $x_{ijt}^u \sim U[5, 10]$, $x_{ijt}^l \sim U[1, 3]$; fixed effects $u_i^u$ and $u_i^l \sim U[0, 0.5]$; $\epsilon_i^u = (\epsilon_{i1}^u, \epsilon_{i2}^u, ..., \epsilon_{iT}^u) \sim N(0, \sigma_u^2 R_i(\alpha))$, $\epsilon_i^l = (\epsilon_{i1}^l, \epsilon_{i2}^l, ..., \epsilon_{iT}^l) \sim N(0, \sigma_l^2 R_i(\alpha))$, where $\sigma_u^2 = \sigma_l^2 = 4, i = 1, 2$; and let the correlation coefficient $\alpha = 0.9, 0.5, 0.3$ to generate data with different degrees of correlation.

4.1.3.   Implementation and evaluation

Four criteria for assessing the performances of different models are shown as follows:

(1) Mean average absolute error (MAE), mean magnitude of relative error (MMER), and root mean square error (RMSE) defined by Ji et al.^[20] are:

(2) The rate of different intervals (RI) defined by Hu and He ^[33] is:

where $w(.)$ is the interval width.

Tables 1-4 show the fitting results when the working matrix is specified correctly, while Tables 5-8 show the fitting results when the working matrix is specified incorrectly. Figures 1–4 display the box plots of four evaluation indices calculated from the GEEs method in different models with the dataset generated based on the CRM model or the min-max model, when the working matrices are specified correctly. The results of the box plots are similar, when the working matrices are specified incorrectly.

● Tables 1 and 2 list the fitting results of $\beta$, and we have the following findings: (i) Regardless of the data generated by which model, the GEEs method performs better than the LSDV method. (ii) As $\alpha$ or $N$ increases, the mean square error of GEEs estimation decreases. It means that the higher the degree of the correlation within the group, the better the performance of the GEEs. This also verifies the consistency of the estimate.

● Tables 3 and 4 show that (i) the GEEs method performs better than the LSDV method, regardless of the degree of correlation within the group of the data and the sample size, according to the four indicators. (ii) The CRM model performs best when data is generated based on centers and the radius. Similarly, the min-max model performs best when data is generated based on minimum and maximum values.

● From Tables 5–8, we can see that the fitting result of the GEEs method is close to that of the LSDV method when the working correlation matrix is incorrectly specified.

● Figures 1 and 2 display the box plots of four evaluation indices calculated from the GEEs method in different models with the dataset generated based on the CRM model, when the working matrices are specified correctly with $N = 20$, $\alpha = 0.9, 0.3$. We can see that the CRM model performs better than the other models.

● Figures 3 and 4 display the box plots of four evaluation indices calculated from the GEEs method in different models with the dataset generated based on the min-max model, when the working matrices are specified correctly with $N = 20$, $\alpha = 0.9, 0.3$. We can see that the min-max model performs better than the other models.

4.2.   Real data analysis

As everyone knows, $PM_{2.5}$ is the main pollutant affecting human health and atmospheric environment quality. With the continuous development of social economy and the acceleration of the industrialization process, the problem of air pollution is also increasing. Studies have shown that $O_3, SO_2, NO_2, CO, PM_{10}$, five air pollutants, have a significant impact on $PM_{2.5}$ concentrations. Exploring the relationship between these air pollutants and $PM_{2.5}$ can provide help for effective prevention and control of air pollution. So, this section selects the 26 cities in China's Yangtze River Delta as the research object, and analyzes the interval-valued data of air pollution for them. The daily data of six air quality indices of the Yangtze River Delta urban agglomeration during the 59 days from January 1 to February 28, 2023, can be obtained from the China Air Quality Online Monitoring platform (https://www.aqistudy.cn/historydata/). Through inspection, there are no outliers or missing values in the data. Then we consider the maximum and minimum values of the daily data of each index as the interval-valued dataset.

The observed data of $PM_{2.5}$ are taken as response variables, and the other five indices are taken as independent variables. The data from the first 54 days are used as the training set, and the remaining 5 days are used as the test set. Substituting the data into the models is shown in Sections 2 and 3. Finally, we demonstrate the fitting performance of the GEEs method and the LSDV method when the working correlation matrix is AR(1). The predictive performance of the two methods are listed in Table 9, and the regression results of the GEEs method are shown in Table 10.

● Table 9 shows that both methods have good fitting effects on actual data in different models. Moreover, we can see that the GEEs method had better fitting performance than the LSDV method in different models. Taking into account the four performance indices, the min-max model performs better than other models, which may be due to the type of experimental data, which is also consistent with the conclusions in the Monte Carlo experiments.

● Table 10 shows that there is a positive correlation between $PM_{10}, NO_2, CO$, and $PM_{2.5}$, while there is a negative correlation between $SO_2, O_3$, and $PM_{2.5}$. Especially, $CO$ has the greatest impact on $PM_{2.5}$.

5.   Conclusions and discussion

In this paper, we proposed an innovative approach which was to apply GEEs to fixed effects panel interval-valued data models, in which the correlation within a group is considered. Monte Carlo simulation experiments test the feasibility of this method and our proposed method was applied to $PM2.5$ forecasting. Experiments show that the proposed method has a better performance than the LSDV method proposed before.

However, the estimation method we proposed has some limitations. For example, the model can be extended to the nonlinear regression model to capture complex nonlinear relationships between variables. In addition, we did not take into account the correlation between the bounds of the panel interval-valued data. For future studies, we will introduce variance modeling or joint modeling to obtain better estimators.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The authors express their thanks to the Editor and all referees for their constructive suggestions, which led to improvement of an earlier version of this paper.

Conflict of interest

The authors declare there are no conflicts of interest.