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.
40,826 characters · 10 sections · 28 citation commands
Estimation of Graphical Models using the $L_1,2$ Norm
{\color{black} The Gaussian graphical model is a graph summarizing the conditional independence relationships among a group of random variables. Suppose that $(X_{1},\dots,X_{p})$ is distributed with $\mathcal{N}(\bm{0},\Sigma_{0})$. If the $(i,j)$-th entry of the precision matrix $\Sigma^{-1}_{0} \equiv \Omega_{0}$ is zero, i.e., $\Omega_{0,ij}=0$, then it is known that $X_{i}$ and $X_j$ are independent conditional on all other $X_{k}$, $k\neq i,j$ (e.g., Chapter 9 of Hastie2015statistical). The Gaussian graphical model is then obtained by identifying the zeros of the precision matrix, and letting nodes $i$ and $j$ be linked if and only if the $(i,j)$-th entry of $\Omega$ is non-zero.
In our empirical illustration, the variables $(X_{1},\dots,X_{p})$ represent the annual investment or the residuals of investment equations of the $p$ number of firms: $(1,\dots,p)$. If two firms $i$ and $j$ are {\em not} linked in this graphical model, it means that the investment decisions of firms $i$ and $j$ are independent of each other, conditional on all other firms $X_{k}$, $k\neq i,j$. Hence the investment decisions of firms $i$ and $j$ do not directly affect each other. On the other hand, if two firms $i$ and $j$ are linked in the graphical model, then the investments of firms $i$ and $j$ have {\em direct} effects on each other, without the mediation of other firms. In recent applications in economics, for instance, giudici2016graphical uses graphical modeling to obtain the network of international financial flows.
It is well known that the inverse of the standard sample covariance estimator, $\hat{\Sigma}^{-1},$ performs poorly particularly in recovering the underlying graphical structure. The reason is that the sample precision matrix produces a dense matrix without any zeros, and hence the resulting graphical model is always a complete network where all nodes are linked to all other nodes.
A widely-used technique to estimate graphical models is called Graphical Lasso (GLASSO), which aims to estimate the precision matrix $\Omega$ while imposing a sparsity constraint on $\Omega$. For an overview of this topic, see Hastie2015statistical,cai2016estimating,fan2016overview. In this paper, we propose an estimation method for graphical models, which works by estimating the precision matrix, $\Omega_0$, taking into account the {\em structure} of sparsity in $\Omega_0$. }
The standard GLASSO approach estimates the precision matrices while controlling for the sparsity of $\Omega_0$, i.e. the number of zeros in $\Omega_0$. However it does not consider when these zeros could be distributed in a certain pattern or structure in $\Omega_0$. For instance, in many economic applications it is reasonable to think that in the graphical model, some nodes have very few links, while some hub-like nodes have many links to other nodes. Such settings are common in economic applications: for instance, we commonly observe the so-called small-world properties in social and economic networks (jackson2008social).
The main feature of our estimation method is to use the $L_{1,2}$-norm penalty, which corresponds to the second moment of the weighted degree of nodes, i.e. $\sum_{k=1}^{p} |\Omega_{ik}|$, for $i=1,\dots, p$. Notice that Lasso uses the $L_{1,1}$-norm penalty, which is proportional to the first moment of the weighted degree distribution. We derive the asymptotic distribution of the new estimator, called Structured Graphical Lasso (SGLASSO) when the sample size $T \rightarrow \infty$ and the dimension of $\Omega_0$, $p$, is fixed.
The main theoretical finding is that {\color{black}when the network size $p$ is fixed}, SGLASSO is asymptotically equivalent to an infeasible GLASSO problem which penalizes an entry of $\Omega$ proportionally according to the influence of the entry, where the influence of an $(i,j)$-th entry is defined as the sum of the true weighted degrees of nodes $i$ and $j$ in the graphical model. Specifically, SGLASSO is asymptotically equivalent to a modified GLASSO where the element-by-element penalty factor for the $(i,j)$-th entry of $\Omega$ is proportional to $(d_{i} + d_{j})$, where $d_{i} = \sum_{k=1}^{p} |\Omega_{0,ik}|$ is the {\em true} weighted degree of node $i$. Note that $d_{i}$ is unknown to us, and hence it is infeasible to be implemented as a GLASSO estimator.
{\color{black}By means of Monte Carlo simulations, we demonstrate} that SGLASSO performs better than GLASSO in finite samples. In terms of estimating the overall precision matrix, SGLASSO achieves lower Kullback-Leibler and Frobenius losses. Moreover, in terms of recovering the structure of the graphical model, we also show that SGLASSO outperforms GLASSO.
Overall, we should use the $L_{1,2}$ norm in cases where we are more interested in the relationships among influential nodes. For instance, in a financial network that exhibits a core-periphery structure, we are more interested in the relationships between the major core banks, and less interested in the relationships between the minor peripheral banks. Knowing how shocks cascades and propagates among influential banks is key to understanding financial contagion (elliott2014financial).
However if we were to remain completely agnostic about the empirical context, there are still reasons for using the $L_{1,2}$ norm. Here, the virtue of the $L_{1,2}$ norm is that it allows us to express lasso penalties that scale automatically according to degree influences (this insight is due to our main theoretical result). In practice, when we use the $L_{1,2}$ norm in conjunction with the $L_{1,1}$ norm, we are allowing for lasso penalties that have a constant factor and an increasing factor over the nodes' degrees. This can only improve fit. The analogy is that the $L_{1,1}$ norm is a linear first-order lasso penalty, and the $L_{1,2}$ norm is a quadratic lasso penalty. Moreover, our simulation shows that using just the $L_{1,2}$ norm can achieve significant improvements over the $L_{1,1}$ norm. This suggests that it is important to allow for lasso penalties that scale according to nodes' degrees, and not restricted to a constant lasso penalty.
The paper is organized as follows. Section 2 introduces the model, the $L_{1,2}$-norm penalty, and a motivating example. Section 3 presents the theoretical results. {\color{black}Section 4 provides a brief summary of the Monte Carlo simulation results. The full simulation results are available in the Supplemental Information. Section 5 contains an empirical illustration.} The Appendix contains technical proofs and derivations.
{\color{black}This paper is related to literature of estimating high-dimensional inverse covariance matrices. In addition to the aforementioned surveys, one can refer to banerjee2008model,friedman2008sparse,yuan2007model,rothman2008sparse,ravikumar2011high.}
A different class of procedure frames the estimation of graphical models as nodewise or pairwise regressions, where the Lasso or Dantzig selector can then be used to achieve variable selection and high-dimensional regularization. The relevant papers belonging to this class are cai2011constrained,cai2016estimating,meinshausen2006high,peng2012partial,ren2015asymptotic. For ultra high-dimensional Gaussian graphical models, a recent approach is the Innovated Scalable Efficient Estimation proposed by fan2016innovated.
lam2009sparsistency study a general class of estimators that nests the GLASSO. In particular, they consider $\operatornamewithlimits{argmax}_{\Omega\succ0} \big\{\log{\det \Omega} - \mathrm{tr}(\Omega \hat{\Sigma}) - \sum_{i \neq j} p_{\lambda}(|\Omega_{ij}|)\big\}$, where $p_{\lambda}(\cdot)$ is a penalty function that depends on a regularization parameter $\lambda$. The SGLASSO estimator does not belong to this class because our penalty term $\|\Omega\|_{1,2}^{2}$ is not additively separable in the entries of $\Omega$. Hence, SGLASSO cannot be analyzed using the framework proposed in lam2009sparsistency.\footnote{It is possible to incorporate other penalty functions into our estimator. Consider: $\hat{\Omega}_{\text{sc-glasso}} := \operatornamewithlimits{argmax}_{\Omega\succ0} \big\{\log{\det \Omega} - \mathrm{tr}(\Omega \hat{\Sigma}) - \|p_\lambda(\Omega)\|^{2}_{1,2} \big\}$, where $p_\lambda(\Omega)$ denotes the matrix (of the same dimension as $\Omega$) such that $p_\lambda(\Omega)_{ij} = p_\lambda(\Omega_{ij})$. For instance, we can then let $p_\lambda(\cdot)$ to be the Smoothly Clipped Absolute Deviation (SCAD) penalty function (fan2001variable).}
Another type of mixed norm that is widely used in the high-dimensional linear regression setting is the $L_{2,1}$-norm, also known as the Group Lasso (friedman2010note,yuan2006model), which is useful when parameters are naturally partitioned into disjoint groups, and we would like to achieve sparsity with respect to whole groups. In contrast, our proposed $L_{1,2}$-norm is distinct from the $L_{2,1}$-norm, and has not been studied in the context of covariance estimation.
Let $\bm{X} = (X_{1},\dots,X_{p}) \in \mathbb{R}^{p}$ be a $p$-dimensional vector distributed according to a multivariate distribution with mean zero and covariance matrix $\Sigma_{0}$. Our goal is to obtain the graphical model corresponding to $\bm{X}$, which is defined as follows.
The (Gaussian) graphical model of $\bm{X}$ is an undirected graph $G = (V,E)$. The set of nodes (or vertices) of this graph is $V = \{1,\dots,p\}$. Each node $i$ corresponds to a random variable $X_{i}$. The set of links (or edges) $E$ is such that $\{i,j\} \notin E$ if and only if $X_{i} \perp X_{j}$ conditional on all $X_{k}, k \in V \setminus \{i,j\}$. That is, the graphical model is defined such that there is {\em no} link between nodes $i$ and $j$ in the graph if and only if $X_{i}$ and $X_j$ are independent conditional on all other nodes besides $i$ and $j$. As such the graphical model $G$ summarizes the pairwise conditional independence relationships among $(X_{1},\dots,X_{p})$. The existence of a link in $G$ amounts to conditional dependence between the two nodes given all other ones.
Estimation of graphical models is based on the following known fact: if $\bm{X}$ is Gaussian, then the non-zeros in its precision matrix (or inverse covariance matrix) $\Omega_{0}:=\Sigma_{0}^{-1}$ corresponds exactly to links in the graphical model (see Hastie2015statistical,whittaker2009). That is, there is a link between nodes $i$ and $j$ if and only if the $(i,j)$-th entry of $\Omega_{0}$ is non-zero. Similarly the lack of a link $\{i,j\}$ between a pair of nodes $i$ and $j$ is equivalent to the $(i,j)$-th entry of $\Omega_{0}$ being zero.
Now suppose we observe $T$ i.i.d samples ${\bm X}_{t}, t=1,\dots,T$ from $\bm{X} = (X_{1},\dots,X_{p})$. We now introduce the GLASSO estimator, as well as our proposed SGLASSO estimator, that allows us to recover the graphical model from the samples. The output of these estimators is a sparse precision matrix $\hat{\Omega}$. From $\hat{\Omega}$, the estimated graphical model is then constructed by including a link between nodes $i$ and $j$ if and only if $\hat{\Omega}_{i,j} \neq 0$.
For a $m \times n$ matrix $A$, define the mixed $L_{p,q}$ norm as follows:
The following $L_{1,2}$ norm in equation ((ref)) will play a crucial role in our estimator:
Our {\bf Structured-GLASSO} (SGLASSO) estimator is defined in equation ((ref)) below, where $\hat{\Sigma}$ is the sample covariance matrix calculated from the sample $({\bf X}_{t})_{t=1}^{T}$.
In comparison, the {\bf GLASSO} estimator (banerjee2008model,friedman2008sparse,yuan2007model,fan2015overview) is defined in equation ((ref)). While SGLASSO uses the $L_{1,2}$ mixed norm as a penalty to the likelihood, GLASSO in equation ((ref)) uses the $L_{1,1}$ norm, which is the familiar lasso penalty (sum of the absolute values of the entries of $\Omega$).
Another definition of the GLASSO is one where only $\Omega^{-}$ is penalized, where $\Omega^{-}$ denotes the matrix $\Omega$ where its' diagonal entries are set to zero. While yuan2007model,rothman2008sparse use equation ((ref)), other authors such as banerjee2008model,friedman2008sparse use the latter. For various expositional reasons, we will penalize $\Omega$ instead of $\Omega^{-}$. However in the Monte Carlo simulation, we will also compare with the variant of GLASSO where the diagonals are not penalized.
Although we are comparing SGLASSO with GLASSO in this paper, we could also combine them together. That is, $\hat{\Omega} := \operatornamewithlimits{argmax}_{\Omega\succ0} \big\{\log{\det \Omega} - \mathrm{tr}(\Omega \hat{\Sigma}) - \lambda_{1T}\|\Omega\|_{1,1} - \lambda_{2T}\|\Omega\|^{2}_{1,2} \big\}$, and use cross-validation procedures to tune the parameters $\lambda_{1T}$ and $\lambda_{2T}$.
We do not require that the data-generating process be Gaussian.\footnote{When the DGP is non-Gaussian, the graphical model estimated using GLASSO or SGLASSO is still useful -- it corresponds to the structure of pairwise conditional correlations (partialling out other variables).} Therefore our estimator can be seen as maximizing a Gaussian quasi-likelihood subject to regularization. When there is no penalty ($\lambda_T = 0$), then $\hat{\Omega}$ defined in both ((ref)) and ((ref)) correspond to the Quasi Maximum Likelihood estimator of the inverse covariance matrix, which is given by the inverse of the empirical sample covariance matrix. As we mentioned in the introduction, it is well-known that the unpenalized sample estimator behaves poorly, and it is unsuitable for graphical modeling when zeros of the estimates are important.
We illustrate the difference between the $L_{1,2}$ and the conventional $L_{1,1}$ lasso norms. Define $\Omega_{1}$ and $\Omega_{2}$ as the precision matrices given by equations ((ref)) and ((ref)) below.
The precision matrices $\Omega_{1}$ and $\Omega_{2}$ give rise respectively to the graphical models in Figure (ref) and Figure (ref) below, where links in the graphs represent non-zero entries in the precision matrices.
Both the standard Lasso and Frobenius norms do not distinguish between $\Omega_{1}$ and $\Omega_{2}$, while our proposed $L_{1,2}$-norm does. The reason is clear: both $\Omega_{1}$ and $\Omega_{2}$ have the same number of zeros. They have the same sparsity but the structure of this sparsity is very different. For the $L_{1,1}$ lasso norm, we have: $L_{1,1}(\Omega_{1}) = L_{1,1}(\Omega_{2}) = \sum_{j=1}^{p} \sum_{i \neq j}^{p}|\frac{1}{\sqrt{5}}| + 5 = \frac{8}{\sqrt{5}}+ 5$. Similarly, for the $L_{2,2}$ Frobenius norm, we have $L_{2,2}(\Omega_{1}) = L_{2,2}(\Omega_{2})$.
Now on the other hand, the $L_{1,2}$-norm can distinguish between these two structures. The $L_{1,2}$ norm evaluated at $\Omega_{1}$ for the star graph is $\|\Omega_{1}\|^{2}_{1,2} \approx 16.16$. While the $L_{1,2}$ norm evaluated at $\Omega_{2}$ for the AR(1) graph is $\|\Omega_{2}\|^{2}_{1,2} \approx 14.96$.
In this section we derive the limiting distribution of SGLASSO when the number of samples $T$ goes to infinity while $p$ is fixed (following yuan2007model). The main finding of this section is that SGLASSO is asymptotically equivalent to an infeasible GLASSO problem which prioritizes the sparsity-recovery of high-degree nodes. More precisely, we show that SGLASSO is asymptotically equivalent to a modified GLASSO estimator where each row of $\Omega$ is penalized proportionally according to its {\em true} degree. The true degree of row (or node) $i$ is defined by $d_{i} = \sum_{j =1}^{p}|\Omega_{0,ij}|$. The term “degree" should be understood as “weighted degree" henceforth.
This asymptotic equivalence result has two implications. First, SGLASSO inherits the sparsity-discovery property of GLASSO where zeros in the precision matrix are estimated precisely to be zero (due to the non-differentiability of the penalty function at zeros). Therefore SGLASSO estimates are also sparse. Secondly, SGLASSO differs from GLASSO in that it prioritizes recovering the zero-pattern of high-degree nodes. Nodes with higher weighted degrees can be viewed as being more influential and important, and we might be more interested in uncovering the relationships among influential agents.
To derive the limiting distributions of our SGLASSO estimator, we assume a high level assumption on the weak limit of the sample covariance matrix.
where the notation $vech$ is the half-vectorization operator that takes only the lower-triangular part of a symmetric matrix.
This high level condition holds when $X_{j,t} X_{k,t}$ have higher moments and the serial dependence is weak. When ${\mathbf X}_t \sim_{i.i.d.} \mathcal{N}(0,\Sigma_0)$, then $\sqrt{T} ( \text{vech}\: (\hat{\Sigma} - \Sigma_0 )) \Rightarrow \text{vech} \: (W) \sim \mathcal{N}(0,\Lambda)$, where $\Lambda$ is such that $\mathrm{Cov}(W_{ij}, W_{i'j'}) = \mathrm{Cov}(X_{i,t}X_{j,t}, X_{i',t}X_{j',t})$.
Now let $D_p$ be a $p \times p$ matrix whose $(i,j)$ entry is $d_{j} \equiv \sum_{k=1}^{p} |\Omega_{0,jk}|$, the true weighted degree of node $j$. Define
to be a variant of the GLASSO estimator. Here, the operator $\circ$ is the element-wise multiplication. The SGLASSO estimator is defined as: $\hat{\Omega}_{\lambda} = \operatornamewithlimits{argmax}_{\Omega\succ0} \{\log{\det \Omega} - \mathrm{tr}(\Omega \hat{\Sigma}) - \lambda_{T} \| \Omega \|^{2}_{1,2} \}$.
Compared to the limit in Theorem (ref), the GLASSO estimator has the following limiting distribution (yuan2007model)\footnote{yuan2007model consider the GLASSO estimator where the diagonals are not penalized, so that $V_{\text{glasso}}(U) = \mathrm{tr}(U \Sigma_{0} U \Sigma_{0}) + \mathrm{tr}(U W) + \lambda_{0} \sum_{i=1}^{p} \sum^{p}_{j \neq i} g(u_{ij})$ in their paper.}: \[ \operatornamewithlimits{argmin}_{U \in \mathbb{R}^{p \times p}} V_{\text{glasso}}(U) \] where
where $\sqrt{T}\lambda_T \rightarrow \lambda_0$.
Comparing equations ((ref)) and ((ref)), we see that the SGLASSO obtains the same limiting distribution as that of a modified GLASSO with element-wise LASSO penalty given by $2\| D_p \circ \Omega\|_{1,1}$. Since $D_p$ is constructed using the true, unknown values of $\Omega_{0}$, SGLASSO cannot be simply implemented as GLASSO.
Moreover from Proposition (ref) below, we can further say that SGLASSO is asymptotically equivalent to a variant of GLASSO where the penalty term is $\sum_{i=1}^p \sum_{j=1}^p (d_i + d_j) | \Omega_{ij} |$. That is, each entry $\Omega_{ij}$ is given the penalty factor equals to the sum of the (true) weighted degrees of nodes $i$ and $j$.
All proofs are relegated to the Appendix. In the Supporting Information, we plot and illustrate the asymptotic distribution in Theorem (ref). We show that with a higher probability, SGLASSO correctly estimates a non-link (a zero entry) belonging to a high-degree node.
The distribution of $\operatornamewithlimits{argmin}_{ U \in \mathbb{R}^{p \times p}}V(U)$ can be simulated after replacing the unknown components $\Sigma_0, \Omega_0,$ and $\Lambda$ with their consistent estimates, and this can be used in inference procedures on $\Omega_0$.
A noticeable feature, however, is that the limit distribution in Theorem (ref) is discontinuous with respect to true DGP. This is because $g_{ij}(u)$ is discontinuous as a function of $\Omega_{0,ij}$ at $\Omega_{0,ij}=0$. This implies that the inference based on this limit would only be valid pointwise, not uniformly in $\Omega_0$.\footnote{See belloni2014inference and van2014asymptotically. We thank one of the referees who pointed this out.}
We compare our estimator against GLASSO by considering 3 different graphical models here. In the Supplemental Information, we consider 8 other models. The first two models are depicted in Figure (ref). The last model is the graphical model calibrated to the empirical application as depicted in Figure (ref). From a given graphical model, we generate the true precision matrix $\Omega_{0}$ such that $\Omega_{0,ij} = 0$ if and only if there is a link between nodes $i$ and $j$, otherwise we set $\Omega_{0,ij} = 0.2$. We set $\Omega_{0,ii} = 1$.
For each model, we draw 1,000 independent datasets from $\mathcal{N}(0,\Omega^{-1}_{0})$. That is, each dataset comprises of $({\bf X}_{t})_{t=1}^{T}$, where ${\bf X}_{t} \in \mathbb{R}^{p}$ is randomly drawn from $\mathcal{N}(0,\Omega^{-1}_{0})$. For the sample size, we consider $T=20$ and $T=50$.
Both estimators involve choosing the $\lambda_T$ tuning parameters. We use a 2-folds cross-validation procedure to tune $\lambda_T$. Specifically, we use the Kullback-Leibler (KL) loss averaged over the two-folds to evaluate predictive accuracies. Equation ((ref)) gives the KL loss between the estimated $\hat{\Omega}$ from the training set versus the estimated $\Omega$ from the validation set:
We report the simulation result in Table (ref). In the table, the (a) columns corresponds to SGLASSO whereas the (b) columns refer to GLASSO. In Columns 1(a) and 1(b), we report the optimal $\lambda_T$ as determined by cross-validations, averaged across the 1,000 replications. In Columns 2(a) and 2(b), we report the Kullback-Leibler loss averaged across 1,000 replications.\footnote{The KL loss between $\hat{\Omega}$ and $\Omega_{0}$ is given by $KL(\hat{\Omega},\Omega_{0}) = \log \det(\Omega_{0}) - \log \det \hat{\Omega} + \mathrm{tr}(\hat{\Omega} \Omega_{0}^{-1}) - p$} In Columns 3(a) and 3(b), we report the average Frobenius loss between $\hat{\Omega}$ and $\Omega_{0}$.
In the last two columns of Table (ref), we report the $F_{1}$ score, which measures the accuracy of graph recovery. The Kullback-Leibler loss and the Frobenius norms may not fully capture how accurately the zeros are recovered. We introduce an additional metric: $F_{1} = \frac{2 \text{precision}\cdot \text{recall}}{\text{precision}+\text{recall}}$, where {\em precision} is the ratio of true positives (TP) to all predicted positives (TP + FP), {\em recall} is the ratio of true positives to all actual positives (TP + FN). Alternatively, the $F_{1}$ score can be written as $F_{1} = \frac{2TP}{2TP + FP + FN}$. The $F_{1}$ score measures the quality of a binary classifier by equally balancing both the precision and the recall of a classifier. The larger the $F_{1}$ score is, the better the classifier is. The $F_{1}$ score is commonly used in machine learning to evaluate binary classifiers. For instance, the Yelp competition uses the $F_1$ score as a metric to rank competing models.\footnote{https://www.kaggle.com/c/yelp-restaurant-photo-classification} The $F_{1}$ score is favored over the metric $Accuracy = \frac{TP+TN}{TP + TN + FP + FN}$ especially in our current setting where the graphical models are sparse. This is because a model that naively predicts all negatives will obtain a high {\em Accuracy} score just because there are many actual negatives, and the TN term dominates the {\em Accuracy} score.
The conclusion from Table (ref) is that SGLASSO achieves significantly lower Kullback-Leibler and Frobenius losses. Moreover in in terms of graph accuracy, SGLASSO also outperforms GLASSO, as indicated by the $F_{1}$ scores.
To abstract away from the effects of cross-validations, we consider the lowest Kullback-Leibler losses that can be achieved by our estimator versus the GLASSO. Specifically, we vary $\lambda_{T}$ and at each $\lambda_{T}$, we compute the KL and Frobenius losses between the true $\Omega_{0}$ and $\hat{\Omega}$. We see from Table (ref) that the superior performance of SGLASSO over GLASSO exists after taking away the randomness due to cross-validations. The lowest possible KL losses achievable by our estimator appears to be smaller, than the corresponding KL losses for the GLASSO.
Theorem (ref) says that the SGLASSO prioritizes recovering the sparsity between nodes that have higher degrees. We now show some numerical evidence that lends support to this. For each model, we examined the pair of nodes $(i,j)$ such that $d_{i} + d_{j}$ is largest and $\Omega_{0,ij} = 0$, where $d_{i} = \sum_{k=1}^{p}|\Omega_{0,ik}|$ is the true weighted degree of node $i$. For instance, for model 1, the $\operatornamewithlimits{argmax}_{i,j} (d_{i}+d_{j})$ such that $\Omega_{0,ij} =0$ would correspond to the pairs of nodes $(1,4)$ and $(2,5)$. For simplicity, if there are multiple pairs of nodes that maximize $d_{i}+d_{j}$, we will pick the pair of nodes that comes first when the adjacency matrix is vectorized.
In each of the model (set $T=20$), we calculate the fraction of times that GLASSO and SGLASSO correctly recover $\Omega_{0,ij} =0$ for the largest $d_{i} + d_{j}$. We plot this result in Figure (ref), which shows that at any given $\lambda$, the probability of correctly recovering $\Omega_{0,ij} =0$ for high-degree nodes is greater when SGLASSO is used, compared to GLASSO.
Computing the SGLASSO estimator is a convex optimization problem. More precisely, the problem consists of minimizing a smooth convex term given by $-\log{\det \Omega} + \mathrm{tr}(\Omega \hat{\Sigma})$, and a non-differentiable convex term corresponding to $\lambda_T\|\Omega\|^{2}_{1,2}$. Since the objective can be formulated as minimizing the sum of a differentiable convex function and a non-differentiable convex function, we can use the {\em proximal gradient method} (beck2009fast) to compute SGLASSO.
In this paper, we use CVX, a Matlab package for specifying and solving convex programs (grant2008cvx). CVX supports convex objective functions that are non-smooth (grant2008graph). CVX also allows us to easily enforce symmetry and positive-definiteness of the matrix minimand, which is needed here. The code is readily available from the authors upon request.
For a real-world empirical illustration, we use the classic Grunfeld investment data (see greene2012econometric or baltagi2008econometric). The goal is to estimate the network dependence structure of firms' investment decision.\footnote{The data can be downloaded freely from greene2012econometric's companion website.}
The data consist of time series of $T = 20$ yearly observations on $p=10$ firms. The three variables are $I_{it} = $ real gross investment of firm $i$ in year $t$, $F_{it} = \text{real market value of the firm}$, $C_{it} = \text{real value of capital stock such as plant and equipment}$.
The investment decision for firm $i$ is modeled as follows: $I_{it} = \beta_{1i} + \beta_{2i}F_{it} + \beta_{3i}C_{it} + \epsilon_{it}$, where $(\epsilon_{1t},\epsilon_{2t},\dots,\epsilon_{pt}) \sim iid (0,\Omega_{0}^{-1})$. The zeros and non-zeros in $\Omega_{0}$ correspond to the dependence structure of firms' investment decisions (the object of interest here).
For each firm $i$, we first estimate the linear equation: $I_{it} = \beta_{1i} + \beta_{2i}F_{it} + \beta_{3i}C_{it} + \epsilon_{it}$. This is also the first step in the estimation of Seemingly Unrelated Regressions. Secondly, we use our proposed SGLASSO to estimate the precision matrix of the first-step residuals $\hat{\epsilon}_{it}$, $i=1,\dots,10$, $t=1,\dots,20$. The tuning parameter is determined using a two-fold cross-validation procedure. The result is shown in Figure (ref). (We also estimated the graphical model corresponding to the observed investment variable $I_{it}$, the result is similar)
The recovered graphical model exhibits a clear core-periphery structure. There are two groups of firms: a group of core firms and a group of periphery firms. The core firms are linked to each other, as well as linked to the periphery firms. The periphery firms however, are {\em not} linked to each other, but are only linked to the core firms. Therefore every node is linked directly to a core firm. The core firms here are firms 1, 2 and 3, which are respectively General Motors, U.S. Steel and General Electric. Using the GLASSO estimator, we recover an identical graphical model -- the finding is robust to different choices of penalty functions.
Finally we cannot recover this core-periphery structure using the sample estimator (without regularization). To show this, we estimate the sample precision matrix of the residual estimates $\hat{\epsilon}_{it}$ from the first step. The graphical model corresponding to this sample precision matrix is the complete graph, where every node is linked to all other nodes -- the sample precision matrix is dense and do not contain any zeros.
Graphical model is a potentially useful tool for economists. Using graphical models, we can obtain the network of dependence among random variables. We introduce Structured Graphical Lasso (SGLASSO) as an estimator of graphical models. Using a classic firms' investment dataset, we find that the dependence structure of firms' investment decision exhibits a core-periphery network. This has some relevant policy implications -- small shocks to those core firms will affect the entire network, creating large aggregate fluctuations as in acemoglu2012network. Indeed per our finding, one of the core firms is General Motors, whose Chapter 11 restructuring posed a systemic risk to the U.S. economy in 2009.
We thank Victor Chernozhukov and two referees for helpful comments and suggestions. The first draft of the paper was written while Moon and Chiong were Associate Director and a postdoctoral fellow of USC Dornsife INET, respectively. Moon acknowledges that this work was supported by the Ministry of Education of the Republic of Korea and the National Research Foundation of Korea (NRF-2017S1A5A2A01023679).