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.
48,179 characters · 12 sections · 26 citation commands
Estimation of High-Dimensional Seemingly Unrelated Regression Models
A SUR comprises multiple individual regression equations that are correlated with each other. In our setup, we assume that there are $N$ regression equations that are observed over periods $t=1,2,...,T$. These regression equations are related in the sense that the regression errors of the equation system are correlated.
The SUR estimator originally proposed by zellner1962efficient is a feasible generalized least squares (FGLS) estimator that is based on an estimator of the inverse of the covariance matrix, the precision matrix ($\Omega := \Sigma^{-1}$), of the SUR equation system. Often this estimator is computed in two steps. In the first step, one estimates each equation by the ordinary least squares (OLS) and computes the residuals. In the second step, one computes the FGLS based on the inverse of the sample covariance matrix of the residuals.
It is well-known that when the number of equations, $N$, is large, the FGLS estimator performs poorly (e.g., greene2003econometric). The main reason is when $N$ is large relative to the number of observations $T$, the sample precision matrix (inverting the sample covariance matrix) performs poorly. For example, when $N/T \to c>0$, the empirical covariance matrix is not consistent and would be rank deficient when $N>T$ (e.g., see johnstone2001distribution and hastie2015statistical). In this paper, we revisit the problem of estimating the classical SUR model when the number of regression equations, $N$, is large and even comparable to $T$.
The problem of estimating high dimensional precision matrices has been widely studied in the machine learning and statistical learning literature. A popular method is to estimate the precision matrix with various regularizations -- imposing various restrictions that the precision matrix is sparse. (e.g., see cai2011constrained, negahban2011estimation, friedman2008sparse, lam2009sparsistency.)
It is well known that the sparsity of the precision matrix has a nice interpretation if the underlying distribution is Gaussian - the set of non-zero entries in the precision matrix correspond to the set of edges in an associated Gaussian Markov random field (GMRF) (see hastie2015statistical). Thus imposing sparsity on the precision matrix corresponds to the assumption that not all regression equations are related to each other. In another words, the graph representing which regression equation is related to which, is sparse.\footnote{liu2009nonparanormal relaxed the results to a more general class of distribution named as non-paranormal distribution.} This restriction is reasonable -- a prevailing result in the literature of social and economic networks is that these graphs are sparse, and the degrees of the nodes grow much slower than the network size $N$ (barabasi2016network).
In this paper, assuming sparsity of the true precision matrix $\Omega$, we propose a new estimator called the FGLasso that works when both $N,T \to \infty$ and under some conditions, $N \geq T$. Unlike FGLS which estimates $\Sigma$ by the OLS residuals and then taking the inverse, we directly estimate the precision matrix $\Omega:=\Sigma^{-1}$ using valid high-dimensional techniques. In particular, we directly estimate the precision matrix using the {\em Graphical Lasso} estimator, which is a popular estimator for the high-dimensional precision matrix (see fan2016overview).
Our results are as follows. We find a set of regularity conditions under which our FGLasso estimator is asymptotically equivalent to the (infeasible) GLS estimator uniformly across equations. Moreover, we show that if the maximum nonzero entries per row in $\Omega$ is bounded or grows much smaller than $N$, FGLasso estimator performs well even when $N > T$. In the Monte-Carlo study, we compare the performance of the FGLasso estimator with the OLS estimator, the GLS estimator, and the FGLS estimator, corroborating our findings that our proposed estimator performs better in a high-dimensional setting. In deriving these results, we build upon ravikumar2011high, which show that under certain regularity conditions, the graphical lasso estimator $\widehat{\Omega}_{gl}$ converges to the true $\Omega$ at the rate $\mathcal{O}_p(\sqrt{\log N/T})$, in terms of the element-wise maximum norm, while preserving the sparsity pattern.
The remainder of this paper is organized as follows. Section 2 discusses the SUR model in details, summarizes the OLS, the GLS, the FGLS, and the FGLasso estimators. In section 3, we discuss the results from ravikumar2011high and present the main theoretical results. Section 4 reports Monte Carlo simulation results\footnote{Computation for the work described in this paper was supported by the University of Southern California’s Center for High-Performance Computing (hpc.usc.edu).} and section 5 concludes. All the technical proofs and additional simulation results are provided in the appendix.
For the convenience, we briefly summarize the notation to be used throughout the paper. We denote $s_{\mathrm{min}}(A)$ and $s_{\mathrm{max}}(A)$ as the maximum and the minimum singular values of real valued matrix $A\in \mathbb{R}^{m \times n}$, respectively. The operator norm and Frobenius norm are defined as ${\left\vert\kern-0.25ex\left\vert A \right\vert\kern-0.25ex\right\vert}_{op} = s_{\mathrm{max}}(A)$ and ${\left\vert\kern-0.25ex\left\vert A \right\vert\kern-0.25ex\right\vert}_F=\sqrt{\sum_{i,j}A_{ij}^2}$, respectively. Let ${\left\vert\kern-0.25ex\left\vert A \right\vert\kern-0.25ex\right\vert}_\infty = \underset{i,j}{\mathrm{max}}|A_{i,j}|$ denote the element-wise maximum norm and ${\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert A \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}_\infty = \underset{i=1,2,...,L}{\mathrm{max}}\sum_{j=1}^{N'}|A_{ij}|$the maximum absolute row sum matrix norm. Let $A'$ denote the transpose of $A$ and $\otimes$ denote the Kronecker product. For a real sequence $\{a_n\}_{n=1}^\infty$ and a positive sequence $\{b_n\}_{n=1}^\infty$, we denote $a_n=\mathcal{O}(b_n)$ if there exists a finite constant $C$ such that $|a_n| \leq C b_n$ as $n \to \infty$, and $a_n = \mathcal{O}_p(b_n)$ if $\mathbb{P}(|a_n| \leq C b_n) \to 1$ as $n \to \infty$. We use notation $\Rightarrow$ and $\xrightarrow{p}$ to denote the convergence in distribution and the convergence in probability, respectively.
Suppose we estimate a system of linear equations:
for $i = 1,\cdots,N$ and $t=1,\cdots,T$. Here $X_{it}=(X_{it,1},X_{it,2},\cdots,X_{it,K_i})'$ is a $K_i$- column vector of the regressors for unit $i$, and $U_{it}$ is the unobserved error term. The heterogeneous regression coefficients $\beta_i \in \mathbb{R}^{K_i \times 1}$ are the parameters of interest.
Stacking the observations over $N$ units, let $Y_t = (Y_{1t},...,Y_{Nt})' \in \mathbb{R}^N $, $U_t = (U_{1t},...,U_{Nt})' \in \mathbb{R}^N $, $X_t = diag(X_{1t},...,X_{Nt}) \in \mathbb{R}^{N \times \sum_{i=1}^N K_i}$, and $\beta = (\beta_1',...,\beta_N')' \in \mathbb{R}^{\sum_{i=1}^N K_i}$. The system of the equations in ((ref)) can be expressed as
Alternatively, stacking the observations in ((ref)) over $t$, we can also express the system of the equations in ((ref)) as
where $Y_i=(Y_{i1},Y_{i2},\cdots,Y_{iT})'\in \mathbb{R}^T$, $X_i=(X_{i1},X_{i2},\cdots,X_{iT})'\in \mathbb{R}^{T \times K}$ and $U_i=(U_{i1},U_{i2},\cdots,U_{iT})'\in \mathbb{R}^T$.
In a matrix form, we can write the model as
where $Y=(Y_1',Y_2',\cdots,Y_N')' \in \mathbb{R}^{NT}$, $U=(U_1', U_2',\cdots, U_N')'\in \mathbb{R}^{NT}$, \[X=
_{TN \times \sum_{i=1}^N K_i} and \quad \beta=
.\]
In this paper, we assume the classical linear system equation assumptions:
The conditions in Assumption (ref) are quite classical in the SUR literature. The first condition excludes the case all the regressors are identical. In this case, it is well known that the OLS estimator becomes efficient and there is no gain of using the information in the equation system. The second condition assumes that the regressors are strictly exogenous and it excludes regressors like lagged dependent variables. This assumption may be restrictive in many applications, but we need this assumption as a technical regularity condition that is required for proving the asymptotic properties of $GLS$ estimator and $FGLasso$ estimator when $N,T\to \infty$ (see details in Appendix A.). The third condition is the homoskedasticity assumption. This condition may be relaxed to allow for conditional heteroskedasticity but at the cost of technical complexity of the asymptotic results of the paper. We assume it just for simplicity in deriving asymptotic results.
Without loss of generality, in the rest of the paper, we assume the number of regressors of each unit is the same, i.e. $K_i=K,$ $i=1,2,\hdots,N$, a constant number. The case in which $K_i's$ are different can be easily extended. The results remain the same as long as $K_i=\mathcal{O}(1)$ for all $i$.
In this section, we first briefly summarize the OLS, the GLS, and the FGLS estimators of $\beta$ in the SUR model. Then we introduce the FGLasso estimator.
The OLS estimator is defined as
It is equivalent to the OLS estimators of individual equations,\\ $\widehat{\beta}_{OLS}=(\widehat{\beta}_{1,OLS}',\widehat{\beta}_{2,OLS}',\cdots,\widehat{\beta}_{N,OLS}')'$, where $\widehat{\beta}_{i,OLS}=\big(X_i' X_i\big)^{-1} \big(X_i' Y_i\big)$ for $i=1,2,...,N$.
In his seminar paper, zellner1962efficient proposed the SUR estimator to improve the OLS estimator by exploiting the correlation in the equation system. Suppose that $\Sigma$ is known. As earlier, define the precision matrix as $\Omega:=\Sigma^{-1}$. Then GLS estimator is defined as
In most applications, however, $\Sigma$ and $\Omega$ are not known. A FGLS estimator (see details in greene2003econometric) is defined by replacing the unknown $\Sigma$ with the consistent estimator. A widely used estimator of $\Sigma$ is $\widehat{\Sigma}=\frac{1}{T}\sum_{t=1}^T\widehat{U}_t\widehat{U}_t'$ and $\widehat{U}_t$ is the $OLS$ residuals, that is, $\widehat{U}_t=Y_t-X_t'\widehat{\beta}_{OLS}$. Then
The FLGS estimator in ((ref)) suffers from two major problems when $N$ is large. Suppose that $T>N$, but both $T$ and $N$ are large and in the same order. It is known that the estimator $\widehat{\Omega}=\widehat{\Sigma}^{-1}$ behaves poorly. Further, $\widehat{\Sigma}^{-1}$ is only well defined when $T \geq N$. When $T$ is less than $N$, $\widehat{\Sigma}$ is rank deficient and therefore not invertible.
Our estimator is motivated by these two issues. Suppose that $\Omega$ is sparse. In this case, we propose FGLasso estimator by replacing $\widehat{\Omega} = \widehat{\Sigma}^{-1}$ in ((ref)) with a graphical lasso estimator, $\widehat{\Omega}_{gl}$, where
${\left\vert\kern-0.25ex\left\vert \Omega \right\vert\kern-0.25ex\right\vert}_{1,\textrm{off}}=\sum_{i\neq j}^N |\Omega_{ij}|$ and $\lambda_n>0$ is a penalization parameter which is often chosen by a cross-validation method (e.g., see friedman2008sparse). More specifically,
In the next section, under a certain restriction on the sparse structure of $\Omega$, we show that $\widehat{\Omega}_{gl}$ is consistent even when $T\leq N$, and $\widehat{\beta}_{FGLasso}$ is consistent and has similar asymptotic properties as $\widehat{\beta}_{GLS}$.
The sample properties of the FGLasso estimator is mainly dependent on the sample properties of $\widehat{\Omega}_{gl}$ in ((ref)). Intuitively, if $\widehat{\Omega}_{gl}$ is close to true $\Omega$ in some metric, then $\widehat{\beta}_{FGLasso}$ would also be close to $\widehat{\beta}_{GLS}$. In section 3.1, we briefly summarize results in ravikumar2011high regarding the properties of $\widehat{\Omega}_{gl}$. In section 3.2, we show the consistency and asymptotic properties of $\widehat{\beta}_{FGLasso}$.
Define $\Gamma^*=\Omega \otimes \Omega$ and let $E(\Omega)$ be the edge set including self-link, i.e $E(\Omega):=\{(i,j)|\Omega_{ij}\neq 0\}$. Denote $S^c$ be the complement of S set. For set $P,P'$, we also use $\Gamma^*_{PP'}$ to denote the $|P|\times |P'|$ matrix with rows and columns of $\Gamma^*$ indexed by $P$ and $P'$ respectively. $\kappa_{\Sigma}:={\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert \Sigma \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}_\infty$ and $\kappa_{\Gamma}:={\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert (\Gamma^*_{SS})^{-1} \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}_\infty$.
The first condition assumes that $D_N$ is an upper bound of the number of nonzero elements in $\Omega$. It usually refers to the maximum degree of nodes in graph theory literature. In this paper, we do not put a restriction that a finite constants bounds $D_N$. We allow $D_N$ to increase slowly to infinity as $N$ increases. The second condition assumes that the distribution of $U_{it}$ has a thin tail like the Gaussian distribution. The third assumption is usually referred as the incoherence condition, which guarantees the exact recovery of $\widehat{\Omega}_{gl}$. The condition that $\kappa_\Sigma$ remains constant in the third assumption indicates that the singular value of true precision matrix $\Omega$ is lower bounded\footnote{Note $s_{min}(\Omega)=\frac{1}{s_{max}(\Sigma)}\geq \frac{1}{{\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert \Sigma \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}_\infty}$.}. The assumption that the rest of the parameters $(\kappa_{\Gamma^*},\alpha,\sigma)$ remain constant is only for simplicity, more detailed results and discussion can be found in Corollary 1 in ravikumar2011high.
The first result of Lemma (ref) guarantees that, the error between each element of $\widehat{\Omega}_{gl}$ and $\Omega$ shrinks uniformly at a rate of $\sqrt{\log N/T}$, meaning that as long as $T$ increases faster than $\log N$, the error will go to zero. The maximum degree $D_N$ plays an important role as it determines the lower bound of sample size $T$. If $D_N$ is bounded or increases much slower than $N$ (for example, $D_N=\log N$), then it is possible that the properties in Lemma (ref) hold for $T<N$.
The second result of Lemma (ref) shows the exact recovery property, meaning that $\widehat{\Omega}_{gl}$ from ((ref)) remains similar sparsity structure of $\Omega$. The non-edge set $E^c(\Omega) \subseteq E^c(\widehat{\Omega}_{gl})$, that is, if $\Omega_{ij}=0$, then $\widehat{\Omega}_{gl,ij}=0$ wp1. Therefore, consider matrix $\Delta_\Omega:=\Omega-\widehat{\Omega}_{gl}$, the maximum nonzero entries per row is at least $D_N$.
In this section, we discuss the consistency and asymptotic properties of $\widehat{\beta}_{FGLasso}$ defined in ((ref)). As the main result, we show that $\widehat{\beta}_{FGLasso}$ and $\widehat{\beta}_{GLS}$ are asymptotically equivalent. For this, we assume the following assumption.
The first assumption in Assumption (ref) regulates the true precision matrix $\Omega$. In this paper, we assume its smallest eigenvalue does not shrink to zero and the largest elements does not explode as $N,T \to \infty$. Conditions (ii) and (iii) ensure all individual OLS estimators, $\widehat{\beta}_{i,OLS}, i=1,...,N$ are well defined.
The uniform upper bound $C_N$ in condition (iv) controls the second moment of $X_t$. In this paper, we do not restrict $C_N$ to be bounded, and allow $C_N \to \infty$ slowly as $N \to \infty$. As we will show later, $C_N$ is one of the key parameters that determine the performance of graphical lasso estimator.
Some examples of the data generating processes of $\{X_{it}\}$ where $C_N=\mathcal{O}_p(1)$ are as follows:
For the rest of the paper, we focus on the general situation thus do not put any constrain on $C_N$, i.e, we allow $C_N\to \infty$ as $N\to \infty$.
The sparsity parameter $D_N$, nonzero entries per row in $\Omega$, determines the minimum requirement of sample size $T$ and the convergence rate. As shown in the proof, the accuracy of $\widehat{\beta}_{Fglasso}$ in terms of $\widehat{\beta}_{GLS}$ depends heavily on the row norm of $\widehat{\Omega}_{gl}-\Omega$, which is bounded by $\mathcal{O}_p \left( D_N \sqrt{\frac{\log N}{T}} \right)$.
If $C_N$ and $D_N$ are bounded, then ${\left\vert\kern-0.25ex\left\vert \widehat{\beta}_{FGLasso}-\widehat{\beta}_{GLS} \right\vert\kern-0.25ex\right\vert}_\infty \leq \mathcal{O}_p \left( \frac{\sqrt{N}\log N}{\sqrt{T}} \right)$. As long as $\sqrt{T}$ grows faster than $\sqrt{N}\log N$, asymptotically, $\widehat{\beta}_{FGLasso}$ is going to perform similarly to $\widehat{\beta}_{GLS}$.
Instead of $D_N=\mathcal{O}_p(1)$, if we assume $D_N$ grows as $N$ but in a much slower rate, say $D_N=\mathcal{O}_p(\log N)$ as suggested by barabasi2016network in network literature, then $\| \widehat{\beta}_{FGLasso}-\widehat{\beta}_{GLS} \|_\infty = \mathcal{O}_p \left(\frac{\sqrt{N}(\log N)^{3}}{\sqrt{T}}\right)$.
Next we will discuss the asymptotic property between $\widehat{\beta}_{Fglasso}$ and true $\beta$. Before we present the main result, Theorem (ref), we first show an asymptotic property of $\widehat{\beta}_{GLS}$.
In particular, let $b=e_i \in \mathbb{R}^{KN}$, the column vector that only $i'th$ element is $1$ and $0$ otherwise. Then for each element of $\widehat{\beta}_{GLS,i}$ $(i=1,2,\cdots,KN)$, we have:
Combining the results from Propositions (ref) and (ref), we deduce that if $T$ grows fast enough compared with $(D_N,N)$, then $\widehat{\beta}_{Fglasso}$ is asymptotically equivalent with $\widehat{\beta}_{GLS}$, therefore the distribution of $\widehat{\beta}_{Fglasso}$ tends to a normal distribution asymptotically. Summerizing this, we provide the following theorem as the main theoretical result of the paper.
Similar to Proposition (ref), let $b=e_i \in \mathbb{R}^{KN}$, Theorem (ref) implies that each element of $\widehat{\beta}_{FGLasso,i}$ $(i=1,2,\cdots,KN)$ satisfies:
In this section, we discuss finite sample properties of the FGLasso estimator using Monte Carlo simulation experiments. For different pairs of $(N,T)$, we generate data and calculate FGLasso as well as OLS, GLS and FGLS estimators. Also, we compare their distances to the true coefficients $\beta$ based on the element-wise maximum norm and root mean square error (RMSE).
The DGP is \[ Y_t=X_t'\beta+U_t \quad t=1,2,....,T, \] where $K=1$, $X_{it} \sim N(0,1)$ for all $t$ and $i$, $\beta_{L\times 1} \sim U[-1,1]$. Further, let $U_{t} \sim N(0, \Omega^{-1})$ where $\Omega \in \mathbb{R}^{N\times N}$ is generated from the following four different designs:
The first two designs generate a sparse precision matrix $\Omega$ with a certain pattern. Specifically, the number of nonzero entries per row is always $3$ in band structure, and is at most $4$ in four-nearest neighbor lattice.
The precision matrix $\Omega$ from the AR(1) design can be seen as a special case of band graph as design 1, the values of the entries exponentially decay as they move away from the diagonal. When $N$ is small, $\Omega$ is relatively dense but becomes a sparse matrix when $N$ is large.
The precision matrix $\Omega$ from the last design is dense. Instead of precision matrix to be sparse, we consider the case that covariance matrix $\Sigma$ is sparse and has a banded structure.
We estimate $\widehat{\Omega}_{gl}$ by solving ((ref)) using algorithm proposed by friedman2008sparse. We choose $\lambda_n$ by the 5-fold cross validation.\footnote{More precisely, in each replication, we divide the $T$ samples into 5 folds and use four of them as the training data set and one as the validation set. With each choice of $\lambda_n$, we estimate the $\widehat{\beta}_{Fglasso}$ estimators using the training data, then plug them into the validation set and calculate the mean squared error. We choose $\lambda_n$ that minimizes the averaged MSE.}
For each experiment, we fix $T=200$ and let $N \in \{ 50,100,200,300,400 \}$\footnote{For four-nearest neighbor lattices design, $N$ can only be square number, so we choose $N=\{49,100,196,289,400\}$. }. For each pair of $\{N,T\}$, we compare $\| \widehat{\beta}_{OLS}-\beta\|$, $\| \widehat{\beta}_{GLS} - \beta\|$, $\|\widehat{\beta}_{FGLS}-\beta\|$ and $\| \widehat{\beta}_{FGLasso}-\beta \|$ in terms of the element-wise maximum $(l_\infty)$ norm and the root mean square error ($RMSE$)\footnote{Here, the $l_\infty$ norm is the element-wise maximum norm $\|\cdot \|_\infty$, and the RMSE is defined as $\| \cdot \|_F / \sqrt{KN}$.}. Table (ref) reports the results$(\times 100)$ on the average of 100 replications, as well as the number of times out of 100 that $\widehat{\beta}_{FGLasso}$ outperforms $\widehat{\beta}_{FGLS}$.
First, the results in Table (ref) confirms that GLS, as the infeasible efficient estimator, performs better than other estimators and has a lower standard deviation. It also confirms the well-known results that when $N\leq T$, FGLS estimator exists and performs closer to GLS when $N$ is small relatively but gets worse as $N$ rises. In the dense design, when the efficiency gain of GLS is little, FGLS behaves worse than OLS.
For the first two designs where $\Omega$ is exactly sparse, the FGLasso estimator outperforms the FGLS estimator even when $N$ is relatively small. Moreover, it maintains the good performance when $N$ increases even when $N>T$. For example, in the band structure where $N=300$ and $T=200$, the $l_\infty$ norm of the $(\widehat{\beta}_{FGLasso}-\beta)$ is $0.2656$ with standard deviation $0.035$, very close to the GLS estimator, which are $0.2242$ and $0.0314$, respectively.
When $N$ large, the FGLasso estimator in the third AR(1) design behaves similarly to the FGLasso in the band design. It is because the entries per row decay exponentially as it moves away from the diagonal, so $\Omega$ has a sparseness structure when $N$ is large. When $N$ is relatively small, both the FGLS and FGLasso estimators perform well and are close to the GLS estimator, but there is no significant evidence of the advantage of the FGLasso over the FGLS in finite samples. For example when $N=50$ and $T=200$, among 100 simulations, the number of times the FGLasso estimator beats the FGLS estimator is 44 times in $l_\infty$ norm, and 51 times in the RMSE.
In the fourth design, when the covariance matrix has a band structure, the efficiency gain of the GLS estimator is limited. Though the FGLasso estimator performs better than the FGLS estimator, their performances are not significantly better compared to the OLS estimator.
In general, when $\Omega$ is sparse, the FGLasso estimator outperforms the FGLS estimator and behaves closely to the infeasible GLS estimator even when $N>T$. When $\Omega$ is not sparse, it requires larger $T$ (or less $N$) for the FGLS estimator to perform well, as the theory predicted.
This paper proposes a new estimator $\widehat{\beta}_{FGLasso}$ in order to deal with high dimensional SUR model. We show that under certain conditions, as $N,T$ goes to infinity at a certain rate, our FGLasso estimator is asymptotically equivalent to the GLS estimator, thus more efficient than the OLS estimator. Further, if the nonzero entries per row $(D_N)$ in the precision matrix grows much slower than the number of the equation $N$, the FGLasso estimator works well even when the number of the equations in the system $(N)$ is greater than sample size ($T$).
The key assumption under which the FGLasso estimator performs well is the sparsity of true precision matrix $\Omega$. As our knowledge, there is no consensus in the literature on how to test the true precision matrix is sparse or not, and we admit that the applicability our estimator to any general economics data remains questionable if the sparsity condition does not hold. However, we believe this is the price to pay in order to recover the $N$ by $N$ precision matrix with a limited sample size.