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,827 characters · 15 sections · 35 citation commands
Bi-integrative analysis of two-dimensional heterogeneous panel data model
\baselineskip=18.0pt
Panel (or longitudinal) data models have been widely-used in economics, finance, and many other fields. Panel models exhibit various advantages in combining useful cross-sectional and time series information in the data. Traditional homogeneous panel data model assumes that the slope coefficients are constant across individuals and periods. However, homogeneous assumption maybe too restrictive in many applications. In practice, both cross-sectional and time domain variations are observed. Many panel datasets cover lots of individuals coming from different backgrounds, such as different experimental methods, distinct crowds or geographic locations (census, tract, county, state, etc.), external classification, observable explanatory categories, nested (hierarchical) or non-nested datasets, which lead to heterogeneity across individuals. In addition, the heterogeneity usually represents some individual characteristics that are unobserved, such as the ability of individuals BH2002 in the labor market, the loan willingness for banks cornett2011. On the other side, time-specified coefficients captures unobserved time-varying behavior such as the historical events in the process of democratizationBM2015, technological progress, institutional transformation, or economic transition. For these and other reasons, it is important to take into account for both cross-sectional and temporal heterogeneity in many applications.
Over the last few years, there is a fast growing literature on panel data model with heterogeneous slope coefficients, along two directions. One direction of research assumes that the regression coefficients are time-varying. In this case, the regression coefficients are assumed to be functions of time trend under the nonparametric framework Li2011,Pei2018. In order to deal with the incidental parameter problem, some researchers assume that there exist unknown common breaks Bai2010,Kim2011 or multiple structural breaks QianSu2016,Li2017 of regression coefficients in prior, so that the difference of regression coefficients are successively pairwise sparse in time-dimension. Most of these works use shrinkage or fused penalty method to detect and estimate multiple change points simultaneously. The other direction assumes that the slope coefficients are heterogeneous across individuals, due to some individual-specific characteristics Bester2016,FGPZ2017. This literature includes complete heterogeneity and group-based heterogeneity. Complete heterogeneity assumes that slope coefficients are different across each individuals, and are usually modelled by random coefficient panel data models Wooldridge2005,MW2008, nonparametric model Boneva2015,Vogt2017 or other complete heterogeneous settings, see Pesaran2006, CPT2011, KHJW2017, among others. Group heterogeneous models assume that individuals can be classified into different groups, where the regression coefficients are the same within each group but heterogeneous across groups. Another approach uses finite mixture models in discrete choice panel data models with parametricKS2009 or nonparametric method BC2010. A third method is clustering or integrating the group membership and estimating parameters simultaneously by solving the penalized objective function added the penalty term of slope coefficients between different individuals, such as the C-lasso Su2016,SuJu2018,Huang2018, Panel-CARDSWangSu2018 and others.
An important issue in practice is that individual heterogeneity and time variation may occur simultaneously. For example, the saving-retention coefficients or saving-investment relations are heterogeneous across countries and periods in the famous Feldstein-Horioka puzzle. For this reason, research attention has been recently shifted to panel data model under two-dimensional heterogeneity. The main challenge in this situation is the increasing number of unknown parameters along with N and T. Thus, dimensional decomposition strategy becomes popular for reducing the number of coefficients along the two-dimensions in parametric model. Baltagi2016 extend the common correlated estimation(CCE) with common structural break in the individual-specific slope coefficients. Smith2018 develops a new Bayesian approach to estimate non-common structural breaks in panel regression models. Neal2018, LuSu2019 decomposed the two-dimensional heterogeneous regression coefficients into three parts additively, i.e., the average component, individual-specific and time-varying component respectively. chernozhukov2018 considered interactive pattern of two-dimensional heterogeneous regression coefficients. Another technique identifying the amounts of unknown coefficients with two-dimensions is to assume block structures, for example, OW2020,RobinOkuiWang2020 consider block-based structural slope coefficients with structural breaks and time-invariant grouped structure simultaneously and utilize fused lasso to detect the true pattern. In addition, Su2019 use nonparametric method to allow slope coefficients as a smooth function of time, synchronously considering heterogeneity across units.
In this paper, we propose a method of estimation and inference of panel data model with a more flexible heterogeneous structure. We argue that the concave pairwise fusion approach, proposed by Ma2016,MaHuang2017, can be extended to conduct a Bi-integrative analysis of Group and Cohort Recovery (BIGCORE) for two-dimensional heterogeneous panel structure models. In particular, the BIGCORE method deals with panel structural model where the regression coefficients have block structure, in which the coefficients of observations within the same block are identical, but distinct across blocks in the rectangular arrangement of the two dimensional heterogeneous coefficients. The two-dimensional heterogeneous panel model is general, not only including the existing homogeneous panel data model, panel structural model with pure grouped structure or structural breaks, but also time-varying grouped structure with multiple structural breaks, diverse multiple structural breaks with time-varying grouped structure, and others. Therefore, our model in this paper is more general than the previous research and has great potential in empirical analysis for panel data with grouped and structural changes. Compared to the existing methods, our approach do not require strong assumptions that the coefficients have specific sparse structure, such as common structural breaks, invariant group membership and others and similarly, our approach also need not to determine the number of groups and change points in prior.
We use the ADMM algorithm to solve the optimization problem with double fused penalties. We establish that the estimators are consistent and asymptotically normal. We prove that the estimators have the oracle property in the sense that it is asymptotically equivalent to the infeasible estimator with known two-dimensional heterogeneous structure. Monte Carlo simulations are conducted and show nice sampling properties of our estimators in finite sample. Finally we illustrate the potential of our methods by an empirical application.
The rest of this paper is organized as follows. In Section 2, the two-dimensional panel structure model and the proposed estimation method are presented. In Section 3, an estimation procedure based on the ADMM algorithm is given to solve the optimization problem. In Section 4, we derive the asymptotic properties of the estimator. Section 5 discusses the determination of the penalty parameters and initial values. Section 6 conducts a Monte Carlo simulation. In section 7, we apply the proposed approach to a real dataset. Finally, we conclude the paper.
Notation. we introduce following notations that will be used throughout this paper. $\otimes$ is the Kronecker product, $\circ$ is the Hadamard product, $vec$ is the vectorization operator, $\gg$ denote much greater, the superscript $\top$ denote the transpose of a matrix, $\|\cdot\|$ stands for the Euclidean norm for vector, $\|\cdot\|_{F}$ denote the Frobenius norm of matrix. $\langle a,b\rangle=a^{\top}b$ be the inner product of two vectors a and b with the same dimension. $A^{+}$ denotes a vector obtained from row sums of matrix A. For a given vector $b=(b_{1},\ldots,b_{t})\in\mathbb{R}^{t}$ and a symmetric matrix $A_{t\times t}$, define $\|b\|_{\infty}=\max_{1\leq s\leq t}|b_{s}|$, $\|A\|_{\infty}=\max_{1\leq i\leq t}\sum_{j=1}^{t}|A_{ij}|$, $\|A\|=\|A\|_{2}=\max_{b\in\mathbb{R}^{t},\|b\|=1}\|Ab\|$ and $\|A\|_{2,\infty}=\max_{1\leq i\leq t}\|A_{i,}\|$, where $A_{i,}$ denotes vector of $i$th row of $A$. $\gamma_{\min}(A)$ and $\gamma_{\max}(A)$ be the smallest and largest eigenvalues of $A$ respectively. $\overset{D}{\rightarrow}$ denotes convergence in distribution.
Giving a panel dataset $\{(y_{it},z_{it}):i=1,\cdots,N;t=1,\cdots,T\}$, where $N$ and $T$ correspond to the total number of individuals and periods respectively, we consider the following heterogeneous regression model with two-way varying coefficients:
where $y_{it}\in\mathbb{R}^{1}$ is the dependent variable, $\mu_{it}$ is the time-varying individual fixed effect, $\boldsymbol{z}_{it}=(z_{it(1)},\cdots,z_{it(P-1)})^{\top}$ is $P-1$ dimensional regressors with slope coefficients $\boldsymbol{\eta}_{it}={(\boldsymbol{\eta}_{it1},\cdots,\boldsymbol{\eta}_{it(P-1)})^{\top}}\in\mathbb{R}^{(P-1)}$ that are potentially heterogeneous in both individual and temporal dimensions and $P$ is fixed. $\epsilon_{it}$'s are independent random errors with mean zero and standard error $\sigma$. In this model, the fixed effect $\mu_{it}$ and the slope coefficients $\boldsymbol{\eta}_{it}$ may vary in both individual and temporal dimensions.
In order to identify the unknown two-dimensional regression coefficients, we assume the following block structure, which combines grouped pattern among individuals with cohort structure across time. The time cohort structure is flexible, it allows the existence of common coefficients between nonadjacent time points and also includes structural breaks as special cases.
Let $\boldsymbol{\beta}_{it}=(\mu_{it},\boldsymbol{\eta}_{it}^{\top})^{\top}$, the true unknown block structure can be characterized in the following form:
where $L$ is the unknown number of blocks with unknown partition of rectangle $\{\mathcal{A}_{l}:1\leq l\leq L\}$.
The block structure on regression coefficients given by (ref) is quite general. Apparently, it includes many classical structures, such as: homogeneous constant with $\boldsymbol{\beta}_{it}=\boldsymbol{\alpha}$ for all $i=1,\cdots,N$ and $t=1,\cdots,T$; grouped pattern corresponding to $\boldsymbol{\beta}_{it}=\boldsymbol{\alpha}_{l}$, $(i,t)\in\mathcal{A}_{l}$ for all $t=1,\cdots,T$; and structural breaks with $\boldsymbol{\beta}_{it}=\boldsymbol{\alpha}_{l}$, $(i,t)\in\mathcal{A}_{l}$ for all $i=1,\cdots,N$ and $t_{l1}\leq t\leq t_{lk_{l}}$, where $t_{lj}$ for $j=1,\cdots,k_{l}$ is the time to event in the vertical ordinate of $\mathcal{A}_{l}$. Furthermore, it also includes some irregular heterogeneous structures depicted Figure ((ref)), which illustrates some values of coefficients matrix, where the vertical ordinate represents the individuals and the horizontal ordinate represents temporal points. The block structure in our paper given by (ref) under two-dimensional heterogeneity also includes the time-varying group memberships with common structural break in Figure ((ref)), constant group memberships with non-common structural in Figure ((ref)), identical group memberships with non-common structural breaks that depicted in Figure ((ref)), which is the same as the block structure on regression coefficients considered by OW2020. The block structure in Figure ((ref)) is more complex: group memberships are time-varying in individual dimension and the cohorts\footnote{The pattern that homogeneous coefficient exits in adjacent or nonadjacent temporal dimension is called cohort structure in our paper, which is similar to the grouped structure. Therefore, structural breaks can be regarded as a special case of cohort structure. } are different across groups. Individuals are divided into three groups, where the first quarter and the last quarter belong to identical group with constant and the second quarter and third quarter consist of different cohorts, where parts of nonadjacent time periods have common values. In reality of economics, there may exist other complex block structure under two-dimensional heterogeneity that also can be included in the setting of (ref). To find out the above block structure, we need to estimate three sets of parameters: the block-specified coefficients $\boldsymbol{\alpha}$, the membership $\mathcal{A}_{l}$'s of individuals and time, and the number of blocks $L$.
In the traditional homogeneous panel data models with individual or time fixed effects, a commonly-used technique to tackle incidental parameter problem is 'difference'. However, this strategy of eliminating the heterogeneous fixed effects is invalid whenever the regression slop coefficients are heterogeneous, no matter one-dimensional (i.e., $\boldsymbol{\eta}_{i}$ or $\boldsymbol{\eta}_{t}$) or two-dimensional heterogeneous (i.e., $\boldsymbol{\eta}_{it}$).
In this paper, to estimate the panel regression model with two-dimensional heterogeneous structure given by ((ref)), we propose a bi-integrative procedure via doubly penalized least square with concave fused penalties. Penalized procedures are commonly used for parameter estimation and sparsity structure recovery.
To estimate the parameters $\boldsymbol{\beta}=\left(\boldsymbol{\beta}_{11}^{\top},\boldsymbol{\beta}_{12}^{\top},\cdots,\boldsymbol{\beta}_{1T}^{\top},\cdots,\boldsymbol{\beta}_{N1}^{\top},\boldsymbol{\beta}_{N2}^{\top},\cdots,\boldsymbol{\beta}_{NT}^{\top}\right)^{\top}$, and recover block structure under the fused sparse assumption $\Vert\boldsymbol{\beta}_{it}-\boldsymbol{\beta}_{jt^{\prime}}\Vert=0$ for $(i,t)$ and $(j,t^{\prime})$ belonging to common block $\mathcal{A}_{l}$ for $l=1,\cdots,L$, we consider the following double penalized least squares objective function:
where $\mathcal{P}_{\lambda}(\cdot)$ and $\mathcal{P}_{\gamma}(\cdot)$are pairwise concave penalty functions, for example, SCAD penaltyFan2001 with tuning parameters $\lambda$ \[ \mathcal{P}_{\lambda}(\kappa)=\lambda\int_{0}^{\kappa}\left(1-x/(\lambda\pi)\right)_{+}dx, \] and MCP penaltyZhang2010 with tuning parameters $\gamma$, \[ \mathcal{P}_{\gamma}(\kappa)=\gamma\int_{0}^{\kappa}\text{min}\{1,(\pi-x/)_{+}/(\pi-1)\}dx, \] where the fixed parameter $\pi$ controls the concavity of the penalty function, $\kappa$ represents the pairwise term between individuals or periods.
Notice that the penalty functions in our objective function are composed of two parts: $\mathcal{P}_{\lambda}(\cdot)$ and $\mathcal{P}_{\gamma}(\cdot)$, where $\mathcal{P}_{\lambda}(\cdot)$ classifies individuals into the grouped structure and $\mathcal{P}_{\gamma}(\cdot)$ is used to integrate observations across different (adjacent and nonadjacent) periods into the cohort structure, apparently including detection of structural breaks. $\lambda$, $\gamma\geq0$ are tuning parameters that control the amount of penalty on $\Vert\boldsymbol{\beta}_{it}-\boldsymbol{\beta}_{jt}\Vert$'s and $\Vert\boldsymbol{\beta}_{it}-\boldsymbol{\beta}_{it^{\prime}}\Vert$'s, respectively and determine an estimation path of the coefficient matrix $\boldsymbol{\beta}$, in which it can shrink $\Vert\boldsymbol{\beta}_{it}-\boldsymbol{\beta}_{jt}\Vert$'s and $\Vert\boldsymbol{\beta}_{it}-\boldsymbol{\beta}_{it^{\prime}}\Vert$'s towards zero with large enough values of $\lambda$ or $\gamma$.
For given $\lambda$ and $\gamma$, we define
and the values of $\lambda$ and $\gamma$ can be selected via a properly constructed Bayesian Information Criterion in the following sections. Specifically, for $\gamma\in[\gamma_{\mathrm{min}},\gamma_{\mathrm{max}}]$, $\lambda\in[\lambda_{\mathrm{min}},\lambda_{\mathrm{max}}]$, let the values of $\gamma$ and $\lambda$ be from a grid $\gamma_{\mathrm{min}}=\gamma_{0}<\ldots<\gamma_{M}=\gamma_{\mathrm{max}}$ and $\lambda_{\mathrm{min}}=\lambda_{0}<\ldots<\lambda_{W}=\lambda_{\mathrm{max}}$, respectively. Then for given $\gamma_{m}$, we compute the solution path $\widehat{\boldsymbol{\beta}}(\gamma_{m},\lambda_{w})$ based on the initial value $\widehat{\boldsymbol{\beta}}(\gamma_{m},\lambda_{w-1}).$ Using $\lambda_{w}$ and $\gamma_{m}$, we can compute the $\widehat{L}(\gamma_{m},\lambda_{w})$ distinct values of $\widehat{\boldsymbol{\beta}}_{it}(\gamma_{m},\lambda_{w})$, corresponding to $\{\widehat{\boldsymbol{\alpha}}_{1},\ldots,\widehat{\boldsymbol{\alpha}}_{\widehat{L}(\gamma_{m},\lambda_{w})}\}$. Then we select optimal $\widehat{\gamma}$ and $\widehat{\lambda}$ minimizing a data-driven criterion BIC defined later in ((ref)), i.e., $(\widehat{\gamma},\widehat{\lambda})=\arg\min_{\gamma_{m},\lambda_{w}}\mathrm{BIC}(\gamma_{m},\lambda_{w})$. Given $\widehat{\gamma}$ and $\widehat{\lambda}$, we can calculate the estimates $\widehat{\boldsymbol{\beta}}=\widehat{\boldsymbol{\beta}}(\widehat{\gamma},\widehat{\lambda})$. Thus all the observations can be separated into $\widehat{L}=\widehat{L}(\widehat{\gamma},\widehat{\lambda})$ blocks accordingly, for example, $\widehat{\mathcal{A}}_{l}=\{(i,t):\widehat{\boldsymbol{\beta}}_{it}=\widehat{\boldsymbol{\alpha}}_{l},1\leq l\leq L\}$, and $\{\widehat{\mathcal{A}}_{1},\ldots,\widehat{\mathcal{A}}_{\widehat{L}}\}$ is a mutually exclusive partition of $\{(i,t):i=1,\cdots,N,t=1,\cdots,T\}$.
The construction of solution path with varying double tuning parameters uses the “bottom up”\ strategy - an important and necessary tactic in the literature of fusion penalty method, because the way of block structure recovery shares similarity as that of dendrogram for agglomerative hierarchical clustering.
Since the objective function does not have a closed-form solution, we use the Alternating Direction Method of Multipliers (ADMM) to solve the optimization problem. Let $\boldsymbol{\rho}_{ij,t}=\boldsymbol{\beta}_{it}-\boldsymbol{\beta}_{jt}$ be the difference of two individual-specified coefficients at a given period, and let $\boldsymbol{\delta}_{i,tt^{\prime}}=\boldsymbol{\beta}_{it}-\boldsymbol{\beta}_{it^{\prime}}$ represent the difference of two period-specified coefficients under a given individual, then the objective function is equivalent to
\[ s_{\cdot}t_{\cdot}\quad\boldsymbol{\rho}_{ij,t}=\boldsymbol{\beta}_{it}-\boldsymbol{\beta}_{jt}\quad and\quad\boldsymbol{\delta}_{i,tt^{\prime}}=\boldsymbol{\beta}_{it}-\boldsymbol{\beta}_{it^{\prime}}, \] where $\boldsymbol{\rho}=\{\boldsymbol{\rho}_{ij,t}^{\top},i<j,t=1,\cdots T\}^{\top}$ and $\boldsymbol{\delta}=\{\boldsymbol{\delta}_{i,tt^{\prime}}^{\top},t<t^{\prime},i=1,\cdots,N\}^{\top}$. Under the constraints, the augmented Lagrangian objective function is given by \[
\] where the dual varibles $\boldsymbol{\nu}=\{\boldsymbol{\nu}_{ij,t}^{\top},i<j,t=1\cdots T\}^{\top}$ and $\boldsymbol{\upsilon}=\{\boldsymbol{\upsilon}_{i,tt^{\prime}}^{\top},t<t^{\prime},i=1\cdots N\}^{\top}$ are Lagrangian multipliers, $\psi$ and $\phi$ are fixed tuning parameters.
The ADMM method iteratively updates $\boldsymbol{\beta}$, $\boldsymbol{\rho}$, $\boldsymbol{\delta}$, $\boldsymbol{\nu}$, and $\boldsymbol{\upsilon}$ based on the following three steps: (1) For given values of $\left(\boldsymbol{\beta},\boldsymbol{\nu},\boldsymbol{\upsilon}\right)$, we update $\boldsymbol{\rho}$, and $\boldsymbol{\delta}$. (2) Then, we update $\left(\boldsymbol{\nu},\boldsymbol{\upsilon}\right)$ given other parameters. (3) Finally, the regression parameters $\boldsymbol{\beta}$ can be updated based on $(\boldsymbol{\rho},\boldsymbol{\delta},\boldsymbol{\nu},\boldsymbol{\upsilon})$.
More specifically, given $\boldsymbol{\beta}^{(s)}$, $\boldsymbol{\nu}^{(s)}$, $\boldsymbol{\upsilon}^{(s)}$ at the $s$th step, we obtain $\boldsymbol{\beta}^{(s+1)}$, $\boldsymbol{\nu}^{(s+1)}$, $\boldsymbol{\upsilon}^{(s+1)}$, $\boldsymbol{\rho}^{(s+1)}$, $\boldsymbol{\delta}^{(s+1)}$ in the $(s+1)$th step, by using the following ADMM iterative algorithm. First, we update $\boldsymbol{\rho}^{(s+1)}$ and $\boldsymbol{\delta}^{(s+1)}$, by solving (ref) and (ref) below, i.e.,
where
and
By arguments similar to Ma2016,MaHuang2017, under (ref) and (ref), the elements $\boldsymbol{\rho}_{ij,t}^{(s+1)}$ of $\boldsymbol{\rho}^{(s+1)}$ and the elements $\boldsymbol{\delta}_{i,tt^{\prime}}^{(s+1)}$ of $\boldsymbol{\delta}^{(s+1)}$ are the minimizers of $\frac{\varphi}{2}\Vert\boldsymbol{\xi}_{ij,t}^{(s)}-\boldsymbol{\rho}_{ij,t}\Vert^{2}+\mathcal{P}_{\lambda}(\Vert\boldsymbol{\rho}_{ij,t}||)$ , $\frac{\phi}{2}\Vert\boldsymbol{\vartheta}_{i,tt^{\prime}}^{(s)}-\boldsymbol{\delta}_{i,tt^{\prime}}\Vert^{2}+\mathcal{P}_{\gamma}\left(\Vert\boldsymbol{\delta}_{i,tt^{\prime}}\Vert\right)$, respectively, where $\boldsymbol{\xi}_{ij,t}^{(s)}=\boldsymbol{\beta}_{it}^{(s)}-\boldsymbol{\beta}_{jt}^{(s)}+\varphi^{-1}\boldsymbol{\nu}_{ij,t}^{(s)}$ and $\boldsymbol{\vartheta}_{i,tt^{\prime}}^{(s)}=\boldsymbol{\beta}_{it}^{(s)}-\boldsymbol{\beta}_{it^{\prime}}^{(s)}+\phi^{-1}\boldsymbol{\upsilon}_{i,tt^{\prime}}^{(s)}$. For different threshold operators $\mathcal{P}_{\lambda}(\cdot)$ and $\mathcal{P}_{\gamma}(\cdot)$, the estimates $\boldsymbol{\rho}_{ij,t}^{(s+1)}$ and $\boldsymbol{\delta}_{i,tt^{\prime}}^{(s+1)}$ are updated based on different formula corresponding to that operator. In particular,
where $\varphi$ is turning parameter and \[ S(w,t)=\left\{
\right. \]
Next, we update $\boldsymbol{\nu}^{(s+1)}$ and $\boldsymbol{\upsilon}^{(s+1)}$ by
and
At last, we update the coefficients $\boldsymbol{\beta}^{(s+1)}$ via
where
Minimizing the objective function ((ref)) with respect to $\boldsymbol{\beta}$ is equivalent to minimizing
where $\Omega=(\mathcal{E}\otimes\boldsymbol{I}_{T})\otimes\boldsymbol{I}_{P}$, $\Phi=(\boldsymbol{I}_{N}\otimes\mathcal{D})\otimes\boldsymbol{I}_{P}$, $\mathcal{E}=\{(e_{i}-e_{j}),i<j\}_{\frac{N(N-1)}{2}\times N}^{\top}$ with $e_{i}$ being the $i$th unit vector whose $i$th element is 1 and the remaining elements are 0 and $\mathcal{D}=\{(e_{t}-e_{t^{\prime}}),t<t^{\prime}\}_{\frac{T(T-1)}{2}\times T}^{\top}$ with $e_{t}$ being the $t$th unit vector whose $t$th element is 1 and the remaining elements are 0. The integrative or fusion matrix $\mathcal{E}$ aims to calculate the difference of coefficients between each pairwise individuals, similarly to fusion matrix $\mathcal{D}$ for temporal dimension. Then, we get
where $\boldsymbol{Y}=\left(y_{11},\cdots,y_{1T},\cdots,y_{N1},\cdots,y_{NT}\right)^{\top}$, $\boldsymbol{X}=\mathrm{diag}(\boldsymbol{X}_{1},\cdots,\boldsymbol{X}_{N})$ and $\boldsymbol{X}_{i}=\mathrm{diag}(\boldsymbol{x}_{i1}^{\top},\cdots,\boldsymbol{x}_{iT}^{\top})$ with $\boldsymbol{x}_{it}=(1,\boldsymbol{z}_{it}^{\top})^{\top}$.
Explicit solution of $\boldsymbol{\beta}$ in ((ref)) involves computational burden caused by calculating the inverse of a $NTP\times NTP$ dimensional matrix, especially with large $N$ and $T$. It is also noted that the design matrix $\boldsymbol{X}$, the fusion matrix $\Omega$ and $\Phi$ contain amounts of sparsity part, motivating us to accelerate the calculation process by by saving memory space and employing some equivalent algebra. Let $\widetilde{\boldsymbol{X}}_{NT\times P}=\left(\boldsymbol{x}_{11},\cdots,\boldsymbol{x}_{1T},\cdots,\boldsymbol{x}_{N1},\cdots,\boldsymbol{x}_{NT}\right)^{\top}$ with $\boldsymbol{x}_{it}$ being $P\times1$ regressors under given individual $i$ and period $t$, $\widetilde{\boldsymbol{\beta}}_{NT\times P}=\left(\boldsymbol{\beta}_{11},\cdots,\boldsymbol{\beta}_{1T},\cdots,\boldsymbol{\beta}_{N1},\cdot,\boldsymbol{\beta}_{NT}\right)^{\top}$ with $\boldsymbol{\beta}_{it}$ being $P\times1$ coefficients under given individual $i$ and period $t$, as the dense regressors and coefficients, rearranging the nonzero element of $\boldsymbol{X}$ and $\boldsymbol{\beta}$. Correspondingly, we set another form of dual variables $\widetilde{\boldsymbol{\nu}}$ and $\widetilde{\boldsymbol{\rho}}$, which are $\frac{N\times(N-1)}{2}\times TP$ dimensional matrices, and $\widetilde{\boldsymbol{\upsilon}}$, $\widetilde{\boldsymbol{\delta}}$ are $\frac{T\times(T-1)}{2}\times NP$ dimensional matrix. Therefore, ((ref)) can be also rewritten as
Let $A=\psi\Omega^{\top}\Omega+\phi\Phi^{\top}\Phi=\left[\psi(\mathcal{E}^{\top}\mathcal{E}\otimes\boldsymbol{I}_{T})+\phi(\boldsymbol{I}_{N}\otimes\mathcal{D}^{\top}\mathcal{D})\right]\otimes\boldsymbol{I}_{P}$, $\mathcal{E}^{\top}\mathcal{E}=N\boldsymbol{I}_{N}-1_{N}1_{N}^{\top}$, $\mathcal{D}^{\top}\mathcal{D}=T\boldsymbol{I}_{T}-1_{T}1_{T}^{\top}$. Applying the Sherman-Morrison-Woodbury formula, we can solve the above matrix inverse by \[ (\boldsymbol{X}^{\top}\boldsymbol{X}+A)^{-1}=A^{-1}-A^{-1}\boldsymbol{X}^{\top}(\boldsymbol{I}_{NT}+\boldsymbol{X}A^{-1}\boldsymbol{X}^{\top})^{-1}\boldsymbol{X}A^{-1}. \] We have $A^{-1}=D\otimes\boldsymbol{I}_{P}$, where
Through setting \[ M=\left(I_{NT}+\left(\widetilde{\boldsymbol{X}}\widetilde{\boldsymbol{X}}^{\top}\right)\circ D\right)^{-1}, \] \[ b^{(s+1)}=\widetilde{\boldsymbol{X}}\circ\boldsymbol{Y}+\mathcal{E}^{\top}(\psi\breve{\boldsymbol{\rho}}^{(s+1)}-\breve{\boldsymbol{\nu}}^{(s+1)})+\mathcal{D}^{\top}(\psi\breve{\boldsymbol{\delta}}^{(s+1)}-\breve{\boldsymbol{\upsilon}}^{(s+1)}), \] where $\breve{\boldsymbol{\rho}}^{(s+1)}$ and $\breve{\boldsymbol{\nu}}^{(s+1)}$ are matrices by rearranging $\widetilde{\boldsymbol{\rho}}^{(s+1)}$ and $\widetilde{\boldsymbol{\nu}}^{(s+1)}$ into a $\frac{N(N-1)}{2}\times TP$ matrix whose rows store the fused values between each individuals, sequentially. Similarly, $\breve{\boldsymbol{\delta}}^{(s+1)}$ and $\breve{\boldsymbol{\upsilon}}^{(s+1)}$ are matrices by rearranging $\widetilde{\boldsymbol{\delta}}^{(s+1)}$ and $\widetilde{\boldsymbol{\upsilon}}^{(s+1)}$ into a $\frac{T(T-1)}{2}\times NP$ matrix whose rows store the fused values between each periods, sequentially. Finally, let \[ B^{(s+1)}=\widetilde{\boldsymbol{X}}\circ\left\{ M\left[\widetilde{\boldsymbol{X}}\circ\left(Db^{(s+1)}\right)\right]^{+}\right\} , \] we get \[ \boldsymbol{\beta}^{(s+1)}=\mathrm{vec}\left\{ \left[D\left(b^{(s+1)}-B^{(s+1)}\right)\right]^{+}\right\} . \]
In order to characterize the block structure on regression coefficients, we may first partition the grouped structure among individuals, and then determine the time structure of breaks in each group, as described in Figure ((ref))(a), we call this the group-cohort pattern. Alternatively, we may capture the block structure using a cohort-group pattern (described in Figure ((ref))(b)) which firstly partitions the structural breaks along the time dimension and then determines the group membership in each cohorts. Fortunately, it does not matter which pattern we select, since they depict exactly the same block structure by different structural matrices.
We firstly introduce the group-cohort pattern in detail below. Let $K$ denote the split number of groups and $\mathcal{G}_{0k}$ denote the individual memberships for the $kth$ group for $k=1,\cdots,K$. Further, $R(k)$ denotes the number of blocks in the $kth$ group and $\mathcal{H}_{0r}(k)$ denote the temporal memberships for the $rth$ block in the $kth$ group $r=1,\cdots,R(k)$. Let $\boldsymbol{\Pi}$ denotes the grouped structure across individuals and $\widetilde{\Pi}=\{\pi_{ik},i=1,\cdots,N\}$ denotes an $N\times K$ matrix with $\pi_{ik}=1$ for $i\in\mathcal{G}_{0k}$ and $\pi_{ik}=0$ for $i\notin\mathcal{G}_{0k}$, indicating the group structure. Then, $\boldsymbol{\Pi}_{NTP\times KTP}=(\widetilde{\Pi}\otimes I_{T})\otimes I_{P}$. Furthermore, we let $\boldsymbol{W}$ denote the structural breaks in each group and for $k=1,\cdots,K,$ set $\widetilde{W}(k)=\{w_{tr}\}$ denote an $T\times R(k)$ matrix with $w_{tr}=1$ for $t\in\mathcal{H}_{0r}(k)$ and $w_{tr}=0$ for $r\notin\mathcal{H}_{0r}(k)$, which depicts the structural breaks under each group. Then, let $R=\sum_{k=1}^{K}R(k)$ and $\widetilde{\boldsymbol{W}}_{KT\times R}=\mathrm{diag}(\widetilde{W}(1),\cdots,\widetilde{W}(K))$ and $\boldsymbol{W}_{KTP\times RP}=\widetilde{\boldsymbol{W}}\otimes I_{P}$. As we can see, the number of blocks partitioned by group-cohort is not smaller than the number of real blocks, at least. Therefore, there is also a structural matrix to depict the relationship between group-cohort and real blocks. We use $\boldsymbol{Q}$ to depict the partitioned structure and $\boldsymbol{\eta}$ is the vector of values for split blocks with $L^{0}$ different values under the group-cohort pattern. It is obvious that some of the split blocks belong to same true block. Therefore, we set $\widetilde{\boldsymbol{Q}}_{R\times L^{0}}=\{q_{rl}\}$, with $r=\{1,\cdots,R\}$ and $l=\{1,\cdots,L^{0}\}$, denote an $R\times L^{0}$ matrix with $q_{rl}=1$ for $\eta_{r}=\alpha_{l}$, which depicts a structural matrix that integrate between spitted blocks and then $\boldsymbol{Q}_{RP\times L^{0}P}=\widetilde{\boldsymbol{Q}}\otimes I_{P}$. .As a result, $\boldsymbol{\beta}^{0}=\boldsymbol{\Pi}\boldsymbol{W}\boldsymbol{Q}\boldsymbol{\alpha}^{0}$ and the design matrix with known structural information $\mathbb{X}=\boldsymbol{X}\boldsymbol{\Pi}\boldsymbol{W}\boldsymbol{Q}$.
Similarly, we can also firstly set the cohorts structural matrix denoted by $\bar{\boldsymbol{W}}$ and then use $\bar{\boldsymbol{\Pi}}$ to describe the grouped structure in each cohorts in cohort-group pattern. Furthermore, $\bar{\boldsymbol{Q}}$ denotes the block integration in cohort-group pattern. Obviously, $\boldsymbol{\beta}^{0}=\Pi\boldsymbol{W}\boldsymbol{Q}\alpha^{0}=\bar{\boldsymbol{W}}\bar{\boldsymbol{\Pi}}\bar{\boldsymbol{Q}}\boldsymbol{\alpha^{0}}$, implying that we only need to select one of patterns. The two dimensional heterogeneous structure is more general than that in OW2020, due to the existence of structural matrices $\boldsymbol{Q}$ or $\bar{\boldsymbol{Q}}$, which depict the integration between splitted blocks and can be viewed as a mediator between different pattern.
To study the theoretical results of the proposed block regression estimator, we first investigate the asymptotic properties of the estimator with known block structure. Although, in practice, $L$ is generally unknown and such an estimator is infeasible. This infeasible procedure provides important information to which we should compare our feasible estimator. Let $\boldsymbol{\beta}^{0}$, $\boldsymbol{\alpha}^{0}$, $\mathcal{A}_{0}$ and $L_{0}$ denote the true values of $\boldsymbol{\beta}$, $\boldsymbol{\alpha}$, $\mathcal{A}$ and $L$, respectively. We also let $|\mathcal{A}_{l}|$ to signify the amount of elements in $\mathcal{A}_{l}$. $\mathcal{A}_{\min}=\min_{1\leq l\leq L}|\mathcal{A}_{l}|$ and $\mathcal{A}_{\max}=\max_{1\leq l\leq L}|\mathcal{A}_{l}|$, respectively represent the true minimum and maximum sample sizes among all blocks.
If the underlying block structure $\mathcal{A}=\{\mathcal{A}_{l}:l=1,\cdots,L^{0}\}$ is known, which is equivalent to know the prior information of matrices $\boldsymbol{\Pi}$ and $\boldsymbol{W}$, $\boldsymbol{Q}$, and notice that $\boldsymbol{\beta}=\boldsymbol{\Pi}\boldsymbol{W}\boldsymbol{Q}\boldsymbol{\alpha}$, it is equivalent to consider $\widetilde{\boldsymbol{\alpha}}$ or $\widetilde{\boldsymbol{\beta}}=\boldsymbol{\Pi}\boldsymbol{W}\boldsymbol{Q}\widetilde{\boldsymbol{\alpha}}$. In this case, the post bi-integrative estimator is defined by:
where $\widetilde{\boldsymbol{\alpha}}=(\widetilde{\alpha}_{1}^{\top},\cdots,\widetilde{\alpha}_{L^{0}}^{\top})^{\top}$.
Due to the block structure information, i.e., $\mathcal{A}$ is generally unknown in advance, the block-oracle estimators are infeasible in practice. However, it can shed light on the theoretical properties of the proposed estimators.
For investigating the statistical properties of the induced minimizer $\widetilde{\boldsymbol{\alpha}}$, we impose the following conditions,
Theorem (ref) states the post bi-integrative estimator $\widetilde{\boldsymbol{\alpha}}$ and the estimator $\widetilde{\boldsymbol{\beta}}$ with known block structure are consistent as both $N$ and $T$ $\rightarrow\infty$. Furthermore, the estimator $\widetilde{\boldsymbol{\beta}}$ is uniformly consistent across the samples.
In practice, the block structure is unknown. In this section, we study the asymptotic properties of our proposed estimator with unknown block structure. We show that, under appropriate conditions, the induced local minimizer of the objective function ((ref)) is asymptotically equivalent to the post bi-integrative estimator under a prior knowledge of block structure $\widetilde{\boldsymbol{\alpha}}$.
Let \[ b_{n}=\min_{\substack{(i,t)\in\mathcal{A}_{l}\\ (j,t^{\prime})\in\mathcal{A}_{l^{\prime}} } }\left\Vert \boldsymbol{\beta}_{it}^{0}-\boldsymbol{\beta}_{jt^{\prime}}^{0}\right\Vert =\min_{l\neq l^{\prime}}\left\Vert \alpha_{l}^{0}-\alpha_{l^{\prime}}^{0}\right\Vert \] be the minimum difference of the coefficients between any two blocks. In addition, we give assumption (C3):
The result in Theorem (ref) implies that if the minimal difference of the coefficients between any two blocks is restricted by a lower bound, our proposed double penalized least square estimator can attain the block-oracle estimator and actually recover the true block structure with probability tending to one. Since the local minimizer $\widehat{\boldsymbol{\alpha}}$ of the objective function just attains the block-oracle estimator $\widetilde{\boldsymbol{\alpha}}$ , we can conclude the following corollary.
The asymptotic distribution of the estimator provides a theoretical foundation for further statistical inference, such as the testing of heterogeneity. Next, we present an asymptotic $\chi^{2}$ test for hypothesis based on the estimators $\widehat{\boldsymbol{\alpha}}$. Specifically, we consider the null $H_{0}:\mathcal{B}\boldsymbol{\alpha}=0$ versus the alternative hypothesis $H_{1}:\mathcal{B}\boldsymbol{\alpha}\neq0$, where $\mathcal{B}$ is a $q\times LP$ matrix and $q=$ rank$(\mathcal{B})$. Many important special cases belong to this hypothesis. For example, $H_{0lj}$: $\alpha_{lj}=0$, $l\in\{1,\ldots,L\}$ and $j\in\{1,\ldots,P\}$, which can be used to test the significance of the $j$th component of coefficients in the $l$th block; The null hypothesis $H_{0}:\alpha_{l}-\alpha_{l^{\prime}}=0$, $l,l^{\prime}\in\{1,\ldots,L\}$ can be used to test the existence of coefficients heterogeneity among blocks.
A standard $\chi^{2}$-test statistic for testing $H_{0}$: $\mathcal{B}\boldsymbol{\alpha}=0$ can be constructed as follows:
where $\widehat{\mathcal{V}}=\widehat{\sigma}^{2}(\mathbb{X}^{\top}\mathbb{X})^{-1}$.
Theorem (ref) provides the asymptotic distribution of the test statistic $\mathcal{T}(\mathcal{B})$ under the null hypothesis $H_{0}$. Therefore, the $100(1-\tau)\%$ confidence interval for $\mathcal{B}\boldsymbol{\alpha}$ is given by \[ \mathbb{R}_{\tau}=\left\{ \iota:(\mathcal{B}\widehat{\boldsymbol{\alpha}}-\iota)^{\top}(\mathcal{B}\widehat{\mathcal{V}}\mathcal{B}^{\top}-\iota)^{-1}(\mathcal{B}\widehat{\boldsymbol{\alpha}})\leq\chi_{q}^{2}(1-\tau)\right\} , \] where $\chi_{q}^{2}(1-\tau)$ is the $(1-\tau)$-quantile of the $\chi^{2}$ distribution with $q$ degrees of freedom.
The proposed ADMM algorithm requires an initialization. Initial value matters in accelerating the convergence of the iteration. In this paper, we propose the ridge fusion criterion to select initial parameters, since it has closed-form solution. Let
which can be written in matrix form
where $\lambda^{\ast}$, $\gamma^{\ast}$ are the tuning parameters and chosen as $\lambda^{\ast}=\gamma^{\ast}=0.001$ in determination of initial values. By minimizing objective function ((ref)), the initial value of $\boldsymbol{\beta}$ is given by \[ \boldsymbol{\beta}^{(1)}=\mathrm{vec}\left\{ \{D^{\ast}(\widetilde{\boldsymbol{X}}\circ\boldsymbol{Y}-\widetilde{\boldsymbol{X}}\circ(M^{\ast}[\widetilde{\boldsymbol{X}}\circ(D^{\ast}(\widetilde{\boldsymbol{X}}\circ\boldsymbol{Y}))]^{+}))\}^{\top}\right\} , \] where \[ D^{\ast}=\left\{ (\lambda^{\ast}N+\gamma^{\ast}T)\boldsymbol{I}_{NT}-\left[\psi(1_{N}1_{N}^{\top})\otimes\boldsymbol{I}_{T}+\phi\boldsymbol{I}_{N}\otimes(1_{T}1_{T}^{\top})\right]\right\} ^{-1} \] \[ M^{\ast}=\left(I_{NT}+\left(\widetilde{\boldsymbol{X}}\widetilde{\boldsymbol{X}}^{\top}\right)\circ D^{\ast}\right)^{-1}. \]
The proposed estimation is based on a penalized procedure that entails choices of tuning parameters $\lambda$ and $\gamma$. \ Unsuitable choices of tuning parameters can produce poor estimates. Motivated by Wang2009,Ma2016, we select the optimal tuning parameters $\widehat{\lambda}$ and $\widehat{\gamma}$ by minimizing the following modified BIC:
where $\mathcal{C}_{NT}$ is a constant or depending on $N$ and $T$. Following Ma2016,MaHuang2017, we select $\mathcal{C}_{NT}=\log(NTP)$ in the Monte Carlo simulation and empirical analysis.
For convenience of analysis, we introduce some additional notations. Let $\mathcal{L}=\{1,2,\cdots,L_{\max}\}$, its three subsets $\mathcal{L}_{0}=\{L\in\mathcal{L}:L=L_{0}\}$, $\mathcal{L}_{\_{}}=\{L\in\mathcal{L}:L<L_{0}\}$, $\mathcal{L}_{+}=\{L\in\mathcal{L}:L>L_{0}\}$, represent cases of the true, under and over-fitting bi-integration, respectively. We establish asymptotic validity of the proposed BIC criterion in the following theorem.
In this section, we perform Monte Carlo simulation with $\mathcal{R}$ replications to investigate the finite-sample performance of the proposed bi-integration procedure with two data generating processes with various heterogeneous block structures on regression coefficients, measured by two aspects that one is the evaluation of the estimated regression coefficients and the other is the accuracy of bi-integration or recovery of block structures. We evaluate the performance of the estimated regression coefficients by root mean square error (RMSE) and its bias, measured by $\frac{1}{\mathcal{R}}\sum_{r=1}^{\mathcal{R}}\sqrt{\frac{1}{NTP}\|\widehat{\boldsymbol{\beta}}^{r}-\boldsymbol{\beta}^{0}\|^{2}}$ and $\frac{1}{\mathcal{R}}\sum_{r=1}^{\mathcal{R}}\left[\frac{1}{NTP}\sum_{i=1}^{N}\sum_{t=1}^{T}\sum_{j=1}^{P}(\widehat{\beta}_{itj}^{r}-\beta_{itj}^{0})\right]$ respectively, where $\widehat{\boldsymbol{\beta}}^{r}$ is the estimated coefficients vector in the $r$th replicate.
We evaluate the estimated numbers of blocks $\widehat{L}$ by the percentage (Per) of $\widehat{L}$ equal to the true number of blocks by the proposed BIGCORE procedure, calculated by $\frac{1}{\mathcal{R}}\sum_{r=1}^{\mathcal{R}}I(\widehat{L}^{r}=L^{0})$, where $\widehat{L}^{r}$ is the calculated number of blocks in the $r$th replicate. We also use the extended rand index(ERI), which measures percentage of correctly membership in each blocks. The Rand Index (RI) is used to evaluate the accuracy of clustering, which lies between 0 and 1, where higher values indicate better performance. Motivated by the formation of RI, we can get individual or period-specified RIs, denoted by $\text{RI}_{t}$ or $\text{RI}_{i}$ and define the ERI(T) and ERI(N) as the average of the whole periods and individuals respectively, i.e., ERI(T) $=\frac{1}{T}\sum_{t=1}^{T}\text{RI}_{t}$ and ERI(N)$=\frac{1}{N}\sum_{i=1}^{N}\text{RI}_{i}$. At last we adopt ERI=$\frac{1}{2}${[}ERI(T)+ERI(N){]} to evaluate the accuracy of BIGCORE procedure.
In this sections, we generate the simulated panel data observations $\{y_{it},x_{it}\}$, $i=1,\cdots,N$ and $t=1,\cdots,T$ by two data generating processes (DGP) with different block structures on regression coefficients and set the sample size as $N=20,40,60$ with $T=20,40,60$. In order to present the wide applicability of the proposed bi-integrating procedure, we consider the complex block structure in the example DGP1 and classical grouped structure in DGP2, respectively.
DGP1:(Block structure)
In this example, we generated observations from a two-dimensional heterogeneous panel data model, \[ y_{it}=\mu_{it}+x_{it}\eta_{it}+\epsilon_{it},\qquad i=1,\cdots,N;\quad t=1,\cdots,T. \] Both the time-varying individual fixed effect $\mu_{it}$ and one-dimensional slope coefficient $\eta_{it}$ have the same structure of time-varying group memberships with common structural break depicted in Figure ((ref)). Then two blocks are considered and coefficient vector of the first block is set as $\boldsymbol{\alpha}_{1}=(-2,3)$ and let that of the second one be $\boldsymbol{\alpha}_{2}=(2,5)$ with the components corresponding to fixed effect $\mu_{it}$ and slope coefficient $\eta_{it}$, respectively. The block-based two-dimensional heterogeneous structure can be depicted by group-cohort pattern, such as, under the case of $N=40$, $T=40$, we consider group structure with $\mathcal{G}_{01}=\{1,\cdots,10,31,\cdots,40\}$, $\mathcal{G}_{02}=\{11,\cdots,20\}$, $\mathcal{G}_{03}=\{21,\cdots,30\}$ and the corresponding cohort structures under the assumed group formation are set by $\mathcal{H}_{01}(1)=\{1,\cdots,40\}$, $\mathcal{H}_{01}(2)=\{1,\cdots,19,30,\cdots,40\}$, $\mathcal{H}_{02}(2)=\{20,\cdots,29\}$, $\mathcal{H}_{01}(3)=\{1,\cdots,9,35,\cdots,40\}$, $\mathcal{H}_{02}(3)=\{10,\cdots,34\}$ and $\mathcal{H}_{01}(4)=\{1,\cdots,40\}$. Therefore, the ratio of number of block-specified observations is about $|\mathcal{A}_{1}|:|\mathcal{A}_{2}|\approx3:1$. Other cases setting the sample sizes in the block structure under different combination of $N$ and $T$ has the similar way. The regressor $x_{it}$ is generated by \[ x_{it}=1+0.5\mu_{it}+\epsilon_{it}, \] where $\epsilon_{it}$ was taken from the standard normal distribution.
We consider the settings of homoscedasticity and heteroscedasticity on the error term by respectively generating $\epsilon_{it}\sim\ N(0,\sigma^{2})$ with $\sigma^{2}=0.5$, $\sigma^{2}=1$ and \[ \epsilon_{it}=\sigma_{it}e_{it},\sigma_{it}=\tau(0.05+0.05x_{it}^{2})^{1/2}, \] where $\tau=2$ or $\tau=1$ and $e_{it}\sim N(0,1)$.
DGP2(Grouped structure)
In this example, we consider the performance of proposed BIGCORE analysis in panel data model with grouped individual fixed effect and grouped slope coefficients. The Datasets are generated as: \[ y_{it}=\mu_{i}+x_{it}\eta_{i}+\epsilon_{it}. \] Here, both the fixed effects and slope coefficients have identical grouped structure by randomly dividing the individuals into three groups with the proportion that $|\mathcal{G}_{1}|:|\mathcal{G}_{2}|:|\mathcal{G}_{3}|=3:3:4$, in which the true coefficients are $\boldsymbol{\alpha}_{1}=\{-2,3\}$, $\boldsymbol{\alpha}_{1}=\{2,6\}$ and $\boldsymbol{\alpha}_{3}=\{6,-1\}$, respectively. The regressor $x_{it}$ are generated as \[ x_{it}=1+0.5\mu_{i}+\epsilon_{it}. \] Lastly, the error term is set as in DGP 1.
In the simulation, we select the number of replicates $\mathcal{R}=100$. The grid of values of both tuning parameters $\lambda$ and $\gamma$ is set in the range of {[}0.1, 1.5{]} with step size 0.1. The fact that increasing grid of tuning parameters apparently improves integrative results is unconsidered here due to reducing computational cost. In order to accelerate the convergence of the proposed ADMM algorithm, we regard the converged value under the given combined tuning parameters as the initial value of the next iteration, instead of adopting identical initial values under different combined tuning parameters in the grid.
After one replicate, Figures ((ref)) vividly presents the performance of our proposed BIGCORE analysis in estimating coefficients and block structure in DGP1 and DGP2 settings under $N=40$, $T=40$ and heteroscedasticity with $\tau=2$ and the figure elucidates that the developed method can achieve expected outcome because it can recover the true block or grouped structure correctly with consistent value of coefficient estimators. Another inevitable fact we should admit is that although Figures ((ref)) actually shows ideal simulated results in one replicate, several worse BIGCORE results still occur occasionally, which is reflected by the misintegration that a little of observations may be wrongly partitioned, such as, it may occur in the replicates with the cases of large standard deviation of error term. Specifically, a little of observations originally belonging to red block are mistakenly bi-integrated into the blue block.
For checking the representation of coefficients estimators, we report the RMSE and Bias of the slope coefficient for examples DGP1 and DGP2 in Tables ((ref)) and Table ((ref)), respectively, which elucidates that (i) both SCAD and MCP penalties present similar BIGCORE behaviors in terms of close values of RMSE and Bias, and both are also close to the oracle results, which is the reason we reject arguing the effective combination of the forms of double concave penalties; (ii) with increasing number of individuals or periods, the values of RMSE and the Bias decrease remarkably in all cases; (iii) the values of RMSE and Bias in DGP1 are relatively larger than that in DGP2. It may be owe to the more complex formation of block structure than that of grouped pattern and the group-specified cohort structure should be integrated in DGP1; (iv) the post estimators is recommended due to that it attains much smaller RMSE and Bias, especially as $N$ or $T$ increases. It is also noted that the estimation performance on slop coefficients under post-MCP is usually the same to that of post-SCAD, which is attributed to the common block structure recovery by both penalties.
For evaluating the accuracy of BIGCORE procedure in estimating the number of blocks, Table ((ref)) and Table ((ref)) reports the percentage of the estimated numbers of blocks equal to the true number of blocks by the SCAD and MCP shrinkage procedures under different cases of DGP1 and DGP2. In all cases, the percentage of correctly selecting the number of blocks increases as $N$ and $T$ are enlarged. The two concave penalties SCAD and MCP procedures have similar performance.
Another results under the evaluation criterion extended Rand index, which is used to measure the bi-integration ability of recovering the true underlying structures, are reported in Table ((ref)), Table ((ref)) and the results show that the extended Rand index are mostly close to one, which indicates the effectiveness of the proposed BIGCORE method. The results also imply that the BIGCORE performance was worsen with serious heteroscedasticity such as larger $\sigma^{2}$.
In this section, we apply the proposed BIGCORE procedure to measure the heterogeneous impact of inputs on the economic output. Based on the classical Solow model about economic growth equation, the economic output is mainly determined by technological progress, Capital and Labor, and we establish the following regression given by
where $\mathrm{GDP}_{it}$ denotes the real gross domestic product, technological progress are usually represented by human capital(Hc) in the empirical literature, Ck is physical capital stock and Ngd denotes population growth plus break even investments of $5\%$. The coefficients represent different economic meanings, such as, the slope coefficient $\beta_{2,it}$ is the elasticity of investment on output. In the whole world, countries with different resource endowments and technical power are at different stages of development. For the developing countries, the elasticity of investment on output is larger than that of developed countries according to the law of economic development. From the perspective of period, the improvement of technological progress on economic development is diminishing marginally or remain stable, and the marginal effect shifts to higher level as the emergence of new technologies. Therefore, the slope coefficients are heterogeneous across countries and may exist structural breaks in the long span.
The original data are available form the Penn World Tables 8 and can be directly obtained from the package xtdcce2 of Stata. The dataset contains panel data with 92 countries and its yearly observations from 1963 until 2007. Same as QianSu2016, the observations are averaged by each 5 years. Thus, we finally get panel data with $N=92$, $T=9$ and $P=4$. In addition, we set a grid of turning parameters $\lambda$ and $\gamma$ from 0.2 to 3 with interval of 0.2. Finally, we get two blocks $\hat{L}=2$ shown in Figure ((ref)) that the blue block is denoted by $\mathcal{A}_{1}$ and the red block is denoted by $\mathcal{A}_{2}$. In the estimation process, optimal $\lambda$ and $\gamma$ are selected by 2.6 and 1.6 respectively, according to the modified BIC. The penalized estimators and post estimators based on the estimated block structure are reported in the Table ((ref)) and the results imply that all the penalized estimators are statistically significant and most post estimators are statistically significant, except for $\beta_{2,it}$ and $\beta_{3,it}$ in block $\mathcal{A}_{2}$. What's more, the standard errors of post estimators are smaller than that of penalized estimators, which is consistent to that of the Monte Carlo simulation. In conclusion, the block heterogeneity is remarkable in Solow model based on the BIGCORE method that the elasticity or marginal coefficient are significant across blocks, owing to different endowments across countries, continuous development and progress.
In this work, we consider a general panel data model with two-dimensional heterogeneous coefficients and propose a novel BIGCORE procedure to discover the assumed block structure and obtain the estimators simultaneously. An ADMM algorithm is developed to iteratively solve the objective function with double concave fused penalties. Simulated studies suggest that our method is effective and show fine performance in conducting BIGCORE analysis by correctly estimating the structure and coefficients. A modified Bayesian information criteria is proposed to get rid of the knotty issue that the double tuning parameters are sensitive to the estimators. However, computational complexity is growing with the increasing number of samples and time periods, causing the burden of the ADMM algorithm. What's more, this work assumes that all coefficients have identical block structure under two-dimensional heterogeneity, which motivates us to plan to extend the BIGCORE analysis to the case of covariate-specified block structure. All or other issues, including the consideration of two-dimensional heterogeneous panel model with interactive effects, are worthy to be studied in the further research.
\phantomsection\addcontentsline{toc}{section}{\refname}