Ridge Estimation of High Dimensional Two-Way Fixed Effect Regression

Ridge Estimation of High Dimensional Two-Way Fixed Effect Regression

Introduction

In this paper, we study the two-way fixed-effect model: \[ y_{ijs}=\mu_{i}+\phi_{j}+u_{ijs},i\in\{1,...,n\},j\in\{1,...,p\},s\in\{0,...,d_{ij}\}, \] when the observed matches $(i,j)$ have a sparse bipartite graph structure. In our running example, $y_{ijt}$ is the log of the annualized earnings of worker $i$ employed by firm $j$ in the $s$th occurrence of that particular match, or the residual from a preliminary regression on observed covariates (AbowdKramarzMargolis1999). In most applications, the worker panel covers a limited number of years (5 or 7 years). The network of firm-worker matches is therefore sparse as workers sample a small number of firms (i.e. $d_{ij}=0$ for most $j$) during the observation period. The number of wage observations per worker ($d_{i\cdot}=\sum_{j}d_{ij}$) is of the order of the number of years of observations, sometimes lower because of unemployment, sometimes greater because of infra-year job mobility. Similar two-way fixed effect structures arise in a variety of settings (see LarremoreClausetJacobs2014). For example, if $i$ indexes authors and $j$ indexes journals (or fields), one can model a paper-level outcome such as citations as the sum of an author effect and an outlet effect when authors publish a small number of papers in a fixed time interval. Likewise, if $i$ indexes pollinator species and $j$ indexes plant species, an interaction-level measure such as fruit set (or seed set) can be modeled as the sum of a pollinator “effectiveness” component and a plant “receptivity” component. For ease of exposition, we refer to $y$ as the wage, $i$ as workers, and $j$ as firms in what follows.

Our parameter of interest is the distribution of fixed effects. Obviously, with fixed or bounded degrees no fixed effect can be consistently identified, but what about their distribution? Can one say something on the limit of the distribution of estimated fixed effects when the size of the network becomes large? Obviously, the answer to this question should depend on the estimator.

If the network is sparse, the OLS estimator of the fixed-effect model faces a weak identification problem. JochmansWeidner2019 show that the variance of the OLS estimator depends on the value of the second lowest eigenvalue of the Laplacian, which determines the connectivity of the network. The Cheeger inequality (Chung1997) bounds the second eigenvalue of the Laplacian and the Cheeger constant measures the existence of a “bottleneck”. ZhangRohe2018 further show that sparse, inhomogeneous graphs (edges are drawn independently with node-dependent probabilities) generate a number of bottlenecks (“dangling nodes”) that increases with the number of nodes. The number of very small eigenvalues of the Laplacian increases correspondingly.

This issue is well studied in the community detection literature, aiming at finding subsets of nodes with a higher probability to connect to each other. For sparse graphs, node clustering often returns a partition with one big cluster containing most of the data and many small clusters. ChaudhuriGrahamTsiatas2012 proposed to regularize the graph by making it less sparse by adding many “weak links”. Concretely, this amounts to replacing the adjacency matrix $A$ by $A+\lambda\mathbf{1}_{n}\mathbf{1}_{n}^{\top}$, for a small $\lambda$. Since then, several other ways of regularizing the graph have been proposed. The solution of QinRohe2013 is of particular interest to us because it relates to the ridge regression estimator. They only add a small number to the degree matrix, not to the whole adjacency matrix. Recent work by DallAmicoCouilletTremblay2021 shows that Qin and Rohe's particular form of regularization (as well as the equivalent one Dall'Amico et al. propose) outperforms other competing algorithms. In this paper, we will build on the asymptotic theory developed in the statistical network literature to advance our understanding of the asymptotic properties of the ridge estimator of the two-way fixed effect model with a sparse bipartite network structure.

Three features of the underlying network are important to incorporate in this analysis. First, the network is large (many workers and firms) and so is the number of fixed effects that we seek to estimate. Second, the network is sparse, as in all practical situations only a small fraction of worker-firm matches are observed out of the total number of possible ones. Third, the asymptotic theory will assume a large number of workers $n$ and firms $p$, but $p/n$ is largely fixed and independent of $n$. This implies that the degrees of worker and firm nodes are more or less fixed and independent of the size of the network.

Ridge regression in high-dimensional setups is a well-studied problem in the statistical literature. The literature focuses on the regression model $y_{i}=x_{i}^{\top}\beta+u_{i}$, with $n$ independent observations and where the number of regressors $p$ grows with the sample size with $p/n\rightarrow\gamma$. The vector of regressors is assumed spherical, i.e. $x_{i}=\Sigma^{1/2}z_{i}$ where $z_{i}$ has independent entries. The vector of parameters $\beta$ is random with mean 0 and variance $p^{-1}\alpha^{2}I_{p}$. The ridge regression estimator is then $\widehat{\beta}=\left(X^{\top}X+\lambda I_{p}\right)^{-1}X^{\top}Y$, for $Y=(y_{i})$ and $X^{\top}=(x_{1}...x_{n})$. Under these assumptions, DobribanWager2018 show that the expected predictive risk converges almost surely to a deterministic limit that is minimized by the optimal parameter $\lambda^{*}=n\gamma\alpha^{-2}\sigma^{2}$, where $\sigma^{2}$ is the variance of residuals.\footnote{See Dicker2013,ElKaroui2018,RichardsMourtadaRosasco2020,WuXu2020,BigotDaboMale2024,PatilDuTibshirani2025 for other references.} This setup is rather different from ours. In our setup the elements of $\beta$ do not shrink to zero. At the same time the linear form $x_{i}^{\top}\beta$ does not diverge because of sparsity: the number of nonzero entries in $x_{i}$ increases with $n$ and $p$.

We will need a model for the underlying network. We assume a Degree-Corrected Stochastic Block model (KarrerNewman2011,LarremoreClausetJacobs2014), which has become the standard model in the community detection literature. Specifically, we first assume that the set of workers and firms can be partitioned into $K$ groups (communities). Then, the number $d_{ij}$ of links between $i$ and $j$ is drawn independently from a Binomial$(T,p_{ij})$, where $p_{ij}$ is the probability of a link. This probability is higher if $i$ and $j$ belong to the same community.\footnote{See also Nimczik2017,AbowdMcKinneySchmutte2019,BonhommeLamadonManresa2019,LentzPiyapromdeeRobin2023. These papers usually allow for a Markovian process of job mobility, but a little dependence between $d_{ij}$ and $d_{ij'}$ should not affect the overall stochastic structure of the network.}

For any formed network of worker and firm links, we then generate a two-way fixed effect model with homoscedastic residuals. Let $\beta=(\mu^{\top},\phi^{\top})^{\top}$ denote the vector of parameters. We consider the ridge estimator $\widehat{\beta}$ with $X=(W,F)$, where $W$ is the matrix of worker dummies and $F$ that of firm dummies, using two different regularization parameters, one for $\mu$ and one for $\phi$. Our main theoretical contribution is to derive deterministic equivalents for the ridge estimator’s first- and second-order errors. Under our random graph model, if the ridge penalties are of order $\ln(n+p)$, then the high dimensional bias vector $\mathbb E[\widehat{\beta}-\beta]$ and its variance matrix are close to deterministic limits in operator matrix norm. These deterministic equivalents are computed from the expected network, where each node has the expected degree and where the adjacency matrix is the expected adjacency matrix. This helps explain why the overall distribution of fixed effects can be accurately approximated even though each fixed effect cannot be consistently estimated because the number of observations per worker and firm does not increase with the network size. A further study of these deterministic equivalents using random matrix theory is beyond the scope of this paper.

We run a set of Monte Carlo experiments showing that ridge regularization substantially improves the recovery of the fixed-effect distribution relative to OLS. Moreover, the $\ln(n+p)$ scaling is a sufficient condition; our simulations indicate that tuning ridge penalties using prediction error (cross-validation) tends to yield rather lower ridge parameters. Following the empirical literature, we consider the decomposition of the variance of the outcome variance into the contributions of the variance of the worker effect, the variance of the firm effect, the contribution of the covariance between worker and firm effects, and the variance of residuals. When the network is very sparse, homoscedastic bias correction fails,\footnote{See AndrewsGillSchankUpward2008, Gaure2014, KlineSaggioSoelvsten2020 and AzkarateAskasuaZerecero2022.} while the ridge estimator yields more accurate variance decompositions. We also report estimations on actual matched employer-employee data (the French DADS panel). The log wage variance decomposition derived from OLS estimates using the largest connected component delivers a negative contribution of the fixed effect covariance. We find that homoscedastic bias corrections are minimal. However, the ridge estimator yields a more reasonable variance decomposition. The worker effect explains a more reasonable share of the total variance of log wages and the covariance term is positive instead of being negative. We also consider an estimation with the largest strongest connected component, which allows to use the Leave-One-Out (LOO) bias correction method of KlineSaggioSoelvsten2020. OLS and the homoscedastic bias correction improve, but very little. Ridge does not change much and the LOO bias correction comes close to ridge.

The layout of the paper is as follows. In the next section, we describe the ridge estimator and relate it to the underlying graph of connections between workers and firms. Section (ref) develops a Stochastic Block Model for the network. Asymptotic theory is developed in Section (ref) (the network) and Section (ref) (the ridge estimator). The proofs are in the Appendix.

The two-way fixed effect model

The two-way fixed effect model can be written in matrix form as \[ Y=W\mu+F\phi+U=X\beta+U,\quad X=(W,F),\quad\beta=(\mu^{\top},\phi^{\top})^{\top}. \] The vector $Y=(y_{ijs})\in\mathbb{R}^{N}$ (with $N=\sum_{i,j}d_{ij}$) is a vector of continuous outcomes. The matrix $W\in\{0,1\}^{N\times n}$ associates outcome observations to worker IDs $i\in\{1,...,n\}$ and the matrix $F\in\{0,1\}^{N\times p}$ associates observations to firm IDs $j\in\{1,...,p\}$. We assume that $U$ is iid with mean 0 and variance $\sigma^{2}I_{n}$.

OLS and bipartite graph interpretation

The OLS estimator $\widehat{\beta}=(\widehat{\mu}^{\top},\widehat{\phi}^{\top})^{\top}$ solves the normal equations $X^{\top}X\widehat{\beta}=X^{\top}Y$ with

and \[ \det\left(X^{\top}X\right)=\det\left(D_{w}\right)\det\left(D_{f}-B^{\top}D_{w}^{-1}B\right)=\det\left(D_{f}\right)\det\left(D_{w}-BD_{f}^{-1}B^{\top}\right). \] Each entry of $B=W^{\top}F=(d_{ij})$ is the number of wage observations per match $(i,j)$. It is a weighted adjacency matrix of the bipartite graph connecting workers and firms. We focus on the case where $B$ is sparse: $d_{ij}=0$ for most potential links. We also have \[ D_{w}=W^{\top}W=\text{diag}(d_{i\cdot}),\quad D_{f}=F^{\top}F=\text{diag}(d_{\cdot j}). \] These are the degree matrices of workers and firms: $d_{i\cdot}=\sum_{j}d_{ij}$ is the number of observations for worker $i$, and $d_{\cdot j}=\sum_{i}d_{ij}$ is the number of observations for firm $j$. We assume that worker and firm degrees are all positive, ensuring that $\det\left(D_{w}\right)\ne0$ and $\det\left(D_{f}\right)\ne0$.

Assuming that no worker or firm has zero degree, the matrix \[ \widetilde{A}_{f}=B^{\top}D_{w}^{-1}B=\left(\sum_{i}\frac{d_{ij}d_{ij'}}{d_{i\cdot}}\right)_{j,j'}. \] is a weighted adjacency matrix of the undirected graph connecting firms through common employees, also called a weighted one-mode projection of the original bipartite graph (see ZhouRenMedoZhang2007). This matrix is sparse if for most firm pairs $(j,j')$, no worker is observed at both firms over the observation window. In practice, for the usual lengths of periods of observation (5 or 7 years), this is typically the case. The matrix $\widetilde{L}_{f}:=D_{f}-B^{\top}D_{w}^{-1}B=F^{\top}M_{W}F$, where $M_{W}=I_{N}-W(W^{\top}W)^{-1}W^{\top}$, is the corresponding Laplacian.

A similar interpretation applies to the matrix \[ \widetilde{A}_{w}=BD_{f}^{-1}B^{\top}=\left(\sum_{j}\frac{d_{ij}d_{i'j}}{d_{\cdot j}}\right)_{i,i'}. \] It is the adjacency matrix of the undirected, weighted graph connecting workers through common employers. The matrix $\widetilde{L}_{w}:=D_{w}-BD_{f}^{-1}B^{\top}=W^{\top}M_{f}W$ is the worker Laplacian.

In the sequel, we shall use the following normalized, symmetric versions of these adjacency matrices: \[ A_{w}=\left(\sum_{j}\frac{d_{ij}d_{i'j}}{\sqrt{d_{i\cdot}d_{i'\cdot}}d_{\cdot j}}\right)_{i,i'}=EE^{\top},\quad A_{f}=\left(\sum_{i}\frac{d_{ij}d_{ij'}}{\sqrt{d_{\cdot j}d_{\cdot j'}}d_{i\cdot}}\right)_{j,j'}=E^{\top}E, \] with \[ E=D_{w}^{-1/2}BD_{f}^{-1/2}=\left(\frac{d_{ij}}{\sqrt{d_{i\cdot}d_{\cdot j}}}\right)_{i,j}. \] We also define the normalized versions of the Laplacians: $L_{w}=I_{n}-A_{w}$ and $L_{f}=I_{p}-A_{f}$. The normalized adjacency and the Laplacian matrices satisfy the following properties.

Lemma (ref) bounds the eigenvalues of the Laplacian matrices between 0 and 1. It also asserts that whatever the connectedness of the graph, 0 is always an eigenvalue of the Laplacian. This is a consequence of the fact that the columns of $W$ and $F$ sum to one. In practice, this implies that the two-way fixed effect regression requires a normalization of the fixed effects $\mu$ and $\phi$. The next lemma identifies the multiplicity of the eigenvalue 0 to the number of connected components of the graph.

To calculate the OLS estimator, one parameter normalization is required for each connected component. In order to estimate the two-way fixed effect model, one must first find connected components, using Tarjan's algorithm for example. Estimation is usually performed on the biggest component.

The ridge estimator

In practice, the biggest connected component may be weakly connected. That is to say, the second lowest eigenvalue of the Laplacian matrix may be close to zero. This is a weak identification issue. The following ridge estimator is a natural solution (Hastie2020):

where $\lambda_{w},\lambda_{f}$ are nonnegative parameters. Denote $D_{w,\lambda}=D_{w}+\lambda_{w}I_{n}$ and $D_{f,\lambda}=D_{f}+\lambda_{f}I_{p}$. The estimator is feasible if \[ \det\left(

\right)=\det\left(D_{w,\lambda}\right)\det\left(\widetilde{L}_{f,\lambda}\right)=\det\left(D_{f,\lambda}\right)\det\left(\widetilde{L}_{w,\lambda}\right)\ne0, \] where $\widetilde{L}_{w,\lambda}:=D_{w,\lambda}-BD_{f,\lambda}^{-1}B^{\top}$ and $\widetilde{L}_{f,\lambda}:=D_{f,\lambda}-B^{\top}D_{w,\lambda}^{-1}B$ are the regularized Laplacians. In which case, blockwise inversion yields \[ \left(

\right)=\left(

\right). \] Let $L_{w,\lambda}=I_{n}-D_{w,\lambda}^{-1/2}BD_{f,\lambda}^{-1}B^{\top}D_{w,\lambda}^{-1/2}$ and $L_{f,\lambda}=I_{p}-D_{f,\lambda}^{-1/2}B^{\top}D_{w,\lambda}^{-1}BD_{f,\lambda}^{-1/2}$ denote normalized versions of the regularized Laplacians. Then, $\widetilde{L}_{w,\lambda}=D_{w,\lambda}^{1/2}L_{w,\lambda}D_{w,\lambda}^{1/2}$ and $\widetilde{L}_{f,\lambda}=D_{f,\lambda}^{1/2}L_{f,\lambda}D_{f,\lambda}^{1/2}$.

The following lemma shows that the ridge estimator is always feasible (even without normalizing the fixed effect parameters) by proving that the eigenvalues of the regularized Laplacian are bounded away from zero. Let $E_{\lambda}=D_{w,\lambda}^{-1/2}BD_{f,\lambda}^{-1/2}.$ The regularized firm Laplacian can be written as $L_{f,\lambda}=I_{p}-E_{\lambda}^{\top}E_{\lambda}$, and the regularized worker Laplacian as $L_{w,\lambda}=I_{n}-E_{\lambda}E_{\lambda}^{\top}$. Hereafter, the matrix norm is the spectral matrix norm ($\left\Vert A\right\Vert =\sqrt{\operatorname{eigmax}(A^{\top}A)}$ for any matrix $A$).

In-sample bias and variance

Suppose that $U$ is independent of $X=(W,F)$, has mean 0 and $\operatorname{\mathbb{V}} U=\sigma^{2}I_{N}$. Then, the ridge estimator is biased with

The variance of $(\widehat{\mu},\widehat{\phi})$ follows with

A stochastic block model

To study the statistical properties of these estimators, we need a model for the graph of worker-firm matches. We model the bipartite graph as a Degree-Corrected Stochastic Block Model (DasguptaHopcroftMcSherry2004, KarrerNewman2011, LarremoreClausetJacobs2014, RazaeeAminiLi2019).

\paragraph*{Node types}

We first divide an arbitrary number of $n$ workers and $p$ firms into $K$ groups. Let $k\in\{1,...,K\}$ index worker types and $\ell\in\{1,...,K\}$ index firm types. We assume the same number of groups for simplicity. Workers initially draw a type $k_{i}=k$ independently with probability $\pi_{w}(k)$, and firms draw a type $\ell_{j}=\ell$ with probability $\pi_{f}(\ell)$. Let $n_{k}$ denote the number of workers of type $k$. Second, we introduce a parameter $\theta_{j}$ to control the expected degree of firm $j$ with the restriction that the $\theta_{j}'s$ sum to one within firm groups: $\sum_{j'}\theta_{j'}\delta_{\ell_{j'}\ell}=1$. Parameter $\theta_{j}$ can therefore be interpreted as the probability of drawing firm $j$ in its group $\ell$ conditional on drawing any firm from this group. For workers, we assume uniform sampling given type (i.e. with probability $1/n_{k_{i}}$). We could introduce another fixed effect $\theta$ to control worker sampling given type. By using a degree correction for one node type and not for the other type, we cover both the standard SBM and its degree-corrected extension.

\paragraph*{Network}

The block-structure of the network is governed by the affinity matrix $C$, with $C(k,\ell)\ge0$. Independently for each worker-firm couple $(i,j)$, we draw $d_{ij}$ edges between $i$ and $j$ from a Bernoulli distribution where the probability of an edge is $p_{ij}=\frac{1}{n_{k_{i}}}\theta_{j}C(k_{i},\ell_{j})$. The expected number of links between workers of type $k$ and firms of type $\ell$ is \[ \operatorname{\mathbb{E}}\left(\sum_{i,j}d_{ij}\delta_{k_{i}k}\delta_{\ell_{j}\ell}\right)=\sum_{i,j}p_{ij}\delta_{k_{i}k}\delta_{\ell_{j}\ell}=C(k,\ell), \] as $\sum_{k}\delta_{k_{i}k}=n_{k_{i}}$ and $\sum_{j}\delta_{\ell_{j}\ell}\theta_{j}=1$. The expected firm degree is \[ \operatorname{\mathbb{E}}\left(d_{\cdot j}\right)=\operatorname{\mathbb{E}}\left(\sum_{i}d_{ij}\right)=\sum_{i}p_{ij}=\sum_{i}\frac{1}{n_{k_{i}}}\theta_{j}C(k_{i},\ell_{j})=\theta_{j}C(\cdot,\ell_{j}), \] where $C(\cdot,\ell)=\sum_{k}C(k,\ell)$, and the expected worker degree is \[ \operatorname{\mathbb{E}}\left(d_{i\cdot}\right)=\operatorname{\mathbb{E}}\left(\sum_{j}d_{ij}\right)=\sum_{j}p_{ij}=\frac{1}{n_{k_{i}}^{0}}\sum_{\ell}\sum_{j}\theta_{j}\delta_{\ell_{j}\ell}C(k_{i},\ell)=\frac{1}{n_{k_{i}}}C(k_{i},\cdot), \] where $C(k,\cdot)=\sum_{\ell}C(k,\ell)$.

For example, in the simulations we use \[ C=c\frac{p}{K}\left[I_{K}+\delta(J_{K}-I_{K})\right], \] where $c$ is a constant, $J_{K}=\mathbf{1}_{K}\mathbf{1}_{K}^{\top}$ and $\delta$ is a tuning parameter: $\delta=0$ means perfect segregation and $\delta=1$ means indifference. With this specification, the number of connections per pair $(k,\ell)$ is proportional to the number of firms. If workers are uniformly assigned to groups, then $n_{k}\simeq n/K$ and the average worker degree is proportional to $p/n=\gamma$. It is therefore fixed and does not grow with $p$ and $n$ if the numbers of workers and firms grow at the same speed.

The realized graph typically features one main connected component containing orders of magnitude more observations than any other components. For identification reasons, it is usual to use it as estimation sample. Then, $n$ and $p$ denote the numbers of workers and firms in the estimation sample, and $n_{k}$ is the number of workers of type $k$ in the sample.

Asymptotic theory for the network

We develop the asymptotic theory for a large network with many firms and workers. Worker and firm types are given, as well as parameters $\theta_{j}$ and $C$. The edge weights $d_{ij}=\text{Bernoulli}(p_{ij})$, with $p_{ij}=\frac{1}{n_{k_{i}}}\theta_{j}C(k_{i},\ell_{j})$, are the only stochastic variables. The following analysis does not assume the selection of the biggest connected component.

Deterministic equivalent network

Following QinRohe2013, we use a calligraphic font to indicate expectations with respect to $B=(d_{ij})$. We thus denote

where $Z_{w}=\left(\delta_{k_{i}k}\right)\in\{0,1\}^{n\times K}$ and $Z_{f}=\left(\delta_{\ell_{j}\ell}\right)\in\{0,1\}^{p\times K}$ are the selection matrices containing the information on the community memberships of workers and firms. Moreover,

Using the normalization $\sum_{j}\theta_{j}\delta_{\ell_{j}\ell}=1$ and the notation $C(k,\cdot)=\sum_{\ell}C(k,\ell)$, we have

And using the notation $C(\cdot,\ell)=\sum_{k}C(k,\ell)$ and the fact that \[ \sum_{i}\tfrac{1}{n_{k_{i}}}C(k_{i},\ell_{j})=\sum_{k}\sum_{i}\frac{\delta_{k_{i}k}}{n_{k_{i}}}C(k,\ell_{j})=\sum_{k}C(k,\ell_{j}), \] we have

Let $\mathfrak{E}_{\lambda}=\mathfrak{D}_{w,\lambda}^{-1/2}\mathfrak{B}\mathfrak{D}_{f,\lambda}^{-1/2}.$ The generic element of $\mathfrak{E}_{\lambda}$ is

denoting $\widetilde{C}(k,\ell)=\frac{C(k,\ell)}{\sqrt{C(k,\cdot)C(\cdot,\ell)}}$, $\omega_{\lambda,k}=\frac{1}{n_{k}}\sqrt{\frac{C(k,\cdot)}{\frac{1}{n_{k_{i}}}C(k,\cdot)+\lambda_{w}}}$ and $\phi_{\lambda,j}=\theta_{j}\sqrt{\frac{C(\cdot,\ell_{j})}{\theta_{j}C(\cdot,\ell_{j})+\lambda_{f}}}$. The regularization reduces the noise impact of worker communities with few employment links per capita ($\frac{1}{n_{k}}C(k,\cdot)$ small) and firm groups with few employees per firm ($\theta_{j}C(\cdot,\ell_{j})$ small).

Then, define the regularized adjacency and Laplacian matrices

The following lemma shows that the matrix $\mathfrak{E}_{\lambda}$ contains the community structure of the random graph and that the analog Laplacian $\mathfrak{L}_{f,\lambda}=I_{p}-\mathfrak{E}_{\lambda}^{\top}\mathfrak{E}_{\lambda}$ has its eigenvalues bounded away from 0.

Asymptotic bounds for the Laplacian

The main tool used to derive asymptotic bounds is a version of Bernstein inequality for random matrices (see ChungRadcliffe2011 and Tropp2015, Theorem 6.1.1). Let $S_{1},...,S_{n}$ be independent, centered random matrices with common dimension $d_{1}\times d_{2}$ and assume that each one is uniformly bounded: \[ \operatorname{\mathbb{E}} S_{k}=0\quad\text{and}\quad\left\Vert S_{k}\right\Vert \le L,\quad\forall k=1,...,n. \] Introduce the sum $Z=\sum_{k=1}^{n}S_{k}$ and let $v(Z)$ denote the matrix variance statistic of the sum: \[ v(Z)=\max\left\{ \left\Vert \operatorname{\mathbb{E}}\left(ZZ^{\top}\right)\right\Vert ,\left\Vert \operatorname{\mathbb{E}}\left(Z^{\top}Z\right)\right\Vert \right\} =\max\left\{ \left\Vert \sum_{k=1}^{n}\operatorname{\mathbb{E}}\left(S_{k}S_{k}^{\top}\right)\right\Vert ,\left\Vert \sum_{k=1}^{n}\operatorname{\mathbb{E}}\left(S_{k}^{\top}S_{k}\right)\right\Vert \right\} . \] Then, \[ \Pr\left\{ \left\Vert Z\right\Vert \ge t\right\} \le(d_{1}+d_{2})\exp\left(-{\frac{1}{2}t^{2}}/{\left(v(Z)+\frac{1}{3}Lt\right)}\right),\quad\forall t\ge0. \] Furthermore, \[ \operatorname{\mathbb{E}} Z\le\sqrt{2v(Z)\ln(d_{1}+d_{2})}+\frac{1}{3}L\ln(d_{1}+d_{2}). \]

We first prove the following concentration inequality that adapts Theorem 4.1 of QinRohe2013 to regularized Laplacians of bipartite graphs. The ridge coefficients $\lambda_{w}$ and $\lambda_{f}$ allow us to choose $M_{w}\vee M_{f}:=\max(M_{w},M_{f})$ large enough to maintain the wedge $\left\Vert E_{\lambda}-\mathfrak{E}_{\lambda}\right\Vert $ less than any chosen value with any specified probability.

The next theorem shows that the previous property passes to the inverse.

We end this section by showing a similar concentration inequality for the inverse of the un-normalized Laplacian matrices,

The ridge regression

This section shows that the concentration inequalities of the preceding section guarantee that the ridge fixed effects converge to well-defined limits, both in terms of bias and variance.

Rewrite the ridge estimator as

using $Y=W\mu+F\phi+U$, where $U$ has independent entries with mean 0 and variance $\sigma^{2}$.

We assume that the parameters $(\mu,\phi)$ are random with \[ \mu=Z_{w}\mu^{*}+U_{w},\quad\phi=Z_{f}\phi^{*}+U_{f}, \] where $Z_{w}$ and $Z_{f}$ are the $n\times K$ and $p\times K$ matrices indicating worker and firm communities, $\mu^{*}$ and $\phi^{*}$ are $K$-vectors of group fixed effects, and $U_{w}$ and $U_{f}$ are vectors of independent components with mean zero and variance $\sigma_{w}^{2}$ and $\sigma_{f}^{2}$. It is further assumed that $U_{w}$ and $U_{f}$ are independent of $U$.

Bias

Given the community structure $Z=(Z_{w},Z_{f})$ and $\beta^{*}=(\mu^{*},\phi^{*})$, and given the network structure $X=(W,F)$, the biases on $\widehat{\mu}$ and $\widehat{\phi}$ are

Define the deterministic bias limits on $\mu$ and $\phi$ as the ones obtained using the expected adjacency matrix $\mathfrak{B}$,

The following theorem proves that $\mathfrak{b}_{\mu,\lambda}$ and $\mathfrak{b}_{\phi,\lambda}$ are asymptotically good predictors of the bias $b_{\mu,\lambda}$ and $b_{\phi,\lambda}$ in MSE respectively, as long as $\lambda_{w},\lambda_{f}$ grow faster than $\ln(n+p)$.

Variance

The conditional variance of $\widehat{\mu}-\mu$ is

Define the deterministic equivalent of the variance matrix as \[ \mathfrak{V}_{w,\lambda}=\sigma^{2}\widetilde{\mathfrak{L}}_{w,\lambda}^{-1}+\widetilde{\mathfrak{L}}_{w,\lambda}^{-1}\left[\left(\lambda_{w}^{2}\sigma_{w}^{2}-\lambda_{w}\sigma^{2}\right)I_{n}+\left(\lambda_{f}^{2}\sigma_{f}^{2}-\lambda_{f}\sigma^{2}\right)\mathfrak{B}\mathfrak{D}_{f,\lambda}^{-2}\mathfrak{B}^{\top}\right]\widetilde{\mathfrak{L}}_{w,\lambda}^{-1}. \] The firm-side analog of these variances are

The next theorem shows that $\mathfrak{V}_{w,\lambda}$ and $\mathfrak{V}_{f,\lambda}$ are asymptotically equivalent to the variances $V_{w,\lambda}$ and $V_{f,\lambda}$ as long as the ridge parameters grow at the same rate, faster than $\ln(n+p)$.

Prediction

Suppose that we draw another network and outcomes from the same DGP, say $\widetilde{X}=(\widetilde{W},\widetilde{F})$ and $\widetilde{Y}=\widetilde{W}\mu+\widetilde{F}\phi+\widetilde{U}$. The prediction $SSE$ is

as the expectations of $\widetilde{W}^{\top}\widetilde{W}$, $\widetilde{F}^{\top}\widetilde{F}$ and $\widetilde{W}^{\top}\widetilde{F}$ are the same as before.

The first non trivial term is

The other terms follow similarly, yielding \[ SSE=b_{\mu,\lambda}^{\top}\mathfrak{D}_{w}b_{\mu,\lambda}+\text{tr}\left[\mathfrak{D}_{w}V_{w,\lambda}\right]+2b_{\mu,\lambda}^{\top}\mathfrak{B}b_{\phi,\lambda}+\text{tr}\left[\mathfrak{B}C_{\lambda}\right]+b_{\phi,\lambda}^{\top}\mathfrak{D}_{f}b_{\phi,\lambda}+\text{tr}\left[\mathfrak{D}_{f}V_{f,\lambda}\right], \] where $C_{\lambda}=\operatorname{\mathbb{E}}\left[\left(\widehat{\mu}-\mu-b_{\mu,\lambda}\right)\left(\widehat{\phi}-\phi-b_{\phi,\lambda}\right)^{\top}\mid X,Z,\beta^{*}\right]$ denotes the covariance matrix. A deterministic equivalent is obtained by replacing biases and variances by their deterministic equivalents.

Consider the first term: \[ \left|b_{\mu,\lambda}^{\top}\mathfrak{D}_{w}b_{\mu,\lambda}-\mathfrak{b}_{\mu,\lambda}^{\top}\mathfrak{D}_{w}\mathfrak{b}_{\mu,\lambda}\right|\le\left\Vert \mathfrak{D}_{w}\right\Vert \left(\left\Vert b_{\mu,\lambda}\right\Vert +\left\Vert \mathfrak{b}_{\mu,\lambda}\right\Vert \right)\left\Vert b_{\mu,\lambda}-\mathfrak{b}_{\mu,\lambda}\right\Vert , \] where $\left\Vert \mathfrak{D}_{w}\right\Vert =\bar{\delta}_{w}$ is the maximum worker expected degree. It is easy to show that under the condition of the two previous theorems, if the ridge parameters go to infinity at the same rate faster that $\ln(n+p)$, this quantity is bounded of order $nt$.

Consider the second term: \[ \left|\text{tr}\left[\mathfrak{D}_{w}V_{w,\lambda}\right]-\text{tr}\left[\mathfrak{D}_{w}\mathfrak{V}_{w,\lambda}\right]\right|=\left|\text{tr}\left[\mathfrak{D}_{w}\left(V_{w,\lambda}-\mathfrak{V}_{w,\lambda}\right)\right]\right|\le\delta_{w}n\left\Vert V_{w,\lambda}-\mathfrak{V}_{w,\lambda}\right\Vert , \] which is of the same oder as the first term.

And so on for all terms. This shows that the Mean Square Error ($SSE$ divided by number of observations $N$) is close to $n$ times the average worker degree, and will converge to a deterministic equivalent limit. Studying this limit and deriving optimal ridge parameters for high dimensional datasets is outside the scope of this paper.

Applications

Simulation results

We start with an initial number of $n_0$ workers and $p_0$ firms. We simulate the model with \[ C=c\frac{p_0}{K}\left[I_{K}+\delta(J_{K}-I_{K})\right], \] where $c$ is a constant, $J_{K}=\mathbf{1}_{K}\mathbf{1}_{K}^{\top}$ and $\delta$ is a tuning parameter: $\delta=0$ means perfect segregation and $\delta=1$ means indifference. The initial numbers of nodes are set such that $n_{0}=3p_{0}$ and the number of communities is $K=5$. The community strength is controlled by $\delta=0.1$. We draw $\theta_{j}$ independently from a Pareto distribution with scale parameter $\alpha=2$ and minimum value $\theta_{min}=1$: $\Pr\{\theta_{j}>x\}=(\theta_{min}/x)^{\alpha}$. Then, separately for each group $\ell$, we normalize $\theta_{j}$ for all firms $j$ belonging to group $\ell$ by dividing $\theta_{j}$ by $\sum_{j'}\theta_{j'}\delta_{\ell_{j'}\ell}$. We select the biggest connected component and $n$ and $p$ denote the numbers of workers and firms that remain.

Lastly, we consider the following choices of wage parameters:

This implies a form of positive sorting on wages (wage fixed effects are positively correlated). The residual variance is $\sigma=2$.

Table (ref) displays the main characteristics of the simulated networks, one first time with $c=1$, a second time with $c=2$. When $c=1$, the network is very sparse with many disconnected components (after removing nodes with no connection). As $c$ increases, the network becomes denser. We therefore use a much greater initial value of firm nodes $p_{0}$ for $c=1$ than for $c=2$ in order to generate a connected graph of roughly the same size $N$ close to $10000$ observations.

We show in Table (ref) the corresponding simulations for the log-wage variance decomposition. The first row (“true”) shows the true variance decomposition (and the correlation between fixed effects). The fixed effect parameters were chosen to deliver a decomposition that looks like usual empirical ones, with a large worker contribution, a large residual variance and limited firm and sorting contributions. Then, we show the variance decomposition obtained with the OLS estimates of the fixed effects, and the debiased variance components (AndrewsGillSchankUpward2008).\footnote{We do not consider other bias correction methods (Gaure2014, KlineSaggioSoelvsten2020, AzkarateAskasuaZerecero2022) because the model is homoscedastic.} The row labelled “ridge” shows the variance decomposition that is obtained using the ridge estimator. We display the shares of $\operatorname{\mathbb{V}}(W\widehat{\mu})$, $\operatorname{\mathbb{V}}(F\widehat{\phi})$, $2\operatorname{Cov}(W\widehat{\mu},F\widehat{\phi})$. The sixth column of Table (ref) shows the out-of-sample MSE. Lastly, we report the ridge regularization parameters obtained by cross-validation. For cross-validation, we simulate a test sample as follows. We use the same nodes and communities, but draw a new $d_{ij}$ for each couple of worker-firm nodes. Then, we draw a new residual $u_{ijs}$. Again, we keep the largest connected component for prediction.

The OLS estimation is always strongly biased, but the bias correction works very well as long as the graph is not too sparse (i.e. for $c=2$ but not for $c=1$). The variance decomposition using the ridge estimator is a lot less affected by sparsity than OLS. There is a tendency of the ridge estimator to underestimate the shares of $\operatorname{\mathbb{V}}(W\widehat{\mu})$, $\operatorname{\mathbb{V}}(F\widehat{\phi})$, $2\operatorname{Cov}(W\widehat{\mu},F\widehat{\phi})$ and overestimate the “residual share”. We report the ridge regularization parameters as a fraction of the average worker and firm degrees.

Next, we turn to the distribution of fixed effects. Figures (ref) and (ref) show the densities of true and estimated fixed effects and their scatterplots. The scatterplots show the estimated versus the true values. The black line is the 45 degree line for a perfect estimation. We see that OLS estimates are considerably more dispersed than the true distribution, while the ridge estimator is only slightly more concentrated. The scatterplots show evidence of a pattern of the bias for ridge. Large values of the true fixed effects are associated with a stronger negative bias, and vice versa.

Real data application

We then propose an application with French matched employer-employee data (DADS Panel, 1995-2001) similar to the data used in AbowdKramarzMargolis1999. We keep only non agricultural, salaried, private sector employees working full time, younger than 55. We keep all employment spells, including those recorded for a fraction of a year. However, we drop matches with firms with fewer than 50 worker/wage observations and workers with less that 5 recorded spells. The wage concept is earnings per day, and we trim wages below the 0.1 percentile and above the 99.9th one. We regress log wages on job tenure, age and age squared, and use the residuals for our analysis.

After selecting the biggest connected component (out of 2663 ones), we end up with a network with $n=94491$ workers, $p=3408$ firms and $N=450125$ wage observations. The average worker degree is $\overline{d}_{w}=4.7637$ and the average firm degree is $\overline{d}_{f}=132.08$. The average number of wage records per match is 3.5903. The network is very sparse with 127457 matches, which is 0.03958% of the $np$ potential ones.

The log wage residual has a variance of 0.2045 which can be decomposed as reported in Table (ref). The variance decomposition obtained using the OLS estimator of the worker and firm effects yields a rather large negative contribution of the fixed effects covariance to wages. Biases on plugin variance contributions (assuming homoscedasticity) are estimated to be small.

Then we show the variance decomposition obtained from ridge regression.\footnote{We show the contributions to the total variance of the worker effect, the firm effect, 2 times the covariance of worker and firm effects, and a residual contribution that is the sum of the residual variance and the covariances between residuals and fixed effects. Contrary to the OLS estimator, the ridge estimator does not produces orthogonal residuals.} We estimate the ridge model on the years 1995-2000 and keep the last year 2001 for cross-validation. By minimizing out-of-sample error, we tune the ridge parameters $\lambda_{w}$ and $\lambda_{f}$ to be 0.055 and 0.060 times the worker and firm degrees.\footnote{The previous $\log(n+p)$ lower bound is a sufficient condition for our concentration arguments; we do not claim cross-validation yields penalties with this scaling. Instead, we report selected penalties as fractions of average degrees and assess their performance in simulations.} We now estimate a positive covariance between worker and firm effects, and with more reasonable contributions of the worker effects and of the regression residuals.

Next, because we worry about residual heteroscedasticity and autocorrelation, we aim to compare ridge to the Leave-One-Out (LOO) bias correction method of KlineSaggioSoelvsten2020. The LOO method requires firms to remain connected after taking one observation out. We use the strongest version leaving one match out, which asymptotically removes second-order biases when residuals are heteroscedastic or autoregressive within matches. OLS and the homoscedastic bias correction are still very close and deliver a negative covariance term, but of smaller magnitude, proof that the LOO sample is less sparse. However, ridge seems little affected by the choice of the sample. The LOO correction delivers a covariance contribution that is very close to zero, thus closer to ridge. These estimates are also very close to those in BabetGodechotPalladino2022 and AzkarateAskasuaZerecero2022, who use exhaustive DADS data after 2002 (52 million observations) instead of the Panel before 2001, and hourly wages instead of daily wages. Babet et al. use the LOO correction and Azkarate-Zerecero develop a related parametric bootstrap correction.

Next, we show in Figure (ref) the distributions of the worker and firm fixed effects. The ridge estimator is a lot more concentrated than the OLS estimator. Figure (ref) shows each ridge estimate as a function of the corresponding OLS one. Estimates of firm fixed effects tend to disagree more than estimates of worker effects.

Conclusion

In this paper, we have studied ridge penalization as an alternative to OLS for the estimation of the fixed effects in a two-way fixed effect regression model like a log wage equation for matched employer-employee data. We developed a Degree Corrected Stochastic Block model for the creation of the network. The outcome variable was generated given the network in a second stage. We derived the asymptotic properties of the estimator of the fixed effects. We were able to find deterministic equivalents of the bias and the variance of the estimator. In the future, further study of these deterministic equivalents will allow us to characterize the optimal ridge parameter.