The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
116,324 characters
Fixed-$T$ Dynamic Spatial Panel Model with Common Shocks
\begin{titlepage}
\maketitle
\begin{abstract}
We study a dynamic spatial panel model with observed regressors, interactive effects,
and contemporaneous and lagged dependence in a large-$N$, fixed-$T$ framework.
The spatial model constitutes an $N$-dimensional simultaneous-equations system.
In this $N$-equation view, the interactive effects introduce $N$ unit-specific loading vectors.
Estimating them individually when $T$ is fixed
creates the type of incidental-parameters problem underlying Nickell bias.
We use the $N$-equation view to account for spatial simultaneity through the spatial Jacobian.
Crucially, however, we view the same model as a $T$-equation system with $N$ observations.
Together with a spatially enriched random-loadings (SERL) specification,
this $T$-equation view yields a quasi-likelihood without unit-specific incidental parameters.
We propose a computationally tractable block-coordinate algorithm.
Simulations show small estimation errors and generally near-nominal coverage probabilities.
Applying the method to U.S. county female labor-force participation,
we find substantial dynamic and spatial dependence
and an important role for education in explaining the rise in female labor-force participation.
\end{abstract}
\noindent Keywords: interactive effects, incidental parameters, Nickell bias
\thispagestyle{empty}
\end{titlepage}
\setcounter{page}{1}
\section{Introduction}
Dynamic panel data often exhibit two distinct forms of dependence.
One is dependence over time, arising through lagged responses and state persistence.
The other is dependence across units, arising through economic, geographic, or network interactions.
In many applications, these two forms of dependence coexist with common shocks that affect all units but with heterogeneous intensities.
A useful empirical framework should therefore accommodate simultaneously dynamic persistence,
spatial interaction, and interactive effects.
This paper studies such a framework in a large-$N$, fixed-$T$ setting.
We consider a dynamic spatial panel model with a contemporaneous spatial lag,
a lagged dependent variable, observed regressors, and an interactive factor structure.
The specification brings together two familiar sources of cross-sectional dependence.
Weak cross-sectional dependence is captured by the spatial autoregressive structure,
while strong cross-sectional dependence is captured by latent common shocks with heterogeneous loadings.
The parameters of primary interest are the spatial coefficient $\rho$, the dynamic coefficient $\phi$,
and the slope coefficients $\beta$.
The econometric challenge is to conduct likelihood-based inference when the time dimension is fixed,
the cross section is large,
and the model contains both simultaneity across units and dynamic feedback over time.
The framework also allows the observed regressors to be correlated with unit heterogeneity and,
more generally, with the latent factor structure through the factor loadings and common shocks.
The model lies at the intersection of several established literatures.
It builds on the foundational spatial-econometrics framework of \textcite{anselin-1988}
and the spatial and dynamic spatial panel models studied by \textcite{lee-2004}, \textcite{yu-2008}, and \textcite{lee-yu-2010},
among others; see also \textcite{baltagi-2021} and \textcite{elhorst-2014}.
Its common-shock component is related to the interactive
effects literature following \textcite{bai-2009}. Within the likelihood
literature, the model is especially close in spirit to \textcite{shi-lee-2017}
and \textcite{bai-li-2021}, who also combine spatial dependence, dynamics, and
latent common components. Those papers adopt large-$N$, large-$T$ asymptotics
and estimate a growing collection of incidental parameters; the resulting QML
estimators require bias correction. The present paper instead focuses on large
$N$ and fixed $T$.
Fixed $T$ makes the treatment of individual heterogeneity central.
In the present model,
the factor loadings play a role analogous to that of unit fixed effects in conventional dynamic panel models.
Eliminating unit effects by a within transformation in the latter setting
produces the fixed-$T$ bias highlighted by \textcite{nickell-1981}.
The analogous challenge here is to avoid estimating the factor loadings unit by unit.
Our central insight is to exploit two complementary representations of the same
model. Viewed cross-sectionally, the spatial model constitutes an
$N$-dimensional simultaneous-equations system. The contemporaneous spatial
transformation is governed by $B(\rho)=I_N-\rho W$, so the conditional
likelihood retains the familiar spatial Jacobian $|B(\rho)|^T$. Treating the
interactive effects directly in this $N$-equation representation, however,
would require estimating $N$ unit-specific loading vectors and would therefore
reintroduce a fixed-$T$ incidental-parameters problem.
Crucially, following the perspective of \textcite{bai-2013,bai-2024},
we also view the same model as a $T$-equation system with $N$ observations,
where each unit contributes a $T$-dimensional time series. Conditional on the initial
observation, the temporal transformation is represented by the lower-triangular
operator $R_T(\phi)$ with unit diagonal and therefore has determinant one; it
contributes no additional parameter-dependent Jacobian term.
The change in representation does not by itself eliminate the loading heterogeneity.
To avoid estimating the loading vectors unit by unit,
we introduce a spatially enriched random-loadings (SERL) specification.
SERL projects the factor loadings on observed controls,
the initial condition,
and their spatial transformations,
thereby separating spatially propagated systematic heterogeneity from a conditionally idiosyncratic residual component.
Combining SERL with the $T$-equation representation replaces the individual loading vectors with a finite-dimensional conditional mean and a low-dimensional covariance matrix.
The resulting quasi-likelihood retains the spatial Jacobian while avoiding unit-specific incidental parameters.
Two closest fixed-$T$ contributions provide useful contrasts.
\textcite{li-yang-2021} study a dynamic spatial panel with correlated random
effects.
Their scalar, time-invariant correlated-random-effects equation is a special case of the random-loading projection underlying SERL.
Their Mundlak specification accommodates dependence between the
regressors and individual heterogeneity, but conditioning on an initial outcome
generated by the dynamic process generally leaves the conditional score
uncentered. They therefore adjust the score using information about when the
process began. They also show that spatial dependence invalidates covariance
estimation based only on raw individual-score outer products and develop a
spatial decomposition of the pooled score. \textcite{li-miao-yang-2026}, by
contrast, extend the fixed-$T$ analysis to interactive fixed effects, treating
the individual loadings as parameters and using adjusted M-estimation together
with a degrees-of-freedom correction for inference. Our approach draws on the
score decomposition of \textcite{li-yang-2021} for robust covariance estimation
but differs from both methods in its treatment of individual heterogeneity and
the initial condition. We project the factor loadings on a finite-dimensional
collection of observed controls, including spatial transformations of the
initial condition and regressor summaries.
Our approach therefore does not estimate the unit-specific fixed effects or,
more generally, the individual factor loadings. We develop a conditional
Gaussian quasi-maximum likelihood estimator under large-$N$, fixed-$T$
asymptotics. The estimator combines the spatial Jacobian with the
low-rank-plus-diagonal covariance structure induced by the interactive effects.
The resulting criterion admits a computationally tractable block-coordinate
implementation. Conditional on the spatial parameter, the dynamic and slope
coefficients are updated by GLS, while the factor and covariance components are
updated using the low-rank structure. The spatial parameter is then updated
through a one-dimensional conditional maximization. This profiling strategy
remains computationally feasible even when the cross section is large and
extends naturally to the additional specifications considered below. Adding
$Wy_{t-1}$ only enlarges the GLS design, whereas adding a spatial-error filter
replaces the scalar outer search with a two-dimensional optimization over the
roots of a second-order spatial polynomial.
We provide an implementation of the estimation and inference procedures
in the \href{https://github.com/jessekelighine/dspserl}{\texttt{dspserl}} R package.
The simulations examine both estimation and inference. In the baseline designs,
estimation error is small and coverage is close to nominal for the regression
slopes, although first-order Wald intervals for the spatial and dynamic
parameters under-cover in some short, persistent panels. A second experiment
generates the initial outcome from a long-running process and shows that the
first spatial lag in the loading projection produces most of the reduction in
bias and improvement in coverage.
In the empirical application to county-level female labor-force participation data,
the preferred second-order specification leaves the education coefficient positive and
economically important. This result is close to the original OLS evidence of
\textcite{fogli-veldkamp-2011} and provides an alternative interpretation to the
decomposition in \textcite{tziolas-elhorst-2023}.
The remainder of the paper is organized as follows.
Section~\ref{sec:model} presents the model, introduces the two Jacobian arguments, and develops the conditional likelihood.
Section~\ref{sec:asymptotics} establishes the large-$N$, fixed-$T$ inferential theory.
Section~\ref{sec:estimation_algorithm} describes the block coordinate estimation algorithm and the associated computational details.
Section~\ref{sec:extensions} develops extensions with a lagged spatial outcome
and spatially correlated errors.
Section~\ref{sec:simulation} reports Monte Carlo evidence,
Section~\ref{sec:application} applies the method to county-level female labor
force participation, and
Section~\ref{sec:conclusion} concludes.
\section{Dynamic Spatial Panel Model with Common Shocks}\label{sec:model}
This section develops the model and likelihood under fixed $T$.
The key features are spatial dependence, dynamic feedback, and interactive effects.
\subsection{Model}
Let $i=1,...,N$ index units and $t=1,...,T$ index time.
For any positive integer $m$,
let $\mathbf{1}_{m}$ denote the $m$-dimensional vector of ones.
Let $W$ be a known $N\times N$ spatial weights matrix.
Typically, the diagonal of $W$ is assumed to be zero, i.e., $w_{ii}=0$
\parencite{anselin-1988,elhorst-2014}.
Consider the dynamic spatial model with interactive effects
\parencite{shi-lee-2017,bai-li-2021}:
\begin{equation}\label{eq:ysystem}
y_t=\rho W y_t+\phi y_{t-1}+X_t\beta + \mathbf{1}_{N} \delta_t+\Lambda f_t+\varepsilon_t,\qquad t=1,...,T,
\end{equation}
where
\begin{equation*}
y_t \in \mathbb R^N, \quad
f_t \in \mathbb R^r, \quad
\Lambda=(\lambda_1,...,\lambda_N)'\in\mathbb R^{N\times r}, \quad
\varepsilon_t \in \mathbb R^N.
\end{equation*}
The time-varying components of $f_t$ represent latent aggregate shocks or
common conditions that affect many or all cross-sectional units. Examples
include macroeconomic fluctuations, policy and regulatory changes,
technological innovations, financial conditions, commodity-price movements,
public-health shocks, and other economy-wide demand or supply disturbances.
The heterogeneous loadings $\lambda_i$ allow the magnitude and direction of
the response to these common shocks to differ across units.
Time-invariant individual heterogeneity can be incorporated by including a constant factor.
Define $B(\rho):=I_N-\rho W$.
Equation~\eqref{eq:ysystem} can be equivalently written as
\begin{equation}\label{eq:Bsystem}
B(\rho)y_t=\phi y_{t-1}+X_t\beta+\mathbf{1}_{N} \delta_t+\Lambda f_t+\varepsilon_t.
\end{equation}
We observe the initial cross section $y_0=(y_{10},...,y_{N0})'$ and condition on it throughout.
Let $x_{it}\in\mathbb R^K$ and stack $X_t$ as the $N\times K$ matrix with $i$th row $x_{it}'$.
Define the observed design or common conditioning $\sigma$-field as
\begin{equation*}
\mathcal C_N
:=
\sigma(y_0,X_1,...,X_T,W).
\end{equation*}
\subsection{Transformation and Jacobian}
Recall that $y_t\in\mathbb R^N$ denotes the cross section observed at time $t$.
In this subsection we also introduce unit-specific time paths such as $y_i\in\mathbb R^T$,
so the time index $t$ and the unit index $i$ refer to objects of different dimensions.
Define the transformed innovations
\begin{equation}\label{eq:utdef}
u_t(\rho,\phi,\beta):=B(\rho)y_t-\phi y_{t-1}-X_t\beta\in\mathbb R^N.
\end{equation}
Then \eqref{eq:Bsystem} implies
\begin{equation}\label{eq:uteq}
u_t=\mathbf{1}_{N} \delta_t + \Lambda f_t+\varepsilon_t,\qquad t=1,...,T.
\end{equation}
To make the dynamic structure explicit at the unit level, define for each $i$
\begin{equation*}
y_i := (y_{i1},...,y_{iT})' \in \mathbb R^T, \quad
X_i\beta := (x_{i1}'\beta,...,x_{iT}'\beta)' \in \mathbb R^T, \quad
\mathbf e_1 := (1,0,...,0)' \in \mathbb R^T,
\end{equation*}
and let $Y=(y_1,...,y_T)$ be the $N\times T$ matrix of observed cross sections.
For any matrix $C$, write $C_{i\cdot}$ for its $i$th row.
Introduce the $T\times T$ lower-triangular time-direction operator
\begin{equation}
R_T(\phi) :=
\begin{pmatrix}
1 & 0 & 0 & \cdots & 0 \\
-\phi & 1 & 0 & \cdots & 0 \\
0 & -\phi & 1 & \cdots & 0 \\
\vdots & \ddots & \ddots & \ddots & \vdots \\
0 & \cdots & 0 & -\phi & 1
\end{pmatrix}, \qquad
\det R_T(\phi) = 1.
\end{equation}
Conditional on the initial condition $y_{i0}$, the dynamic recursion can be written as
\begin{equation}\label{eq:utdef_unit}
u_i(\rho,\phi,\beta) = R_T(\phi) y_i - \rho (WY)_{i\cdot}' - X_i\beta - \phi y_{i0}\mathbf e_1.
\end{equation}
Equation~\eqref{eq:utdef_unit} isolates the time-direction transformation.
Since $R_T(\phi)$ is lower triangular with ones on the diagonal,
the dynamic filter contributes no parameter-dependent Jacobian term once we condition on $y_{i0}$.
The only nontrivial Jacobian term in the likelihood comes from the contemporaneous spatial filter $B(\rho)=I_N-\rho W$.
\paragraph{Jacobian.}
Consider the change of variables
\begin{equation*}
(y_1,...,y_T)\longmapsto(u_1,...,u_T)
\end{equation*}
given by \eqref{eq:utdef}, conditioning on $y_0$. Stack
\begin{equation*}
y:=(y_1',...,y_T')'\in\mathbb R^{NT},
\qquad
u:=(u_1',...,u_T')'\in\mathbb R^{NT}.
\end{equation*}
The Jacobian matrix $\mathcal J:=\partial u/\partial y'$ is block lower triangular:
\begin{equation}\label{eq:jacobian_matrix}
\mathcal J \;=\;
\begin{pmatrix}
B(\rho) & 0 & 0 & \cdots & 0 \\
-\phi I_N & B(\rho) & 0 & \cdots & 0 \\
0 & -\phi I_N & B(\rho) & \cdots & 0 \\
\vdots & \ddots & \ddots & \ddots & \vdots \\
0 & \cdots & 0 & -\phi I_N & B(\rho)
\end{pmatrix}.
\end{equation}
Hence
\begin{equation}\label{eq:jacobian}
|\mathcal J|=\prod_{t=1}^T |B(\rho)| = |B(\rho)|^T.
\end{equation}
This calculation makes clear that there are two distinct structural transformations.
The first is the time-direction operator $R_T(\phi)$ in \eqref{eq:utdef_unit},
whose determinant is one.
The second is the contemporaneous cross-sectional operator $B(\rho)$,
which generates the nontrivial Jacobian contribution.
Accordingly, the conditional likelihood contains the term $T\log|B(\rho)|$ and no additional Jacobian term involving $\phi$.
\subsection{Factors}
Let the factor path matrix be
\begin{equation*}
F:=\begin{pmatrix} f_1'\\ \vdots\\ f_T'\end{pmatrix}\in\mathbb R^{T\times r}.
\end{equation*}
For each unit $i$, stack innovations and idiosyncratic errors over time:
\begin{equation*}
u_i:=(u_{i1},...,u_{iT})'\in\mathbb R^T,\qquad
\varepsilon_i:=(\varepsilon_{i1},...,\varepsilon_{iT})'\in\mathbb R^T.
\end{equation*}
Passing from the time-stacked vector $(u_1',...,u_T')'$ to the unit-stacked collection $\{u_i\}_{i=1}^{N}$
is only a reordering of coordinates.
Equivalently, it is a fixed permutation of the $NT$ entries and therefore has absolute determinant of one.
Thus, no further Jacobian term arises when the likelihood is re-written as a product over units rather than over time.
Then \eqref{eq:uteq} is equivalent to
\begin{equation}\label{eq:ui_basic}
u_i = \delta + F\lambda_i + \varepsilon_i,\qquad i=1,...,N,
\end{equation}
where $\lambda_i\in\mathbb R^r$ is the $i$th row of $\Lambda$ as a column vector and $\delta=(\delta_1,...,\delta_T)'\in\mathbb R^T$.
\begin{remark}[Notation]
Objects indexed by $t$ are cross-sectional $N\times1$ vectors,
for example $y_t,u_t,\varepsilon_t\in\mathbb R^N$.
Objects indexed by $i$ are unit-specific $T\times1$ vectors, for example $y_i,u_i,\varepsilon_i\in\mathbb R^T$.
Thus, the same letter may denote different objects depending on whether it is indexed by $t$ or by $i$.
We keep this notation because the intended dimension is usually clear from the index and from the context.
\end{remark}
\subsection{Time-Heteroskedastic Idiosyncratic Errors}
We assume the idiosyncratic errors $\varepsilon_i$ are iid, allowing heteroskedasticity over time but no serial correlation.
The errors $\varepsilon_i$ are also independent of factor loadings $\Lambda$.
\begin{assumption}[Idiosyncratic Errors]\label{ass:varepsilon}
Conditional on $\mathcal C_N,$
$\{\varepsilon_i\}_{i=1}^N$
are identically distributed and independent across $i$, and, for some $\zeta>0$ and $C_\varepsilon <\infty$,
\[ \begin{aligned}
&\mathbb E(\varepsilon_i\mid\mathcal C_N)=0,
\qquad
\mathbb E(\varepsilon_i\varepsilon_i'\mid\mathcal C_N)
=
D_\varepsilon,
\qquad
D_\varepsilon=\mathrm{diag}(\sigma_1^2,\ldots,\sigma_T^2),\\
&(\varepsilon_1,\ldots,\varepsilon_N)
\perp\!\!\!\perp
\Lambda
\mid\mathcal C_N, \qquad
\sup_{N\ge1}
\mathbb E\!\left[
\|\varepsilon_i\|^{4+\zeta}
\,\middle|\,
\mathcal C_N
\right]
\le C_\varepsilon
\qquad\text{a.s.}
\end{aligned}
\]
\end{assumption}
\subsection{Spatially Enriched Random Loadings}
We model the factor loadings using a \emph{spatially enriched random-loadings} (SERL) specification.
The construction adapts the correlated-random-effects approach of
\textcite{mundlak-1978},
\textcite{chamberlain-1982},
and \textcite{wooldridge-2005}
to factor loadings in a dynamic spatial model.
In a static spatial Durbin panel,
\textcite{debarsy-2012} uses a Mundlak projection for a time-invariant individual effect with own and spatially lagged regressor averages.
SERL extends this construction to vector-valued factor loadings in a dynamic model and allows the controls to include the initial outcome and higher-order spatial transformations.
It allows the conditional loading means to depend on observed covariates,
the initial condition,
and information propagated through the spatial structure.
Let $Q_X\in\mathbb R^{N\times p}$ collect observed summaries of the regressors.
A Mundlak-type specification may use time averages of the regressors,
whereas a Chamberlain-type specification may use the complete observed regressor history.
For a fixed integer $L\geq0$, define
\begin{equation}\label{eq:correct_projection}
z_i
=
\Big[
(Q_X)_{i\cdot}', (WQ_X)_{i\cdot}', \cdots, (W^LQ_X)_{i\cdot}',
(y_0)_i, (Wy_0)_i, \cdots, (W^Ly_0)_i
\Big]'
\in\mathbb R^{(L+1)(p+1)}.
\end{equation}
Because $Q_X$ is constructed from the observed regressor history and
$z_i$ is constructed from $Q_X$, $y_0$, and $W$,
both $Q_X$ and $z_i$ are $\mathcal C_N$-measurable.
We specify
\begin{equation}\label{eq:proj}
\lambda_i = A z_i + \eta_i,
\end{equation}
where $Az_i$ captures the component of the loading heterogeneity
systematically related to the observed information and $\eta_i$ represents
the idiosyncratic loading heterogeneity.
We refer to this as the SERL specification. It has a natural interpretation in terms of heterogeneous responses to common shocks. The loading $\lambda_i$ determines how unit $i$ responds to the common shocks $f_t$, and this response heterogeneity may be systematically related to observed characteristics, initial conditions, and their spatial transformations. The decomposition
$
\lambda_i=Az_i+\eta_i
$
separates such systematic heterogeneity from unit-specific heterogeneity. The first component, $Az_i$, is predictable from the common conditioning information $\mathcal{C}_N$ and captures systematic or coordinated responses associated with observed characteristics, initial conditions, and their spatial transformations.
The second component, $\eta_i$, captures purely unit-specific heterogeneity in the response to the common shock that is not systematically related to the information in $\mathcal{C}_N$, and is assumed to be independent across $i$ conditional on $\mathcal{C}_N$.
In this sense,
SERL adapts the correlated-random-effects idea to factor loadings in a spatial environment. Assumption~\ref{ass:SERL}
formalizes this decomposition.
\begin{assumption}[Spatially enriched random loadings]\label{ass:SERL}
For some fixed integer $L\geq0$ and coefficient matrix
$A\in\mathbb R^{r\times(L+1)(p+1)}$,
the loading vectors satisfy \eqref{eq:proj}.
Conditional on $\mathcal C_N$,
the vectors $\{\eta_i\}_{i=1}^N$ are identically distributed and independent across $i$, with
\[
\mathbb E(\eta_i\mid\mathcal C_N)=0,
\qquad
\mathbb E(\eta_i\eta_i'\mid\mathcal C_N)=\Sigma_\eta, \qquad \sup_{N\ge1}
\mathbb E\!\left[
\|\eta_i\|^{4+\zeta}
\,\middle|\,
\mathcal C_N
\right]
\le C_\eta< \infty
\qquad\text{a.s.}
\]
where $\Sigma_\eta$ is positive definite.
\end{assumption}
We refer to $L$ as the \emph{spatial enrichment order}.
The case $L=0$ uses only own-unit regressor summaries and the own initial outcome.
For $L>0$,
neighboring information enters through $WQ_X,\ldots,W^LQ_X$ and $Wy_0,\ldots,W^Ly_0$.
For the theoretical analysis, $L$ is treated as fixed. When the conditional loading mean arises from a richer spatial process,
a finite value of $L$ may instead be viewed as an approximation, whose adequacy can be assessed by sensitivity
to additional spatial transformations. Accordingly, in the empirical application, we examine alternative
enrichment orders to assess the sensitivity of the structural estimates to the choice of $L$.
If the decomposition is written as $\lambda_i=a+Az_i+\eta_i$ with $a\in\mathbb R^r$ as an intercept,
the common term $Fa$ is absorbed by the unrestricted time effect $\delta$.
Therefore, we omit $a$ without loss of generality.
Appendix~\ref{sec:joint-initial-condition} develops an alternative approach
that models $y_0$ jointly with the sample outcomes and projects the loadings only on regressor-based controls.
\subsection{Likelihood Conditional on \texorpdfstring{$\mathcal C_N$}{C_N}}
Substituting \eqref{eq:proj} into \eqref{eq:ui_basic} yields
\begin{equation*}
u_i = \delta + F(Az_i+\eta_i)+\varepsilon_i = \delta + FAz_i + F\eta_i + \varepsilon_i.
\end{equation*}
Define
\begin{equation} \label{eq:ei}
e_i := F \eta_i +\varepsilon_i
\end{equation}
By Assumptions~\ref{ass:varepsilon} and \ref{ass:SERL},
conditional on $\mathcal C_N$,
the vectors $\{e_i\}_{i=1}^N$ are identically distributed and independent across $i$, with
\begin{equation} \label{eq:conditional_moments}
\mathbb E(e_i\mid\mathcal C_N)=0, \qquad
\mathbb E(e_ie_i'\mid\mathcal C_N)
=
F\Sigma_\eta F'+D_\varepsilon,
\qquad
\sup_{N\ge1}
\mathbb E\!\left[
\|e_i\|^{4+\zeta}
\,\middle|\,
\mathcal C_N
\right]
\le C_e, \quad \text{a.s.}
\end{equation}
where $C_e<\infty$.
Hence,
\begin{equation}\label{eq:cond_mean}
\mathbb E(u_i\mid\mathcal C_N) = \delta + FAz_i.
\end{equation}
and
\begin{equation}\label{eq:cond_var}
\mathrm{Var}(u_i\mid\mathcal C_N) = F\Sigma_\eta F' + D_\varepsilon \;=:\; \Sigma_u.
\end{equation}
For the moment, let $\alpha$ denote the collection of unknown parameters entering the working likelihood. A precise parameterization of $\alpha$ in terms of free parameters is introduced in Section~\ref{sec:asymptotics} after imposing a normalization on the factor path $F$ to remove its rotational indeterminacy.
No Gaussian distribution is imposed on $\eta_i$, $\varepsilon_i$, or $e_i$. The Gaussian specification used below is a working likelihood based on the conditional mean and covariance in \eqref{eq:cond_mean}--\eqref{eq:cond_var}.
Combining the determinant-one time transformation in \eqref{eq:utdef_unit},
the spatial Jacobian in \eqref{eq:jacobian}, and
the determinant-one permutation from time stacking to unit stacking,
one obtains a conditional likelihood with a single Jacobian term, namely
$|B(\rho)|^T$. Let $Y=(y_1,\ldots,y_T)$. Motivated by the conditional moment restrictions,
we use the conditional Gaussian working likelihood
\begin{equation}\label{eq:lik}
\begin{aligned}
L_N(\alpha;Y\mid\mathcal C_N)
&=
|B(\rho)|^T
\prod_{i=1}^N
(2\pi)^{-T/2}
|\Sigma_u|^{-1/2} \times \\
& \exp\!\left(
-\tfrac12 (u_i-\delta-FAz_i)'\Sigma_u^{-1}(u_i-\delta-FAz_i)
\right),
\end{aligned}
\end{equation}
where $u_i=(u_{i1},\ldots,u_{iT})'$ is obtained by rearranging $u_t$ computed from observed data in \eqref{eq:utdef}.
The corresponding working log-likelihood, up to an additive constant, is given by
\begin{equation}\label{eq:loglik}
\ell(\alpha)
=
T\log|B(\rho)|
-\frac{N}{2}\log|\Sigma_u|
-\frac12 \sum_{i=1}^N (u_i-\delta-FAz_i)'\Sigma_u^{-1}(u_i-\delta-FAz_i).
\end{equation}
\section{Inferential Theory under Large \texorpdfstring{$N$}{N} and Fixed \texorpdfstring{$T$}{T}}\label{sec:asymptotics}
This section states the large-$N$, fixed-$T$ theory for the Gaussian QML
estimator. Gaussianity is not required. The key complication is that the
spatial transformation makes the observed-data score cross-sectionally
dependent, even when the structural residuals are independent across units.
Consequently, inference must use spatially corrected score contributions rather
than raw individual score outer products.
All expectations, variances, and probability statements in this section are
understood conditionally on $\mathcal C_N$ unless otherwise indicated.
\subsection{Parameterization, Normalization, and Sample Criterion}\label{subsec:qml_setup}
The key parameters are
\begin{equation*}
\theta := (\rho,\ \phi,\ \beta')'
\end{equation*}
and the factor path is identified only up to rotation. We assume that its top
$r\times r$ block has full rank and use the normalized representation
\begin{equation}\label{eq:F_normalization_theory}
F=
\begin{pmatrix}
I_r\\ F_2
\end{pmatrix},
\qquad F_2\in\mathbb R^{(T-r)\times r}.
\end{equation}
Write $D_\varepsilon=\mathrm{diag}(\sigma_1^2,...,\sigma_T^2)$ and
$\sigma^2=(\sigma_1^2,...,\sigma_T^2)'$. The nuisance vector contains only
free coordinates,
\begin{equation*}
\varphi
=
\left(
\delta',\operatorname{vec}(A)',\operatorname{vec}(F_2)',
\mathrm{vech}(\Sigma_\eta)',\sigma^{2\prime}
\right)',
\qquad
\alpha=(\theta',\varphi')'.
\end{equation*}
Thus, only $F_2$ is included in $\alpha$, and $\Sigma_\eta$ is represented by its unique elements. The parameter space restricts $\Sigma_\eta$ to be positive definite
and every element of $\sigma^2$ to be positive. For each unit $i$, define the transformed innovations
\begin{equation*}
u_i(\theta):=(u_{i1}(\theta),...,u_{iT}(\theta))'\in\mathbb R^T
\quad\text{where}\quad
u_t(\theta)=B(\rho)y_t-\phi y_{t-1}-X_t\beta
\end{equation*}
for $t=1,...,T$.
The nuisance mean is
\begin{equation*}
\mu_i(\varphi):=\delta+FAz_i\in\mathbb R^T,
\end{equation*}
and the inverse of the working covariance matrix is
\begin{equation*}
M(\varphi):=\Sigma_u(\varphi)^{-1}
\quad\text{where}\quad
\Sigma_u(\varphi):=D_\varepsilon+F\Sigma_\eta F'.
\end{equation*}
Define residuals $e_i(\alpha):=u_i(\theta)-\mu_i(\varphi)$.
Note
$e_i(\alpha_0)=F\eta_i+\varepsilon_i$, which is $e_i$ defined in \eqref{eq:ei}, and satisfies
\eqref{eq:conditional_moments}.
The per-unit Gaussian pseudo log-likelihood is
\begin{equation*}
\ell_i(\alpha)
:=
-\frac12\log|\Sigma_u(\varphi)|
-\frac12\,e_i(\alpha)'M(\varphi)e_i(\alpha),
\end{equation*}
and the sample criterion, normalized by $N$, is
\begin{equation}\label{eq:ell_N}
\ell_N(\alpha)
:=
\frac{T}{N}\log|B(\rho)|+\frac1N\sum_{i=1}^N \ell_i(\alpha).
\end{equation}
Let $\widehat\alpha=(\widehat\theta',\widehat\varphi')'$ be an interior maximizer of $\ell_N(\alpha)$,
and hence a solution to the corresponding first-order conditions.
\subsection{Consistency and Asymptotic Normality}\label{subsec:asymptotic_normality}
Define the ordinary unit score by
\begin{equation*}
s_i(\alpha)
:=
\nabla_\alpha\ell_i(\alpha)
+
\frac{T}{N}\nabla_\alpha\log|B(\rho)|.
\end{equation*}
The normalized pooled score therefore satisfies
\begin{equation*}
S_N(\alpha)
:= \frac1N \sum_{i=1}^N s_i(\alpha)
= \nabla_\alpha\ell_N(\alpha).
\end{equation*}
The spatial filter makes the ordinary unit scores cross-sectionally dependent even when the structural residuals are independent.
In the baseline model,
the dependence requiring correction is due to the $\rho$ and $\phi$ components,
since their reduced-form regressors depend on structural residuals from other units.
Consequently,
\begin{equation}\label{eq:Omega_raw_scores}
\mathrm{Var}_N\left(\frac1{\sqrt N}\sum_{i=1}^N s_i(\alpha_0)\right)
=
\frac1N\sum_{i=1}^N\sum_{j=1}^N
\mathrm{Cov}_N\!\big(s_i(\alpha_0),s_j(\alpha_0)\big),
\end{equation}
where, for notational simplicity, we write
\begin{equation*}
\mathbb E_N(\,\cdot\,):=\mathbb E(\,\cdot\mid\mathcal C_N),
\qquad
\mathrm{Var}_N(\,\cdot\,):=\mathrm{Var}(\,\cdot\mid\mathcal C_N),
\qquad
\mathrm{Cov}_N(\,\cdot\,,\,\cdot\,):=\mathrm{Cov}(\,\cdot\,,\,\cdot\mid\mathcal C_N).
\end{equation*}
The diagonal-only outer product based on the raw scores generally omits the cross-unit terms in \eqref{eq:Omega_raw_scores}.
The cross-sectional dependence only enters the $\rho$ and $\phi$ components of the score,
so the remaining components already have the required unitwise form.
Following \textcite{li-yang-2021},
we therefore rearrange the cross-unit terms in the $\rho$ and $\phi$ scores while preserving the pooled score.
Let $\{g_i(\alpha)\}_{i=1}^{N}$ denote the resulting rearranged contributions,
satisfying
\begin{equation}\label{eq:g_score_decomposition}
N S_N(\alpha)
=
\sum_{i=1}^N s_i(\alpha)
=
\sum_{i=1}^N g_i(\alpha),
\end{equation}
and that at $\alpha_0$,
the rearranged contributions form a martingale-difference array with respect to the filtration
\begin{equation*}
\mathcal F_i
:=
\mathcal C_N\vee\sigma(e_1,...,e_i).
\end{equation*}
Appendix~\ref{sec:spatial_score_construction} gives the exact component-by-component construction of $g_i$.
\begin{assumption}[Regularity]\label{ass:regularity}
\begin{enumerate}
\item
The parameter space is compact and $\alpha_0$ is an interior point.
The criterion converges uniformly in probability over the parameter space to a limiting criterion uniquely maximized at $\alpha_0$.
\item
The criterion is twice continuously differentiable on a convex neighborhood $\mathcal A_0$ of $\alpha_0$.
\item
The row and column sums of the spatial weights matrices are uniformly bounded,
and $B(\rho)^{-1}$ and the dynamic-spatial reduced-form operator are uniformly bounded over the parameter space.
\item
On $\mathcal A_0$,
the eigenvalues of $\Sigma_u(\varphi)$ are uniformly bounded away from zero and infinity,
and its first two derivatives are uniformly bounded.
The unit-to-unit reduced-form response blocks entering $g_i(\alpha)$,
together with their first two parameter derivatives,
admit a common summable envelope whose block row and column sums are uniformly bounded.
The conditionally nonrandom quantities entering the score contributions,
including the observed regressors, projection controls, initial-outcome terms, and reduced-form means,
together with their first two derivatives,
have uniformly bounded cross-sectional average moments of order $4+\zeta$.
\item
There is a matrix function $H(\alpha)$,
continuous at $\alpha_0$,
such that
\begin{equation*}
\sup_{\alpha\in\mathcal A_0}
\left\|
-\nabla_{\alpha\alpha'}^2\ell_N(\alpha)-H(\alpha)
\right\|
\overset{p}\longrightarrow0,
\end{equation*}
where $H(\alpha_0)$ is finite and positive definite.
\end{enumerate}
\end{assumption}
These are standard identification and smoothness requirements for
M-estimation, augmented by stability conditions for the spatial reduced form.
With fixed $T$,
the dimension of $\alpha$ is fixed,
so consistency and the asymptotic-linearization argument follow standard finite-dimensional QMLE reasoning,
as in \textcite{white-1982}.
What is nonstandard here is the cross-sectional dependence of the raw spatial scores;
the rearrangement into $g_i$ and the martingale central limit theorem below provide the required score limit.
\begin{assumption}[Predictable score covariance]\label{ass:score_clt}
For the rearranged contributions $g_i(\alpha_0)$ and the filtration
$\mathcal F_i=\mathcal C_N\vee\sigma(e_1,...,e_i)$,
the predictable quadratic variation satisfies
\[
\frac1N\sum_{i=1}^N
\mathbb E\!\left[
g_i(\alpha_0)g_i(\alpha_0)'
\mid\mathcal F_{i-1}
\right]
\overset{p}\longrightarrow\Omega(\alpha_0),
\]
where $\Omega(\alpha_0)$ is finite and positive definite.
\end{assumption}
\begin{remark}[Lindeberg Condition]
The moment bounds established in Lemma~\ref{lem:score-moment-bounds} of Appendix~\ref{sec:spatial_score_construction},
based on Assumptions~\ref{ass:varepsilon}, \ref{ass:SERL}, and \ref{ass:regularity},
show that,
for some $\delta>0$,
\begin{equation*}
\frac1N\sum_{i=1}^N
\mathbb E_N\|g_i(\alpha_0)\|^{2+\delta}
=O_p(1).
\end{equation*}
Since $\mathcal C_N\subseteq\mathcal F_{i-1}$,
Markov's inequality implies that
\begin{equation*}
\frac1N\sum_{i=1}^N
\mathbb E\!\left[
\|g_i(\alpha_0)\|^{2+\delta}
\,\middle|\,\mathcal F_{i-1}
\right]
=O_p(1).
\end{equation*}
Consequently,
for every $\epsilon>0$,
the conditional Lindeberg condition follows:
\begin{equation*}
\begin{aligned}
\frac1N\sum_{i=1}^N
\mathbb E\!\left[
\|g_i(\alpha_0)\|^2
\mathbf 1\{\|g_i(\alpha_0)\|>\epsilon\sqrt N\}
\middle|\mathcal F_{i-1}
\right]
\leq
\frac{1}{\epsilon^\delta N^{\delta/2}}
\frac1N\sum_{i=1}^N
\mathbb E\!\left[
\|g_i(\alpha_0)\|^{2+\delta}
\,\middle|\,\mathcal F_{i-1}
\right]
=o_p(1).
\end{aligned}
\end{equation*}
\end{remark}
Because the rearranged contributions form a martingale-difference array,
the Lindeberg condition,
Assumption~\ref{ass:score_clt},
and Theorem~3.2 and Corollary~3.1 of \textcite[pp.~58--59]{hall-heyde-1980} yield
\begin{equation*}
\sqrt N S_N(\alpha_0)
=
\frac1{\sqrt N}\sum_{i=1}^N g_i(\alpha_0)
\Rightarrow
\mathcal N\!\left(0,\Omega(\alpha_0)\right).
\end{equation*}
The uniform-convergence and unique-maximizer conditions in Assumption~\ref{ass:regularity},
together with the standard argmax theorem,
give $\widehat\alpha\overset{p}\longrightarrow\alpha_0$.
A mean-value expansion of the first-order condition around $\alpha_0$ then yields
\begin{equation*}
\sqrt N(\widehat\alpha-\alpha_0)
=
H(\alpha_0)^{-1}\sqrt N S_N(\alpha_0)+o_p(1).
\end{equation*}
This gives the following result.
\begin{theorem}[Consistency and asymptotic normality]\label{thm:qml_asymptotics}
Under Assumptions~\ref{ass:varepsilon}--\ref{ass:score_clt},
\begin{equation}\label{eq:alpha_CLT}
\widehat\alpha\overset{p}\longrightarrow\alpha_0,
\qquad
\sqrt N(\widehat\alpha-\alpha_0)
\Rightarrow
\mathcal N\Big(0,H(\alpha_0)^{-1}\Omega(\alpha_0)(H(\alpha_0)^{-1})'\Big).
\end{equation}
\end{theorem}
Because the criterion is a Gaussian
pseudo-likelihood, the information equality need not hold. Estimating
$\Omega(\alpha_0)$ therefore requires separate analysis.
\subsection{Estimation of the Asymptotic Covariance}\label{subsec:sandwich_general}
The martingale representation in \eqref{eq:g_score_decomposition} identifies the score covariance with the probability limit of its predictable quadratic variation:
\begin{equation*}
\Omega(\alpha_0)
=
\operatorname*{plim}_{N\to\infty}
\frac1N\sum_{i=1}^N
\mathbb E\!\left[
g_i(\alpha_0)g_i(\alpha_0)'
\mid\mathcal F_{i-1}
\right].
\end{equation*}
While the full outer product of the rearranged contributions is consistent for this covariance,
we further refine the estimator to remove terms that have zero conditional expectation under the maintained martingale structure.
The full construction of the estimator $\widehat\Omega_M(\widehat\alpha)$ is given in Appendix~\ref{sec:spatial_score_construction}.
Proposition~\ref{prop:feasible_score_covariance} proves the asymptotic equivalence
\begin{equation*}
\widehat\Omega_M(\alpha_0)
-
\frac1N\sum_{i=1}^N g_i(\alpha_0)g_i(\alpha_0)'
=o_p(1).
\end{equation*}
It also proves the plug-in replacement
$\widehat\Omega_M(\widehat\alpha)-\widehat\Omega_M(\alpha_0)=o_p(1)$,
and hence $\widehat\Omega_M(\widehat\alpha)\overset{p}\longrightarrow\Omega(\alpha_0)$.
For improved finite-sample performance,
our preferred implementation applies a pair-only HC1-style correction.
Using the decomposition in \eqref{eq:martingale_covariance_estimator},
write
$\widehat\Omega_M=\widehat\Omega_\dagger+\widehat\Omega_v$,
where $\widehat\Omega_\dagger$ is the one-unit component and
$\widehat\Omega_v$ is the pair component.
We define
\begin{equation}\label{eq:Omega_martingale_pair_hc1}
\widehat\Omega_{M,\mathrm{pair}}(\widehat\alpha)
:=
\widehat\Omega_\dagger(\widehat\alpha)
+
\frac{N}{N-p_\alpha}
\widehat\Omega_v(\widehat\alpha),
\end{equation}
where $p_\alpha:=\dim(\alpha)$ counts all estimated structural and nuisance parameters.
This adjustment targets the component associated with the spatial and dynamic score rearrangement while leaving the one-unit component unchanged.
We use it as a simple finite-sample correction and do not claim that it makes the pair component exactly unbiased.
The effective sample size is $N$ because each cross-sectional unit contributes one fixed-length time path to the score decomposition.
Since $T$ and $p_\alpha$ are fixed as $N\to\infty$,
the multiplier converges to one and does not change the covariance estimator's probability limit.
We use $\widehat\Omega_{M,\mathrm{pair}}$ for our reported simulation and application inference below,
except where a table explicitly labels otherwise.
Let $\widehat H:=-\nabla_{\alpha\alpha'}^2\ell_N(\widehat\alpha)$.
The preferred feasible sandwich covariance estimator for the full estimator is
\begin{equation}\label{eq:Valpha_hat}
\widehat{\mathrm{Var}}(\widehat\alpha)
=
\frac1N\widehat H^{-1}\widehat\Omega_{M,\mathrm{pair}}(\widehat\alpha)(\widehat H^{-1})'.
\end{equation}
The uncorrected benchmark replaces $\widehat\Omega_{M,\mathrm{pair}}$ by $\widehat\Omega_M$.
\subsection{Profile Inference for the Structural Parameters}\label{subsec:profile_sandwich}
Partition the negative Hessian conformably with
$\alpha=(\theta',\varphi')'$ and define
\begin{equation}\label{eq:H_schur}
H_{\theta\mathbin{\cdot}\varphi}
=
H_{\theta\theta}
-H_{\theta\varphi}H_{\varphi\varphi}^{-1}H_{\varphi\theta},
\qquad
P=H_{\theta\varphi}H_{\varphi\varphi}^{-1}.
\end{equation}
The profiled corrected contribution is
\begin{equation}\label{eq:g_profile}
g_{\theta\mathbin{\cdot}\varphi,i}(\alpha_0)
=
g_{\theta,i}(\alpha_0)-P g_{\varphi,i}(\alpha_0),
\end{equation}
with covariance
\begin{equation*}
\Omega_{\theta\mathbin{\cdot}\varphi}
=
\operatorname*{plim}\frac1N\sum_{i=1}^N
\mathbb E[g_{\theta\mathbin{\cdot}\varphi,i}g_{\theta\mathbin{\cdot}\varphi,i}'].
\end{equation*}
It follows that
\begin{equation}\label{eq:theta_CLT}
\sqrt N(\widehat\theta-\theta_0)
\Rightarrow
\mathcal N\!\left(
0,
H_{\theta\mathbin{\cdot}\varphi}^{-1}
\Omega_{\theta\mathbin{\cdot}\varphi}
(H_{\theta\mathbin{\cdot}\varphi}^{-1})'
\right).
\end{equation}
For feasible inference,
let $\widehat P$ be the sample analogue of $P$ in \eqref{eq:H_schur}.
The preferred profiled covariance is the corresponding transformation of the pair-corrected martingale covariance estimator:
\begin{equation*}
\widehat\Omega_{\theta\mathbin{\cdot}\varphi}
=
\begin{bmatrix}
I & -\widehat P
\end{bmatrix}
\widehat\Omega_{M,\mathrm{pair}}(\widehat\alpha)
\begin{bmatrix}
I \\ -\widehat P'
\end{bmatrix}.
\end{equation*}
Then
\begin{equation}\label{eq:Vtheta_hat}
\widehat{\mathrm{Var}}(\widehat\theta)
=
\frac1N
\widehat H_{\theta\mathbin{\cdot}\varphi}^{-1}
\widehat\Omega_{\theta\mathbin{\cdot}\varphi}
(\widehat H_{\theta\mathbin{\cdot}\varphi}^{-1})'.
\end{equation}
\section{Estimation Algorithm}\label{sec:estimation_algorithm}
This section outlines a practical block algorithm to maximize the conditional Gaussian (QML)
log-likelihood \eqref{eq:loglik} when $T$ is fixed and the number of factors $r$ is small.
The parameters of interest are $(\rho,\phi,\beta)$,
while $(\delta,F,A,\Sigma_\eta,D_\varepsilon)$ are nuisance parameters.
We exploit three features:
(i) for given $(F,\Sigma_u)$, $(\delta,A)$ can be \emph{concentrated out in closed form};
(ii) under the current values of $(\delta,F,A)$, the Gaussian working model for the random loading heterogeneity $\eta_i$
yields simple EM/ECM updates for $(F,\Sigma_\eta,D_\varepsilon)$;
(iii) given $\rho$, the parameters $(\phi,\beta)$ admit a closed-form generalized least squares update and a one-dimensional conditional
maximization step in $\rho$ can be implemented using the Jacobian term $\log|B(\rho)|$.
This yields an inner-outer loop structure:
the outer loop performs a scalar conditional maximization in $\rho$, and the inner loop is a block coordinate/ECME iteration
for $(\phi,\beta)$ and $(\delta,A,F,\Sigma_\eta,D_\varepsilon)$ given $\rho$.
\subsection{Block Coordinate/ECME Iteration: Inner Loop}
Let $u_t(\rho,\phi,\beta)$ be defined in \eqref{eq:utdef} and stack $U=(u_1,...,u_N)\in\mathbb R^{T\times N}$,
$Z=(z_1,...,z_N)\in\mathbb R^{q\times N}$. Given current parameter values, iterate the following steps
until convergence, e.g., change in concentrated log-likelihood below a tolerance.
\paragraph{Inputs to Steps 1--4.}
Given current $(\rho,\phi,\beta)$, form the innovation (or transformed residual) matrix
\begin{equation*}
U(\rho,\phi,\beta):=(u_1,...,u_N)\in\mathbb R^{T\times N},\qquad
u_t(\rho,\phi,\beta):=B(\rho)y_t-\phi y_{t-1}-X_t\beta,\ \ t=1,...,T,
\end{equation*}
and let $u_i(\rho,\phi,\beta)\in\mathbb R^T$ denote the $i$th column of $U(\rho,\phi,\beta)$.
In Steps~1--4 below we treat $U(\rho,\phi,\beta)$ as the \emph{data input} and update only nuisance parameters
$(\delta,A,F,\Sigma_\eta,D_\varepsilon)$.
\paragraph{Step 1 (Input: the transformed residual matrix; concentrate out time effects $\delta$ and projection matrix $A$).}
Given current $(F,\Sigma_\eta,D_\varepsilon)$, compute
\begin{equation*}
\Sigma_u=D_\varepsilon+F\Sigma_\eta F',\qquad M:=\Sigma_u^{-1}.
\end{equation*}
This step can be computed efficiently using Woodbury identity; see Section~\ref{sec:woodbury}. Let
\begin{equation*}
\bar u:=\frac1N U(\rho,\phi,\beta)\,\mathbf{1}_{N},\quad \bar z:=\frac1N Z\,\mathbf{1}_{N},\quad
U_c:=U(\rho,\phi,\beta)-\bar u\,\mathbf{1}_{N}',\quad Z_c:=Z-\bar z\,\mathbf{1}_{N}'.
\end{equation*}
Then update $A$ by the profiled GLS formula
\begin{equation*}
A \leftarrow \widehat A(F;U)
=
\big(F'MF\big)^{-1}\,F'MU_c\,Z_c'\,\big(Z_cZ_c'\big)^{-1},
\end{equation*}
and update the free time effect by
\begin{equation*}
\delta \leftarrow \widehat\delta(F,A;U)=\bar u - F A\,\bar z.
\end{equation*}
With this ordering, $\widehat A(F;U)$ depends on $(u_i-\bar u)$ and $(z_i-\bar z)$.
See Section~\ref{sec:profile_delta_A} for details on the concentration of $\delta$ and $A$.
The displayed updates are analogous to first estimating a regression slope and then its intercept.
Equivalently, one may first concentrate out $\delta$, solve for $A$, and substitute the resulting $A$ back into the expression for $\delta$;
concentrating out $(\delta,A)$ jointly yields the same formulas.
\paragraph{Step 2 (E-step for the random loading heterogeneity $\eta_i$ given the transformed residual matrix).}
For each $i$, define
\begin{equation*}
e_i := u_i(\rho,\phi,\beta)-\delta-FAz_i \in\mathbb R^T.
\end{equation*}
Treating $\eta_i$ as missing data, the posterior is Gaussian with
\begin{align*}
V &:= \mathrm{Var}(\eta_i\mid u_i,z_i)
=
\big(\Sigma_\eta^{-1}+F'D_\varepsilon^{-1}F\big)^{-1},\\
m_i &:= \mathbb E(\eta_i\mid u_i,z_i)
=
V\,F'D_\varepsilon^{-1}e_i,
\qquad i=1,...,N.
\end{align*}
(Here $V$ is common across $i$ and is only $r\times r$.)
\paragraph{Step 3 (M-step for factor covariance $\Sigma_\eta$ and idiosyncratic variance matrix $D_\varepsilon$).}
Update the loading covariance by
\begin{equation*}
\Sigma_\eta \leftarrow V + \frac1N\sum_{i=1}^N m_i m_i'.
\end{equation*}
Let $r_i:=e_i-Fm_i$ denote the posterior mean of $\varepsilon_i$.
Since $\mathbb E(\varepsilon_i\varepsilon_i'\mid u_i)=r_i r_i' + FVF'$, the diagonal variances update elementwise as
\begin{equation*}
\sigma_t^2 \leftarrow \frac1N\sum_{i=1}^N r_{it}^2 \;+\; f_t'Vf_t,\qquad t=1,...,T,
\end{equation*}
where $f_t'\in\mathbb R^{1\times r}$ is row $t$ of $F$, and set $D_\varepsilon\leftarrow\mathrm{diag}(\sigma_1^2,...,\sigma_T^2)$.
\paragraph{Step 4 (M-step for factor path $F$; impose normalization for factor path after updating).}
Define
\begin{equation*}
s_i:=Az_i+m_i\in\mathbb R^r,\qquad S:=(s_1,...,s_N)\in\mathbb R^{r\times N}.
\end{equation*}
The closed-form ECM update for $F$ is
\begin{equation*}
\widetilde F \leftarrow \big(U(\rho,\phi,\beta)-\delta\,\mathbf{1}_{N}'\big)\,S'\,\big(SS' + NV\big)^{-1}.
\end{equation*}
Since $F$ is only identified up to rotation, we can post-process $\widetilde F$ by any invertible $r\times r$ matrix
$H$ without changing the likelihood.
Impose the normalization $F_{1:r,:}=I_r$ by the rotation
\begin{equation*}
H := \big(\widetilde F_{1:r,:}\big)^{-1},\qquad
F \leftarrow \widetilde F H,\qquad
A \leftarrow H^{-1}A,\qquad
\Sigma_\eta \leftarrow H^{-1}\Sigma_\eta (H^{-1})'.
\end{equation*}
\begin{remark}[Normalization with Time-Invariant Individual Heterogeneity]
If time-invariant individual heterogeneity is included,
the first factor is normalized to be constant over time,
$f_{t1}=1$,
so that
\[
F=[\mathbf{1}_{T},G],
\qquad
G\in\mathbb R^{T\times(r-1)},
\]
where $r$ is the total number of factors,
including the constant factor.
In this case,
the identity-block normalization in \eqref{eq:F_normalization_theory} is replaced by
\[
F_{1:r,:}
=
\begin{pmatrix}
1 & 0_{1\times(r-1)}\\
\mathbf{1}_{r-1} & I_{r-1}
\end{pmatrix},
\]
or equivalently,
\[
G_{1,\cdot}=0,
\qquad
G_{2:r,\cdot}=I_{r-1}.
\]
This normalization preserves the constant first factor while removing the remaining location and rotation indeterminacy.
The factor update and subsequent rotation must then be restricted to preserve the constant column.
The corresponding loading follows the same SERL representation as the other components of $\lambda_i$;
its systematic component is captured by $Az_i$ and its residual variation by $\Sigma_\eta$,
without estimating a separate parameter for every unit.
When no constant factor is imposed,
we retain the normalization $F_{1:r,:}=I_r$ used above.
\end{remark}
\paragraph{Step 5 (joint GLS update for the dynamic and slope coefficients, conditional on $\rho$).}
Fix $\rho$ and the current nuisance parameters $(\delta,F,A,\Sigma_\eta,D_\varepsilon)$,
hence $\Sigma_u=D_\varepsilon+F\Sigma_\eta F'$ and $M:=\Sigma_u^{-1}$.
Define the spatially transformed dependent variable
\begin{equation*}
\widetilde y_t(\rho):=B(\rho)y_t = (I_N-\rho W)y_t,\qquad t=1,...,T,
\end{equation*}
and for each unit $i$ stack
\begin{equation*}
\widetilde y_i(\rho):=\big(\widetilde y_{i1}(\rho),...,\widetilde y_{iT}(\rho)\big)'\in\mathbb R^T,\qquad
y_{i,-1}:=(y_{i0},y_{i1},...,y_{i,T-1})'\in\mathbb R^T,
\end{equation*}
\begin{equation*}
X_i:=\begin{pmatrix}x_{i1}'\\ \vdots\\ x_{iT}'\end{pmatrix}\in\mathbb R^{T\times K},\qquad
R_i:=[\,y_{i,-1}\ \ X_i\,]\in\mathbb R^{T\times (1+K)},\qquad
\vartheta:=\begin{pmatrix}\phi\\ \beta\end{pmatrix}\in\mathbb R^{1+K}.
\end{equation*}
Then the mean equation implied by \eqref{eq:Bsystem} and the SERL specification is
\begin{equation*}
\widetilde y_i(\rho)=R_i\vartheta+\delta+FAz_i+v_i,\qquad \mathrm{Var}(v_i\mid\mathcal C_N)=\Sigma_u.
\end{equation*}
\emph{Centered (profile-$\delta$) GLS.}
Since $\delta$ is a free $T\times 1$ time effect common across $i$, it is convenient to
difference out $\delta$ by cross-sectional demeaning. Let
\begin{equation*}
\bar{\widetilde y}(\rho):=\frac1N\sum_{i=1}^N \widetilde y_i(\rho),\quad
\bar R:=\frac1N\sum_{i=1}^N R_i,\quad
\bar z:=\frac1N\sum_{i=1}^N z_i,
\end{equation*}
and define centered objects
\begin{equation*}
\widetilde y_{i,c}(\rho):=\widetilde y_i(\rho)-\bar{\widetilde y}(\rho),\qquad
R_{i,c}:=R_i-\bar R,\qquad
z_{i,c}:=z_i-\bar z.
\end{equation*}
Then the (conditional) GLS update for $\vartheta=(\phi,\beta')'$ is the explicit closed form
\begin{equation}\label{eq:theta_GLS}
\vartheta \leftarrow \widehat\vartheta(\rho)
=
\left(\sum_{i=1}^N R_{i,c}'MR_{i,c}\right)^{-1}
\left(\sum_{i=1}^N R_{i,c}'M\big(\widetilde y_{i,c}(\rho)-FAz_{i,c}\big)\right).
\end{equation}
Equivalently, one may use the uncentered version with $\delta$ explicitly present; \eqref{eq:theta_GLS}
is exactly the profile-$\delta$ GLS solution because $M$ is common across $i$.
Given $\widehat\vartheta(\rho)$, form the transformed innovations
\begin{equation*}
u_t(\rho,\widehat\vartheta):=\widetilde y_t(\rho)-\widehat\phi(\rho)\,y_{t-1}-X_t\widehat\beta(\rho),
\qquad t=1,...,T,
\end{equation*}
and form $U\in\mathbb R^{T\times N}$ with $t$th row $u_t(\rho,\widehat\vartheta)'$.
The remaining nuisance updates, concentrating out $(\delta,A)$ and the ECM steps for
$(F,\Sigma_\eta,D_\varepsilon)$, proceed exactly as in Steps~1--4.
\subsection{Conditional Scalar Update for \texorpdfstring{$\rho$}{rho}: Outer Loop}
\paragraph{Step 6 (conditional scalar update for $\rho$).}
Since $\rho$ enters nonlinearly through $B(\rho)$ and the Jacobian term $T\log|B(\rho)|$,
there is no closed form for $\rho$.
Because $\rho$ is scalar,
we update it by a one-dimensional conditional maximization step.
This step holds the current values of
$(\phi,\beta,\delta,A,F,\Sigma_\eta,D_\varepsilon)$ fixed.
For any trial value of $\rho$,
form $U(\rho)\in\mathbb R^{T\times N}$ by stacking $u_t(\rho)'$ as its rows,
and define
\begin{equation*}
E(\rho)
=
U(\rho)-\delta\mathbf{1}_{N}'-FAZ,
\end{equation*}
where $u_t(\rho)=B(\rho)y_t-\phi y_{t-1}-X_t\beta$.
Restrict the numerical search to a compact interval
$[\rho_{\min},\rho_{\max}]$ on which $B(\rho)=I_N-\rho W$ is nonsingular.
Up to terms that do not depend on $\rho$ within this conditional step,
the criterion to maximize is
\begin{equation}\label{eq:rho_profile_obj}
\ell_c^{\text{cond}}(\rho)
=
T\log|B(\rho)|
-\frac12\mathrm{tr}\big(ME(\rho)E(\rho)'\big).
\end{equation}
Maximize \eqref{eq:rho_profile_obj} over $[\rho_{\min},\rho_{\max}]$ by a standard one-dimensional routine,
e.g., golden-section/Brent's method, or by a coarse grid followed by local refinement.
After updating $\rho$,
rerun the inner loop to update the remaining parameters.
\begin{remark}[Computing $\log|B(\rho)|$]
Throughout,
$\log|B(\rho)|$ denotes the log absolute determinant $\log|\det B(\rho)|$.
If the eigenvalues of $W$ are $\{\lambda_j\}_{j=1}^N$,
possibly including complex conjugate pairs when $W$ is nonsymmetric,
then
\begin{equation*}
\log|\det B(\rho)|
=
\sum_{j=1}^N\log|1-\rho\lambda_j|,
\quad
\frac{\partial}{\partial\rho}\log|\det B(\rho)|
=
-\operatorname{Re}\sum_{j=1}^N \frac{\lambda_j}{1-\rho\lambda_j}
=
-\mathrm{tr}\big(B(\rho)^{-1}W\big),
\end{equation*}
on any connected admissible region on which $B(\rho)$ is nonsingular.
When $W$ has a real spectrum and the admissible region ensures $1-\rho\lambda_j>0$ for every $j$,
the first expression reduces to $\sum_j\log(1-\rho\lambda_j)$.
For a nonsymmetric $W$,
conjugate eigenvalue pairs combine to give a real log determinant.
The derivative is useful for derivative-based updates in $\rho$,
and the eigenvalue representation makes evaluating $\ell_c^{\text{cond}}(\rho)$ fast.
When feasible, the eigenvalues $\{\lambda_j\}_{j=1}^N$ can be computed once
and stored for fast evaluation of $\log|B(\rho)|$ and its derivative at any $\rho$.
For very large $N$,
the log absolute determinant can instead be evaluated by sparse LU factorization or trace approximations.
\end{remark}
\subsection{Summary of the Estimation Algorithm}
Algorithm~\ref{alg:profile_rho_description} summarizes the inner-outer loop:
the outer loop performs a conditional scalar maximization in $\rho$,
while the inner loop alternates between GLS updates of $(\phi,\beta)$
and ECM-style updates of the nuisance parameters $(\delta,A,F,\Sigma_\eta,D_\varepsilon)$ given the current innovations.
One possible initialization strategy is discussed in Section~\ref{sec:initialization}.
\begin{algorithm}[h]
\caption{Block coordinate estimation with a conditional outer search over $\rho$}\label{alg:profile_rho_description}
\begin{algorithmic}[1]
\Require Initial value for $\rho$, starting values for $(\phi,\beta)$, and starting values for $(\delta,A,F,\Sigma_\eta,D_\varepsilon)$
\Repeat
\State Given current $\rho$, form transformed innovation matrix $U(\rho,\phi,\beta)$
\Repeat
\State Update $(\delta,A)$ by concentration given the current factor path and covariance parameters
\State Compute the conditional moments of the loading heterogeneity
\State Update $(\Sigma_\eta,D_\varepsilon)$
\State Update factor path $F$ and renormalize
\State Update $(\phi,\beta)$ by GLS
\Until{inner convergence}
\State Update $\rho$ by maximizing the conditional criterion over the admissible parameter space,
holding all current non-$\rho$ parameters fixed
\Until{outer convergence}
\State \Return $(\widehat\rho,\widehat\phi,\widehat\beta,\widehat\delta,\widehat A,\widehat F,\widehat\Sigma_\eta,\widehat D_\varepsilon)$
\end{algorithmic}
\end{algorithm}
\section{Model Extensions}\label{sec:extensions}
The baseline formulation is deliberately parsimonious, but two additions are
especially relevant in applications: a lagged spatial outcome and a spatial
autoregressive error. This section shows that both can be accommodated without
changing the fixed-$T$ logic of the estimator.
It also introduces an unrestricted second-order QML that is useful when
the outcome and error filters are difficult to distinguish empirically.
\subsection{Lagged Spatial Outcome}\label{subsec:extension-lagged-spatial-outcome}
Add $Wy_{t-1}$ to the structural equation:
\begin{equation}\label{eq:extension-spatiotemporal}
B(\rho)y_t
=
\phi y_{t-1}+\xi Wy_{t-1}+X_t\beta
+\mathbf{1}_{N}\delta_t+\Lambda f_t+\varepsilon_t.
\end{equation}
Conditional on $y_0$, the Jacobian from $(y_1',...,y_T')'$ to the
innovations is block lower triangular, with $B(\rho)$ on every diagonal block
and $-(\phi I_N+\xi W)$ on the first subdiagonal. Hence its determinant remains
$|B(\rho)|^T$; the lagged spatial outcome changes the transformed innovation,
but contributes no additional determinant term. The criterion in
\eqref{eq:ell_N} therefore remains valid after replacing $u_t(\theta)$ by
\begin{equation*}
u_t(\theta_\xi)
=
B(\rho)y_t-\phi y_{t-1}-\xi Wy_{t-1}-X_t\beta,
\quad\text{where}\quad
\theta_\xi=(\rho,\phi,\xi,\beta')'.
\end{equation*}
For fixed $\rho$ and nuisance parameters, the GLS block simply adds
$Wy_{t-1}$ to the design matrix. Because the $t=1$ equation contains $Wy_0$,
the projection controls should also contain the corresponding information; the
$L=1$ specification in \eqref{eq:correct_projection} does so directly.
The dynamic reduced form is stable when
\begin{equation}\label{eq:extension-spatiotemporal-stability}
\varrho\!\left[B(\rho)^{-1}(\phi I_N+\xi W)\right]<1,
\end{equation}
where $\varrho(\cdot)$ denotes spectral radius. For asymptotic analysis, the
inverse of the associated unit-stacked dynamic-spatial operator must also have
uniformly bounded row and column sums. Under these conditions and the analogues
of Assumptions~\ref{ass:SERL}--\ref{ass:score_clt}, the
proof of Theorem~\ref{thm:qml_asymptotics} applies to $\theta_\xi$: the estimator
is consistent and $\sqrt N$-asymptotically normal, with the same profiled
sandwich form as in \eqref{eq:Vtheta_hat}. Appendix~\ref{sec:extension-scores}
gives the additional corrected score contribution and the modified reduced-form
response blocks.
\subsection{Spatial Errors and the Second-Order QML}
\label{subsec:extension-spatial-errors}
Consider an outcome equation with a spatial autoregressive error,
\begin{align}
B(\rho)y_t
&=
\phi y_{t-1}+\xi Wy_{t-1}+X_t\beta+WX_t\gamma
+\mathbf{1}_{N}\delta_t+v_t,
\label{eq:extension-spatial-error-outcome}\\
B(\psi)v_t&=\Lambda f_t+\varepsilon_t.
\label{eq:extension-spatial-error-process}
\end{align}
\begin{remark}[Placement of Spatial Error Filter]
Specification~\eqref{eq:extension-spatial-error-process} treats the common
and idiosyncratic innovations as entering before spatial propagation, so
$B(\psi)^{-1}$ applies to both. An alternative specification places
$\Lambda f_t$ outside the separate spatial-error process. This alternative placement is studied by
\textcite{li-yang-2021}.
We maintain
\eqref{eq:extension-spatial-error-process}. Our specification treats the composite disturbance $v_t$
itself as following a spatial autoregressive process, in parallel with the spatial autoregressive specification for the outcome.
In addition, under the alternative placement,
non-Gaussian robust inference generally requires Li and Yang's hybrid
higher-moment correction rather than inference based solely on the score rearrangement used here.
\end{remark}
Applying the error filter $B(\psi)$ to the outcome equation gives the product
$B(\psi)B(\rho)$. This product is a polynomial of degree two in $W$.
To introduce a common notation for its two roots, initially label
$\kappa_1:=\rho$ and $\kappa_2:=\psi$, and write
$\kappa:=(\kappa_1,\kappa_2)'$.
Because the two filter components commute,
their ordering is immaterial.
Write
\begin{equation}\label{eq:extension-second-order-filter}
B_*(\kappa)y_t
:=(I_N-\kappa_1W)(I_N-\kappa_2W)y_t
=(I_N-\rho_1W-\rho_2W^2)y_t,
\end{equation}
where $\rho_1=\kappa_1+\kappa_2$ and
$\rho_2=-\kappa_1\kappa_2$. Here $\rho$ remains the outcome spatial
coefficient from the baseline model, whereas the subscripted coefficients
$(\rho_1,\rho_2)$ describe the product polynomial. We refer to the resulting
filter as \emph{second-order} because it contains both $W$ and $W^2$.
Filtering of the spatial error also generates $W^2y_{t-1}$ and $W^2X_t$ terms.
Instead of imposing at the outset the nonlinear restrictions that link all of these coefficients
to $\psi$, we consider the more general transformed equation
\begin{equation}\label{eq:extension-second-order-model}
B_*(\kappa)y_t
=
(\phi I_N+\xi_1W+\xi_2W^2)y_{t-1}
+X_t\beta_0+WX_t\beta_1+W^2X_t\beta_2
+\mathbf{1}_{N}\delta_t+\Lambda f_t+\varepsilon_t.
\end{equation}
We call \eqref{eq:extension-second-order-model} the \emph{unrestricted
second-order model}: ``second-order'' refers to the inclusion of $W^2$, while
``unrestricted'' means that the coefficients on the $I_N$, $W$, and $W^2$
terms are estimated separately rather than constrained to originate from
\eqref{eq:extension-spatial-error-outcome}--\eqref{eq:extension-spatial-error-process}.
The structural spatial-error model is a restricted special case. Expanding
$B(\psi)$ times the right-hand side of
\eqref{eq:extension-spatial-error-outcome} gives
\begin{equation}\label{eq:extension-spatial-error-restrictions}
\begin{aligned}
\{\kappa_1,\kappa_2\}&=\{\rho,\psi\},
&\xi_1&=\xi-\psi\phi,
&\xi_2&=-\psi \xi,\\
\beta_0&=\beta,
&\beta_1&=\gamma-\psi\beta,
&\beta_2&=-\psi\gamma.
\end{aligned}
\end{equation}
Moreover,
$B(\psi)\mathbf{1}_{N}\delta_t=(1-\psi)\mathbf{1}_{N}\delta_t$ because
$W\mathbf{1}_{N}=\mathbf{1}_{N}$; this rescaling is absorbed into the freely estimated time
effect in \eqref{eq:extension-second-order-model}. The factor innovation
$\Lambda f_t+\varepsilon_t$ therefore retains the baseline
low-rank-plus-diagonal covariance structure.
The normalized conditional criterion for \eqref{eq:extension-second-order-model}
is
\begin{equation}\label{eq:extension-second-order-likelihood}
\ell_{N,2}(\alpha_2)
=
\frac{T}{N}\sum_{j=1}^2\log|B(\kappa_j)|
-\frac{1}{2N}\sum_{i=1}^N
\left[\log|\Sigma_u|+e_{2i}'\Sigma_u^{-1}e_{2i}\right],
\end{equation}
where $e_{2i}$ is obtained from the transformed innovation in
\eqref{eq:extension-second-order-model}. Conditional on $(\kappa_1,\kappa_2)$,
all dynamic and control coefficients enter the same GLS block as before; the
two log-determinant terms correspond to the two components of $B_*(\kappa)$.
Thus, incorporating the spatial-error filter replaces the baseline scalar
outer search over $\rho$ with a two-dimensional optimization over
$(\kappa_1,\kappa_2)$, while leaving the inner GLS profiling steps unchanged.
Dynamic stability requires
\begin{equation}\label{eq:extension-second-order-stability}
\varrho\!\left[
B_*(\kappa)^{-1}(\phi I_N+\xi_1W+\xi_2W^2)
\right]<1.
\end{equation}
The criterion is unchanged when $\kappa_1$ and $\kappa_2$ are interchanged.
Without imposing
\eqref{eq:extension-spatial-error-restrictions}, it does not identify which
root is the outcome coefficient $\rho$ and which is the error coefficient
$\psi$; only the symmetric filter coefficients $(\rho_1,\rho_2)$ are invariant
to relabeling. This distinction is useful in short panels, where direct
separation of $\rho$ and $\psi$ may be weak.
For theoretical purposes,
we label the roots by imposing $\kappa_1<\kappa_2$ and assume that the true root gap $\kappa_{20}-\kappa_{10}$ is bounded away from zero.
This ordering makes the root parameterization locally unique but does not assign structural roles to the two roots.
The equal-root case is excluded because the transformation
$(\kappa_1,\kappa_2)\mapsto(\rho_1,\rho_2)$ has Jacobian determinant $\kappa_2-\kappa_1$
and is therefore singular when the roots coincide.
Suppose this root-separation condition holds,
the ordered roots lie in the interior of the admissible set,
and the bounded-inverse, conditional-moment, smoothness, nonsingular
Hessian, and score-CLT conditions used for
Theorem~\ref{thm:qml_asymptotics} hold for the enlarged operator. Then the QML
estimator of the ordered-root parameter vector
$(\kappa_1,\kappa_2,\phi,\xi_1,\xi_2,\beta_0',\beta_1',\beta_2')'$ is consistent
and $\sqrt N$-asymptotically normal, with the profiled sandwich covariance in
\eqref{eq:Vtheta_hat}. Inference for the reported symmetric coefficients
$(\rho_1,\rho_2)$ follows by the delta method. The spatial score decomposition
must use the second-order reduced form for every endogenous direction;
Appendix~\ref{sec:extension-scores} records this modification. The same argument
covers the restricted spatial-error QML when the restrictions in
\eqref{eq:extension-spatial-error-restrictions} are imposed and locally
identified.
\section{Simulation}\label{sec:simulation}
We report two complementary experiments based on the model in
\eqref{eq:ysystem}. The first evaluates estimation and inference over a broad
parameter and network grid when the initial outcome is generated directly from
the factor loadings. The second considers the more demanding case in which the
initial outcome is inherited from a long-running spatial process and compares
alternative spatial enrichment orders in the SERL specification.
All estimators and covariance matrices reported in the Monte Carlo experiments are computed using
the companion R package \href{https://github.com/jessekelighine/dspserl}{\texttt{dspserl}}.
\subsection{Baseline Estimation and Inference}
\label{subsec:simulation-baseline}
\subsubsection{Design}
The baseline experiment has $N$ units, $T$ time periods, $K$ regressors, and
$r$ common factors. The spatial weights matrix $W$ is a row-normalized
$k$-nearest-neighbor graph constructed from random locations, with $k=8$.
The network is redrawn in each Monte Carlo replication using a seed separate
from the remaining DGP draws and is fixed over time within that replication.
The main-text experiment stores $W$ as an ordinary $N\times N$ array;
we refer to this as a \emph{dense-matrix} design.
This contrasts with the \emph{sparse-matrix} design presented in Appendix~\ref{subsec:aux-sparse-logdet},
which stores $W$ as a sparse matrix and utilizes sparse-matrix methods for the log-determinant.
For each $t$, the regressors are i.i.d.\ $x_{it}\sim\mathcal N(0,I_K)$ and
$\bar x_i$ denotes their time average. We generate loadings directly as
$\lambda_i\sim\mathcal N(0,I_r)$ and set the initial condition by
\begin{equation*}
y_{i0}=\mathbf{1}_{r}'\lambda_i+\nu_i,
\end{equation*}
where $\nu_i\sim\mathcal N(0,1)$ is independent of $\lambda_i$.
Thus $y_{i0}$ is informative about the conditional mean of $\lambda_i$, while
$\bar x_i$ remains independent of the loadings.
Joint Gaussianity gives
\begin{equation*}
\mathbb E(\lambda_i\mid y_{i0})
=
\frac{\mathbf{1}_{r}}{r+1}y_{i0},
\qquad
\mathrm{Var}(\lambda_i\mid y_{i0})
=
I_r-\frac{\mathbf{1}_{r}\mathbf{1}_{r}'}{r+1}.
\end{equation*}
Time effects $\delta_t$ follow a deterministic cycle,
idiosyncratic variances are time-varying,
and common shocks are generated as $f_t\sim\mathcal N(0,I_r)$.
We set $z_i=(\bar x_i',y_{i0})'$ in the SERL specification \eqref{eq:proj}. In
this DGP, population projection parameters are
$A
=
\left[0_{r\times K},\frac{\mathbf{1}_{r}}{r+1}\right],$ and
$\Sigma_\eta
=
I_r-\frac{\mathbf{1}_{r}\mathbf{1}_{r}'}{r+1}.
$
Both are nuisance parameters implied by the DGP.
Given $(W,\rho)$, we solve for $y_t$ using the spatial filter
$B(\rho)=I_N-\rho W$:
\begin{equation}
y_t=B(\rho)^{-1}\big(\phi y_{t-1}+X_t\beta+\mathbf{1}_{N}\delta_t+\Lambda f_t+\varepsilon_t\big),
\qquad t=1,...,T,
\end{equation}
with $y_{i0}$ initialized as above and
$\varepsilon_t\sim\mathcal N(0,\sigma_t^2 I_N)$.
The implementation factors $B(\rho)$ once by LU decomposition and reuses that
factorization to solve the system in every period; it does not form
$B(\rho)^{-1}$ explicitly.
We fix $r=2$ and $K=2$ with $\beta=(0.8,-0.3)$.
To hold persistence associated with the largest eigenvalue of the row-normalized
$W$ constant across values of $\rho$, we set
\begin{equation}
\phi=0.5(1-\rho).
\end{equation}
Because the largest eigenvalue of $W$ is one,
the corresponding dynamic coefficient is $\phi/(1-\rho)=0.5$
in every design.
Accordingly, $\rho+\phi$ equals $0.60$, $0.75$, or $0.90$
and remains below one.
The time effects and idiosyncratic variances follow deterministic schedules,
\begin{equation}
\delta_t=0.3\sin\left(2\pi \frac{t}{T}\right),
\qquad
\sigma_t^2=5+5 \frac{t}{T}.
\end{equation}
Common shocks are $f_t\sim\mathcal N(0,I_r)$.
For the bias and dispersion experiment, we vary $(N,T,\rho)$ according to the grid
\begin{equation}
\begin{aligned}
T&\in\{5,10,20\}, & N&\in\{500,1000\},\\
\rho&\in\{0.2,0.5,0.8\},
& \phi&=0.5(1-\rho).
\end{aligned}
\end{equation}
for a total of $18$ parameter configurations.
For each configuration we estimate $(\rho,\phi,\beta)$ by the block-coordinate QML procedure.
We evaluate $\log|B(\rho)|$ using the eigenvalues of $W$.
Appendix~\ref{subsec:aux-sparse-logdet} repeats these $18$ configurations
using a 30-term Hutchinson trace approximation based on 25 Rademacher vectors.
We also compare this approximation with a richer trace calculation
and exact sparse LU factorization on identical panels.
Both experiments initialize the spatial search with the same truth-independent
nine-point grid over $[-0.95,0.95]$.
The inner loop permits at most 200 iterations with tolerance $10^{-8}$,
and the scalar optimizer for $\rho$ uses tolerance $10^{-10}$.
Each design requests 1000 Monte Carlo replications.
For this experiment, the outer block-coordinate loop permits at most 100
iterations.
A fit enters the bias and dispersion calculations if it either converges
strictly or has a converged inner loop and a final $\rho$ stationarity gap no
larger than $2\times10^{-6}$.
Between 995 and 1000 estimates per design meet this criterion.
The coverage experiment uses the same parameter grid.
Each design requests 1000 Monte Carlo replications.
The outer block-coordinate loop permits at most 300 iterations and requires
the strict $10^{-8}$ convergence criterion.
For each strictly converged fit,
we construct nominal $95\%$ intervals using the martingale-sample covariance estimator.
We report the uncorrected version as a benchmark and the pair-only HC1 correction in \eqref{eq:Omega_martingale_pair_hc1} as the preferred estimator.
\subsubsection{Results}
Table~\ref{tab:sim-bias-sd-dense} reports bias, multiplied by $100$,
and the standard deviation (SD) for $\rho$, $\phi$, and $\beta$
across the 18 dense-matrix designs.
The biases are generally small,
although $\hat\rho$ has a systematic downward bias.
In the original parameter units,
its average bias is $-0.0032$ and its largest absolute bias is about $0.0051$.
The SD generally decreases with larger $N$ and $T$,
and $\rho$ is estimated more precisely when its true value is larger,
consistent with stronger spatial dependence providing more information about
$\rho$. In particular, the SD of $\hat\rho$ decreases with $N$ and $T$ when
$\rho=0.8$.
\begin{sidewaystable}[t]
\centering
\scriptsize
\csvreader[
respect underscore=true,
column names={
T=\Tval,
N=\Nval,
knn_type=\knntype,
rho=\truerho,
phi=\truephi,
beta1=\truebetaone,
beta2=\truebetatwo,
bias_x100_rho=\BIASrho,
sd_rho=\SDrho,
bias_x100_phi=\BIASphi,
sd_phi=\SDphi,
bias_x100_beta1=\BIASbetaone,
sd_beta1=\SDbetaone,
bias_x100_beta2=\BIASbetatwo,
sd_beta2=\SDbetatwo,
},
filter strcmp={\knntype}{dense},
tabular=rrrrrr|rrrrrrrr,
table head=
\toprule
\multicolumn{1}{c}{$T$} &
\multicolumn{1}{c}{$N$} &
\multicolumn{1}{c}{$\rho^{\text{true}}$} &
\multicolumn{1}{c}{$\phi^{\text{true}}$} &
\multicolumn{1}{c}{$\beta_1^{\text{true}}$} &
\multicolumn{1}{c}{$\beta_2^{\text{true}}$} &
\multicolumn{2}{c}{$\hat\rho$} &
\multicolumn{2}{c}{$\hat\phi$} &
\multicolumn{2}{c}{$\hat\beta_1$} &
\multicolumn{2}{c}{$\hat\beta_2$} \\
\cmidrule(lr){7-8}\cmidrule(lr){9-10}
\cmidrule(lr){11-12}\cmidrule(lr){13-14}
&&&&&&
\multicolumn{1}{c}{$100\times\mathrm{Bias}$} &
\multicolumn{1}{c}{SD} &
\multicolumn{1}{c}{$100\times\mathrm{Bias}$} &
\multicolumn{1}{c}{SD} &
\multicolumn{1}{c}{$100\times\mathrm{Bias}$} &
\multicolumn{1}{c}{SD} &
\multicolumn{1}{c}{$100\times\mathrm{Bias}$} &
\multicolumn{1}{c}{SD} \\
\midrule,
table foot = \bottomrule
]{code/output/simulation_bias_sd.csv}{}{
\Tval & \Nval &
\num[group-digits=false,round-precision=1]{\truerho} &
\num[group-digits=false,round-precision=1]{\truephi} &
\num[group-digits=false,round-precision=1]{\truebetaone} &
\num[group-digits=false,round-precision=1]{\truebetatwo} &
\num[group-digits=false,round-precision=4]{\BIASrho} &
\num[group-digits=false,round-precision=5]{\SDrho} &
\num[group-digits=false,round-precision=4]{\BIASphi} &
\num[group-digits=false,round-precision=5]{\SDphi} &
\num[group-digits=false,round-precision=4]{\BIASbetaone} &
\num[group-digits=false,round-precision=5]{\SDbetaone} &
\num[group-digits=false,round-precision=4]{\BIASbetatwo} &
\num[group-digits=false,round-precision=5]{\SDbetatwo}
}
\caption{\textsc{Dense-matrix Design}.
Bias and standard deviations.
Each row is defined by $T$, $N$, and the true parameter values.
The remaining columns report Monte Carlo bias, multiplied by 100,
and the unscaled sample standard deviation of $\hat\rho$, $\hat\phi$,
$\hat\beta_1$, and $\hat\beta_2$.
Each of the 18 designs requests 1000 replications;
a fit is accepted if it either converges strictly or has a converged inner loop and a final $\rho$ stationarity gap no larger than $2\times10^{-6}$.
Between 995 and 1000 estimates per design enter the calculations,
for 17,984 accepted estimates in total.
The log determinant is evaluated from the eigenvalues of $W$.}
\label{tab:sim-bias-sd-dense}
\end{sidewaystable}
Table~\ref{tab:sim-coverage-dense} reports coverage for the dense-matrix designs without a finite-sample correction
and with the preferred pair-only HC1 correction in \eqref{eq:Omega_martingale_pair_hc1}.
Coverages for the regression slopes and $\rho$ are close to the nominal level,
while coverage for $\phi$ is lowest when $T=5$.
Appendix~\ref{subsec:aux-symmetric} holds $\rho+\phi$ fixed while varying its spatial and dynamic components.
The experiment confirms that first-order Wald coverage for $\phi$ can be inaccurate when each unit contributes only five transitions,
whereas increasing the time dimension to $T=20$ eliminates under-coverage.
We view this as an understandable short-panel limitation rather than a failure of the large-$N$ theory.
The pair-only HC1 correction modestly improves coverage for $\rho$ without materially changing the other results,
so we use it as our preferred finite-sample covariance estimator and retain the uncorrected columns as a benchmark.
We also conduct the complete experiment under the sparse-matrix design.
Appendix~\ref{subsec:aux-sparse-logdet} reports the corresponding bias, dispersion, and coverage results.
The parameters remain well estimated and dispersion is similar to that reported in Table~\ref{tab:sim-bias-sd-dense},
but coverage for $\rho$ is generally lower.
The diagnostic using identical panels shows that much of this difference is attributable to approximation error
in the log determinant.
\begin{sidewaystable}[t]
\centering
\scriptsize
\csvreader[
respect underscore=true,
column names={
2=\Tval,
3=\Nval,
4=\knntype,
5=\truerho,
6=\truephi,
7=\truebetaone,
8=\truebetatwo,
9=\uncorrectedrho,
10=\uncorrectedphi,
11=\hcspatialrho,
12=\hcspatialphi,
13=\spatialbetaone,
14=\spatialbetatwo
},
filter strcmp={\knntype}{"dense"},
tabular=rr|rrrr|rr|rr|rr,
table head=
\toprule
&&&&&&
\multicolumn{2}{c}{No correction} &
\multicolumn{2}{c}{Pair-only HC1} &
\multicolumn{2}{c}{Regression slopes} \\
\cmidrule(lr){7-8}\cmidrule(lr){9-10}\cmidrule(lr){11-12}
$T$ & $N$ &
$\rho^{\mathrm{true}}$ & $\phi^{\mathrm{true}}$ &
$\beta_1^{\mathrm{true}}$ & $\beta_2^{\mathrm{true}}$ &
$\rho$ & $\phi$ & $\rho$ & $\phi$ & $\beta_1$ & $\beta_2$ \\
\midrule,
table foot=\bottomrule
]{code/output/simulation_coverage_spatial_comparison.csv}{}{
\Tval & \Nval &
\num[group-digits=false,round-precision=1]{\truerho} &
\num[group-digits=false,round-precision=1]{\truephi} &
\num[group-digits=false,round-precision=1]{\truebetaone} &
\num[group-digits=false,round-precision=1]{\truebetatwo} &
\num[group-digits=false,round-precision=3]{\uncorrectedrho} &
\num[group-digits=false,round-precision=3]{\uncorrectedphi} &
\num[group-digits=false,round-precision=3]{\hcspatialrho} &
\num[group-digits=false,round-precision=3]{\hcspatialphi} &
\num[group-digits=false,round-precision=3]{\spatialbetaone} &
\num[group-digits=false,round-precision=3]{\spatialbetatwo}
}
\caption{\textsc{Dense-matrix Design}.
Coverage of nominal $95\%$ confidence intervals.
The first two coverage columns use the martingale-sample covariance estimator in \eqref{eq:martingale_covariance_estimator} without a finite-sample correction,
while the next two use the pair-only HC1 correction in \eqref{eq:Omega_martingale_pair_hc1}.
Because the correction leaves coverage for $\beta_1$ and $\beta_2$ unchanged in every design,
the regression-slope coverage is reported once.
Each of the 18 designs requests 1000 Monte Carlo replications;
between 979 and 1000 replications per design produce strictly converged estimates and complete sandwich inference,
for 17,930 successful replications in total.
The log determinant is evaluated from the eigenvalues of $W$.}
\label{tab:sim-coverage-dense}
\end{sidewaystable}
\subsection{Long-Running Initial Conditions and SERL Enrichment}
\label{subsec:simulation-initial-condition-projection}
The preceding experiments generate $y_0$ directly from the factor loadings.
We now consider the more demanding case in which $y_0$ is itself an outcome
from a long-running dynamic spatial process. This design evaluates the spatial
enrichment in \eqref{eq:correct_projection}. In particular, it asks whether a
small number of spatial transformations of the initial condition and regressor
summaries can capture the information about the loadings that propagates
through the spatial multiplier.
\subsubsection{Design}
For each replication, we draw $\lambda_i\sim\mathcal N(0,I_2)$ and generate a
stationary regressor component according to
\begin{equation*}
h_{is}=0.8h_{i,s-1}+\sqrt{1-0.8^2}\,v_{is},
\qquad v_{is}\sim\mathcal N(0,I_2).
\end{equation*}
The observed regressors are
\begin{equation*}
x_{is}=0.5\lambda_i+\sqrt{1-0.5^2}\,h_{is},
\end{equation*}
so that both the initial outcome and the regressors are informative about the
individual loadings. This makes the projection exercise more demanding than a
design in which the loadings are related only to $y_0$.
The regressor process is initialized from its stationary distribution.
Starting from zero,
we run the outcome process for 100 presample periods and take its terminal
presample value as $y_0$.
We then retain the next $T$ observations as the estimation sample.
The structural parameters are $\phi=0.3$ and $\beta=(0.8,-0.3)'$,
the idiosyncratic innovations have unit variance,
and the time effects are zero.
There are two common factors whose realizations are independently standard
normal over the presample and sample periods.
We use a sparse row-normalized $k$-nearest-neighbor network with $k=8$ and
consider
\begin{equation*}
T\in\{5,10\},\qquad N\in\{500,1000\},\qquad
\rho\in\{0.2,0.5\}.
\end{equation*}
Unlike the experiments initialized directly at $y_0$, a long-running process
must satisfy a stability restriction during the presample recursion. For a
row-normalized $W$, the leading persistence is $\phi/(1-\rho)$,
which equals $0.375$ and $0.6$ for $\rho=0.2$ and $\rho=0.5$, respectively.
For each $N$, the network is held fixed across replications and shared across
the corresponding $(T,\rho)$ designs.
For each $T$, the factor path is held fixed across replications and shared
across the corresponding $(N,\rho)$ designs.
Each of the eight designs contains 500 Monte Carlo panels.
Let $Q_X$ contain the unit-specific time averages of the two observed
regressors. The full projection of order $L$ uses
\begin{equation*}
z_i^{(L)}=\big[(Q_X)_i',(WQ_X)_i',...,(W^LQ_X)_i',
(y_0)_i,(Wy_0)_i,...,(W^Ly_0)_i\big]',
\end{equation*}
and we consider $L=0,1,2,3$. These specifications contain, respectively,
$3$, $6$, $9$, and $12$ controls. We also consider a parsimonious $y$-only
$L=2$ specification that includes $Q_X$, $y_0$, $Wy_0$, and $W^2y_0$, but
does not include spatial transformations of $Q_X$. The same simulated panel
is used for all five projection specifications, so comparisons across them are
paired.
Each control is centered and standardized across units before estimation;
this rescaling does not change the span of the projection.
The log determinant is evaluated exactly by sparse LU factorization.
This avoids confounding the projection comparison with the stochastic approximation error documented in Appendix~\ref{subsec:aux-sparse-logdet}.
Inference uses the martingale-sample covariance estimator with the preferred pair-only HC1 correction in \eqref{eq:Omega_martingale_pair_hc1}.
\subsubsection{Results}
Table~\ref{tab:sim-initial-projection-summary} summarizes projection quality
and estimator performance across the eight designs.
Under
the conventional $L=0$ projection, the average maximum absolute correlation
between the loading residual and the first omitted spatial transformation of
$y_0$ is $0.167$, while the corresponding diagnostic for the regressor
summaries is $0.116$. Adding the first spatial lag reduces these diagnostics to
$0.018$ and $0.032$, respectively.
It also reduces the RMSE of $\hat\rho$ from $0.0281$ to $0.0204$ and the RMSE of $\hat\phi$ from $0.0238$ to $0.0181$.
Average coverage rises from $0.786$ to $0.944$ for $\rho$ and from $0.852$ to $0.939$ for $\phi$.
Coverage of $\beta_1$ rises from $0.910$ to $0.931$, while coverage of $\beta_2$ rises from $0.941$ to $0.947$.
Higher-order spatial transformations continue to reduce the omitted-variable
diagnostics, but they provide no further improvement in estimation. Relative
to full $L=1$, the full $L=2$ and $L=3$ projections have slightly larger bias
and RMSE for $\rho$ and $\phi$, consistent with a modest finite-sample cost
from estimating additional projection coefficients. The $y$-only $L=2$
projection performs similarly to full $L=1$ and has slightly smaller aggregate
bias and RMSE. It nevertheless leaves more residual correlation with omitted
spatial transformations of $Q_X$ than the full projection. Thus, the spatial
transformations of $y_0$ account for most of the estimation gain in this DGP,
while the spatial transformations of $Q_X$ improve the quality of the loading
projection without materially changing the structural estimates.
Table~\ref{tab:sim-initial-projection-coverage} reports design-specific
coverage. The benefit of the richer projection is largest when spatial
dependence is strong. Averaging over the four designs with $\rho=0.5$, moving
from full $L=0$ to full $L=1$ raises coverage from $0.665$ to $0.944$ for
$\rho$ and from $0.768$ to $0.937$ for $\phi$.
For $(T,N,\rho)=(5,500,0.5)$,
the corresponding changes are from $0.642$ to $0.918$ and from $0.750$ to $0.936$.
For the most difficult design,
$(T,N,\rho)=(5,1000,0.5)$, coverage rises from $0.312$ to $0.935$ for $\rho$
and from $0.496$ to $0.910$ for $\phi$.
Overall, moving from $L=0$ to $L=1$ substantially improves coverage
and supports the practical value of low-order SERL enrichment,
particularly when spatial dependence is stronger.
\begin{sidewaystable}[t]
\centering
\csvreader[
respect underscore=true,
column names={
1=\projectionid,
3=\ncontrols,
7=\corry,
8=\corrx,
9=\biasrho,
10=\rmserho,
11=\coveragerho,
12=\biasphi,
13=\rmsephi,
14=\coveragephi,
17=\coveragebetaone,
20=\coveragebetatwo
},
tabular=lrrr|rrr|rrr|rr,
table head=
\toprule
&&&& \multicolumn{3}{c}{$\rho$} &
\multicolumn{3}{c}{$\phi$} &
\multicolumn{2}{c}{Coverage} \\
\cmidrule(lr){5-7}\cmidrule(lr){8-10}\cmidrule(lr){11-12}
Projection & $q$ & $\operatorname{Corr}_y$ &
$\operatorname{Corr}_x$ & Bias & RMSE & Coverage &
Bias & RMSE & Coverage & $\beta_1$ & $\beta_2$ \\
\midrule,
table foot=\bottomrule
]{code/output/simulation_initial_condition_projection_table_summary.csv}{}{
\ifcase\projectionid Full $L=0$
\or Full $L=1$
\or Full $L=2$
\or Full $L=3$
\or $y$-only $L=2$
\fi & \ncontrols &
\num[round-precision=3]{\corry} &
\num[round-precision=3]{\corrx} &
\num[round-precision=3]{\biasrho} &
\num[round-precision=3]{\rmserho} &
\num[round-precision=3]{\coveragerho} &
\num[round-precision=3]{\biasphi} &
\num[round-precision=3]{\rmsephi} &
\num[round-precision=3]{\coveragephi} &
\num[round-precision=3]{\coveragebetaone} &
\num[round-precision=3]{\coveragebetatwo}
}
\caption{Long-running initial conditions:
projection and estimation performance averaged across 8 network designs.
The full projection includes spatial transformations of both
$Q_X$ and $y_0$ through order $L$; the $y$-only specification includes
$Q_X$ without spatial transformations and spatial transformations of $y_0$
through order two. The column $q$ is the number of controls.
$\operatorname{Corr}_y$ and $\operatorname{Corr}_x$ are the mean maximum
absolute correlations between the loading residuals and the next omitted
spatial transformations of $y_0$ and $Q_X$, respectively.
The regressors contain a factor loading component with coefficient $0.5$.
The eight designs combine $T\in\{5,10\}$, $N\in\{500,1000\}$, and
$\rho\in\{0.2,0.5\}$. Bias and RMSE are computed relative to the true
parameter value in each design and pooled across strictly converged fits.
Each design has 500 Monte Carlo replications.
All designs evaluate the log determinant exactly by sparse LU
factorization.
Inference uses the martingale-sample covariance estimator with the pair-only
HC1 correction. Depending on the projection, 3910--3918 of the 4000 attempted
fits converged strictly and produced spatial-sandwich standard errors.}
\label{tab:sim-initial-projection-summary}
\end{sidewaystable}
\begin{sidewaystable}[t]
\centering
\csvreader[
respect underscore=true,
column names={
2=\Tval,
3=\Nval,
4=\truerho,
5=\Lzerorho,
6=\Lzerophi,
7=\Lonerho,
8=\Lonephi,
9=\Ltworho,
10=\Ltwophi,
11=\Lthreerho,
12=\Lthreephi,
13=\yonlyrho,
14=\yonlyphi
},
tabular=rrr|rr|rr|rr|rr|rr,
table head=
\toprule
&&& \multicolumn{2}{c}{Full $L=0$} &
\multicolumn{2}{c}{Full $L=1$} &
\multicolumn{2}{c}{Full $L=2$} &
\multicolumn{2}{c}{Full $L=3$} &
\multicolumn{2}{c}{$y$-only $L=2$} \\
\cmidrule(lr){4-5}\cmidrule(lr){6-7}\cmidrule(lr){8-9}
\cmidrule(lr){10-11}\cmidrule(lr){12-13}
$T$ & $N$ & $\rho^{\mathrm{true}}$ &
$\rho$ & $\phi$ & $\rho$ & $\phi$ & $\rho$ & $\phi$ &
$\rho$ & $\phi$ & $\rho$ & $\phi$ \\
\midrule,
table foot=\bottomrule
]{code/output/simulation_initial_condition_projection_table_coverage.csv}{}{
\Tval & \Nval & \truerho &
\num[round-precision=3]{\Lzerorho} &
\num[round-precision=3]{\Lzerophi} &
\num[round-precision=3]{\Lonerho} &
\num[round-precision=3]{\Lonephi} &
\num[round-precision=3]{\Ltworho} &
\num[round-precision=3]{\Ltwophi} &
\num[round-precision=3]{\Lthreerho} &
\num[round-precision=3]{\Lthreephi} &
\num[round-precision=3]{\yonlyrho} &
\num[round-precision=3]{\yonlyphi}
}
\caption{Long-running initial conditions: coverage of nominal $95\%$
confidence intervals for $\rho$ and $\phi$. Each design contains 500 Monte
Carlo panels, and the same panel is estimated under all five projection
specifications. The eight dynamically stable designs use
$T\in\{5,10\}$, $N\in\{500,1000\}$, and $\rho\in\{0.2,0.5\}$.
All designs evaluate the log determinant exactly by sparse LU
factorization.
Coverage is computed from the 461--500 strictly converged fits in each cell;
there were no spatial-sandwich failures among those fits.
Inference uses the martingale-sample covariance estimator with the pair-only HC1 correction.
The full and $y$-only projections are defined in the text and summarized in
Table~\ref{tab:sim-initial-projection-summary}.}
\label{tab:sim-initial-projection-coverage}
\end{sidewaystable}
\section{Application: Female Labor-Force Participation}\label{sec:application}
We apply the fixed-$T$ estimator to the county-level female labor-force
participation data studied by \textcite{fogli-veldkamp-2011} and reexamined by
\textcite{tziolas-elhorst-2023}. The data contain 3,074 US counties observed at
decennial intervals from 1940 through 2000. Female labor-force participation
rose from 18.49 percent to 54.69 percent over this period. The observed controls
are the urban and rural-farm population shares, average education, population
density, and manufacturing wages. The economic question is whether changes in
participation diffuse across neighboring counties, potentially with a
one-decade lag.
We estimate all empirical specifications using the companion R package \href{https://github.com/jessekelighine/dspserl}{\texttt{dspserl}}.
\subsection{Empirical Specification and Samples}
Let $y_t$ be the $N\times1$ vector of county female labor-force participation
rates in decade $t$, and let $W$ be a row-normalized spatial-weights matrix, so
that $Wy_t$ contains weighted averages of neighboring counties' participation
rates. Writing $X_t$ for the non-spatial controls, we organize the application
around the encompassing outcome equation
\begin{equation}\label{eq:application-model}
y_t
=
\rho Wy_t+\phi y_{t-1}+\xi Wy_{t-1}
+X_t\beta+WX_t\gamma
+\mathbf{1}_{N}\delta_t+\Lambda f_t+\varepsilon_t.
\end{equation}
The vector $\gamma$ contains the spatial Durbin coefficients on the spatial
lags of the controls. When a spatially lagged control is excluded, its
corresponding element of $\gamma$ is set to zero. The spatially lagged
controls can be stacked with $X_t$ and therefore enter the regressor block of
the baseline model.
The coefficient $\rho$ captures contemporaneous outcome dependence, whereas
$\xi$ captures dependence on neighboring counties' outcomes one decade
earlier. We first impose $\xi=0$, yielding the baseline model analyzed above,
and then estimate $\xi$ freely using the extension in
Section~\ref{subsec:extension-lagged-spatial-outcome}.
The lagged-neighbor term is present in the empirical model of
\textcite{fogli-veldkamp-2011}, who omit the contemporaneous outcome lag.
\textcite{tziolas-elhorst-2023} place both $Wy_t$ and $Wy_{t-1}$ in a broader
dynamic spatial Durbin model that also contains a spatial autoregressive error.
In this empirical application, the common shocks $f_t$ may represent nationwide
developments affecting the benefits and costs of female employment, including
changes in wages and discrimination,
the diffusion of household appliances,
the declining physical demands of jobs,
improved control over fertility,
and changing social norms.
Following \textcite{tziolas-elhorst-2023}, we consider three covariate and
spatial-weight specifications. The first includes wages and spatial lags of all
five controls. Requiring a balanced panel leaves 1,569 counties, and $W$ is the
row-normalized six-nearest-neighbor matrix constructed from their coordinates.
The second excludes wages but retains spatial lags of the other four controls.
The third also removes the spatially lagged controls. The last two
specifications contain 3,066 counties and use the original binary-contiguity
matrix after subsetting and row normalization. Because the 1940 outcome is the
initial condition, each estimation sample has six periods, covering
1950--2000.
All main application tables use two common factors.
This parsimonious choice is natural with only six estimation periods and matches the two-factor specification of \textcite{tziolas-elhorst-2023}.
Appendix~\ref{sec:application-factor-rank} reports three-factor estimates of the second-order model as a robustness check.
With six estimation periods, $r=3$ is the largest factor rank
for which the covariance component is estimable: under the factor
normalization, $r=3$ uses all 21 distinct elements of a $6\times6$ covariance
matrix, whereas $r=4$ would require 24 covariance parameters and is therefore
not feasible. The projection sets $L=1$ in
\eqref{eq:correct_projection}: $z_i$ contains time averages of the non-spatial
controls, their first spatial lags, $y_{i,1940}$, and $(Wy_{1940})_i$. The
controls are centered and standardized. We evaluate the spatial Jacobian
$\log|B(\rho)|$ exactly by sparse LU factorization.
This avoids the additional numerical dispersion from the faster trace approximation documented in Appendix~\ref{subsec:aux-sparse-logdet}.
All standard errors for our application estimates use the martingale-sample covariance estimator with the pair-only HC1 correction in \eqref{eq:Omega_martingale_pair_hc1}.
\subsection{Baseline Model}
The baseline sets $\xi=0$ in \eqref{eq:application-model}. It is therefore an
instance of \eqref{eq:ysystem}, with included spatially lagged controls treated
as additional observed regressors, and is covered by
Theorem~\ref{thm:qml_asymptotics}. Because $W\mathbf{1}_{N}=\mathbf{1}_{N}$, a spatially
uniform component evolves with dynamic coefficient $\phi/(1-\rho)$, and its
long-run multiplier has denominator $1-\rho-\phi$. We therefore call
$s:=\rho+\phi$ the \emph{stability sum} in the baseline model. When $\rho<1$
and $\phi\geq0$, stability of this component is equivalent to $s<1$; full
dynamic stability requires the spectral-radius condition in
\eqref{eq:extension-spatiotemporal-stability}.
To compare the estimated covariate effects with the historical change in
participation, consider a permanent common increase $\Delta x_k$ in control
$k$, holding the other controls, time effects, and common factors fixed. The
change between the initial and new steady states satisfies
\begin{equation*}
\left[(1-\phi)I_N-\rho W\right]\Delta y_k^*
=
(\beta_k I_N+\gamma_k W)\mathbf{1}_{N}\Delta x_k.
\end{equation*}
Thus, premultiplying by $N^{-1}\mathbf{1}_{N}'$ gives the implied change in average
participation. Let $\Delta x_k$ equal the observed 1940--2000 change in the
sample mean of control $k$, and let
$\Delta\mathrm{LFP}=54.69-18.49=36.20$ denote the corresponding percentage-point
increase in mean female labor-force participation. Dividing the predicted
average change by $\Delta\mathrm{LFP}$ gives the reported contribution:
\begin{equation*}
100\times\frac{1}{N}\mathbf{1}_{N}'
\left[(1-\phi)I_N-\rho W\right]^{-1}
(\beta_k I_N+\gamma_k W)\mathbf{1}_{N}
\frac{\Delta x_k}{\Delta\mathrm{LFP}}.
\end{equation*}
The changes $\Delta x_k$ are those reported by
\textcite{tziolas-elhorst-2023}.
\begin{table}[htbp]
\centering
\caption{Female labor-force participation: baseline QMLE estimates}
\label{tab:application-female-lfp-baseline}
\scriptsize
\setlength{\tabcolsep}{3pt}
\begin{tabular}{lrrrrrr}
\toprule
& \multicolumn{2}{c}{Wages and $WX$} & \multicolumn{2}{c}{No wages; with $WX$} & \multicolumn{2}{c}{No wages or $WX$} \\
\cmidrule(lr){2-3}\cmidrule(lr){4-5}\cmidrule(lr){6-7}
Variable & Estimate (SE) & Contribution & Estimate (SE) & Contribution & Estimate (SE) & Contribution \\
\midrule
$y_{i,t-1}$ & 0.372 (0.019) & -- & 0.467 (0.035) & -- & 0.364 (0.016) & -- \\
$Wy_t$ & 0.521 (0.017) & -- & 0.479 (0.016) & -- & 0.506 (0.011) & -- \\
Urban population & 0.011 (0.005) & -10.00 & 0.017 (0.005) & -12.41 & -0.001 (0.003) & -0.52 \\
Farm population & -0.079 (0.015) & 67.55 & -0.081 (0.006) & 36.47 & -0.045 (0.007) & 38.95 \\
Education & 0.555 (0.111) & 76.14 & 0.638 (0.079) & 49.47 & 0.680 (0.089) & 69.98 \\
Density/1000 & 0.084 (0.027) & -0.20 & 0.134 (0.043) & -0.24 & -0.053 (0.049) & -0.05 \\
Wages/1000 & -0.062 (0.017) & -17.27 & -- & -- & -- & -- \\
$W$ Urban population & -0.034 (0.008) & -- & -0.031 (0.005) & -- & -- & -- \\
$W$ Farm population & 0.016 (0.019) & -- & 0.063 (0.007) & -- & -- & -- \\
$W$ Education & 0.051 (0.147) & -- & -0.438 (0.089) & -- & -- & -- \\
$W$ Density/1000 & -0.263 (0.038) & -- & -0.243 (0.091) & -- & -- & -- \\
$W$ Wages/1000 & 0.001 (0.037) & -- & -- & -- & -- & -- \\
\midrule
Sum of contributions & -- & 116.22 & -- & 73.29 & -- & 108.36 \\
Stability sum & 0.893 & -- & 0.946 & -- & 0.870 & -- \\
Observations & 9414 & -- & 18396 & -- & 18396 & -- \\
\bottomrule
\end{tabular}
\begin{minipage}{0.98\textwidth}
\scriptsize Notes: All QMLE specifications use two common factors. The baseline model sets the coefficient on $Wy_{t-1}$ and the spatial-error coefficient to zero. Standard errors use the martingale-sample covariance estimator with the pair-only HC1 correction in \eqref{eq:Omega_martingale_pair_hc1}. Contributions are percentages of the 1940--2000 increase in female labor-force participation.
\end{minipage}
\end{table}
Table~\ref{tab:application-female-lfp-baseline} shows that the contemporaneous spatial coefficient ranges from $0.479$ to $0.521$, while the own lag ranges from $0.364$ to $0.467$.
The resulting stability sums are $0.893$, $0.946$, and $0.870$, relatively large but below one.
The observed controls account for $116.2\%$, $73.3\%$, and $108.4\%$ of the historical increase across the three columns.
In the specification without wages or $WX$, the fall in farm population and increase in education produce the largest contributions.
The contribution sum of $108.4\%$ means that, holding the time effects and common factors fixed, the observed covariate changes imply a long-run increase in female labor-force participation of about $39.2$ percentage points, exceeding the observed increase of $36.2$ percentage points.
\subsection{Extension with a Lagged Spatial Outcome}
We next estimate \eqref{eq:application-model} with $\xi$ unrestricted, thereby
combining the delayed-neighbor channel of \textcite{fogli-veldkamp-2011} with
the contemporaneous outcome lag in \textcite{tziolas-elhorst-2023}.
Section~\ref{subsec:extension-lagged-spatial-outcome} extends the likelihood,
estimator, and asymptotic result, while Appendix~\ref{sec:extension-scores}
gives the corresponding spatial score decomposition. Thus, $Wy_{t-1}$ is not
treated as an ordinary external regressor.
With $\xi$ unrestricted, the stability sum becomes
$s:=\rho+\phi+\xi$. A spatially uniform component now evolves with dynamic
coefficient $(\phi+\xi)/(1-\rho)$, and its long-run multiplier has denominator
$1-s$. When $\rho<1$ and $\phi+\xi\geq0$, stability of this component is
equivalent to $s<1$; full stability continues to require the spectral-radius
condition in \eqref{eq:extension-spatiotemporal-stability}.
For comparison with the published decomposition, the contribution of control
$k$ is now
\begin{equation*}
100\times\frac{1}{N}\mathbf{1}_{N}'
\left[(1-\phi)I_N-(\rho+\xi)W\right]^{-1}
(\beta_k I_N+\gamma_k W)\mathbf{1}_{N}
\frac{\Delta x_k}{\Delta\mathrm{LFP}}.
\end{equation*}
For the specification without wages or spatially lagged controls, we also
report $L=0$ and $L=2$ estimates.
\begin{table}[htbp]
\centering
\caption{Lagged-spatial-outcome extension: published BC-QMLE and proposed QMLE}
\label{tab:application-female-lfp-extension}
\scriptsize
\setlength{\tabcolsep}{3pt}
\begin{tabular}{lrrrr}
\toprule
& \multicolumn{2}{c}{Estimate (standard error)} & \multicolumn{2}{c}{Contribution (\%)} \\
\cmidrule(lr){2-3}\cmidrule(lr){4-5}
Variable & Published BC-QMLE & Proposed QMLE & Published BC-QMLE & Proposed QMLE \\
\midrule
\multicolumn{5}{l}{\emph{Wages and spatially lagged controls included}} \\
$y_{i,t-1}$ & 0.153 (0.012) & 0.436 (0.051) & -- & -- \\
$Wy_t$ & -0.074 (0.075) & 0.530 (0.015) & -- & -- \\
$Wy_{t-1}$ & 0.504 (0.026) & -0.073 (0.043) & -- & -- \\
Urban population & 0.017 (0.004) & 0.009 (0.005) & -2.02 & -9.81 \\
Farm population & -0.081 (0.007) & -0.075 (0.013) & 54.19 & 55.49 \\
Education & 0.017 (0.069) & 0.592 (0.113) & -0.51 & 88.99 \\
Density/1000 & 0.031 (0.064) & 0.076 (0.026) & -0.18 & -0.18 \\
Wages/1000 & -0.053 (0.015) & -0.062 (0.017) & -8.09 & -10.28 \\
$W$ Urban population & -0.036 (0.009) & -0.032 (0.007) & -- & -- \\
$W$ Farm population & -0.120 (0.015) & 0.022 (0.020) & -- & -- \\
$W$ Education & -0.033 (0.127) & 0.119 (0.158) & -- & -- \\
$W$ Density/1000 & -0.664 (0.127) & -0.239 (0.039) & -- & -- \\
$W$ Wages/1000 & -0.059 (0.035) & 0.025 (0.043) & -- & -- \\
$W$ error & 0.595 (0.073) & -- & -- & -- \\
Sum of contributions & -- & -- & 43.40 & 124.21 \\
Stability sum & 0.582 & 0.893 & -- & -- \\
Observations & 9414 & 9414 & -- & -- \\
\addlinespace
\multicolumn{5}{l}{\emph{Wages excluded; spatially lagged controls included}} \\
$y_{i,t-1}$ & 0.148 (0.010) & 0.595 (0.206) & -- & -- \\
$Wy_t$ & -0.340 (0.089) & 0.509 (0.021) & -- & -- \\
$Wy_{t-1}$ & 0.381 (0.017) & -0.150 (0.164) & -- & -- \\
Urban population & 0.011 (0.003) & 0.017 (0.005) & 1.77 & -2.57 \\
Farm population & -0.070 (0.005) & -0.054 (0.099) & 21.27 & 5.35 \\
Education & 0.060 (0.058) & 0.656 (0.972) & 1.88 & 68.67 \\
Density/1000 & 0.093 (0.077) & 0.022 (0.452) & -0.10 & -0.09 \\
$W$ Urban population & 0.021 (0.007) & -0.019 (0.058) & -- & -- \\
$W$ Farm population & -0.083 (0.012) & 0.052 (0.016) & -- & -- \\
$W$ Education & 0.053 (0.100) & -0.416 (0.331) & -- & -- \\
$W$ Density/1000 & -0.820 (0.162) & -0.060 (0.587) & -- & -- \\
$W$ error & 0.830 (0.090) & -- & -- & -- \\
Sum of contributions & -- & -- & 24.81 & 71.36 \\
Stability sum & 0.189 & 0.953 & -- & -- \\
Observations & 18396 & 18396 & -- & -- \\
\addlinespace
\multicolumn{5}{l}{\emph{Wages and spatially lagged controls excluded}} \\
$y_{i,t-1}$ & 0.127 (0.009) & 0.389 (0.043) & -- & -- \\
$Wy_t$ & 0.087 (0.052) & 0.509 (0.011) & -- & -- \\
$Wy_{t-1}$ & 0.454 (0.017) & -0.025 (0.034) & -- & -- \\
Urban population & 0.006 (0.003) & -0.002 (0.003) & 0.89 & -0.69 \\
Farm population & -0.090 (0.005) & -0.042 (0.008) & 30.76 & 37.43 \\
Education & 0.004 (0.055) & 0.702 (0.098) & 0.17 & 73.99 \\
Density/1000 & 0.089 (0.075) & -0.051 (0.048) & 0.03 & -0.05 \\
$W$ error & 0.432 (0.051) & -- & -- & -- \\
Sum of contributions & -- & -- & 31.86 & 110.67 \\
Stability sum & 0.668 & 0.873 & -- & -- \\
Observations & 18396 & 18396 & -- & -- \\
\addlinespace
\bottomrule
\end{tabular}
\begin{minipage}{0.98\textwidth}
\scriptsize Notes: The published BC-QMLE columns report the Tziolas--Elhorst bias-corrected QML results transcribed from their replication output. The proposed QMLE specifications use two common factors, and their standard errors use the martingale-sample covariance estimator with the pair-only HC1 correction in \eqref{eq:Omega_martingale_pair_hc1}. Contributions use the 1940--2000 changes reported in their Table~1. Their specifications include a spatial autoregressive error coefficient; the proposed specifications do not.
\end{minipage}
\end{table}
\begin{table}[htbp]
\centering
\caption{Robustness to the SERL enrichment order with two common factors}
\label{tab:application-projection-robustness}
\small
\begin{tabular}{lrrr}
\toprule
Parameter & $L=0$ & $L=1$ & $L=2$ \\
\midrule
$Wy_t$ & 0.471 (0.011) & 0.509 (0.011) & 0.507 (0.011) \\
$y_{i,t-1}$ & 0.523 (0.021) & 0.389 (0.043) & 0.378 (0.038) \\
$Wy_{t-1}$ & -0.142 (0.016) & -0.025 (0.034) & -0.017 (0.030) \\
Urban population & 0.002 (0.003) & -0.002 (0.003) & -0.001 (0.003) \\
Farm population & -0.021 (0.006) & -0.042 (0.008) & -0.044 (0.008) \\
Education & 0.692 (0.099) & 0.702 (0.098) & 0.677 (0.097) \\
Density/1000 & -0.001 (0.024) & -0.051 (0.048) & -0.049 (0.047) \\
\midrule
Stability sum & 0.852 & 0.873 & 0.869 \\
\bottomrule
\end{tabular}
\begin{minipage}{0.90\textwidth}
\small Notes: Standard errors use the martingale-sample covariance estimator with the pair-only HC1 correction in \eqref{eq:Omega_martingale_pair_hc1}. For $L=0$, the smallest estimated eigenvalue of $\Sigma_\eta$ is effectively zero, so those standard errors should be interpreted cautiously because the interior covariance condition is nearly binding.
\end{minipage}
\end{table}
Table~\ref{tab:application-female-lfp-extension} compares the proposed QMLE
with the corresponding bias-corrected QML estimates of
\textcite{tziolas-elhorst-2023}.
The estimated coefficient $\xi$ is negative in all three specifications, at $-0.073$, $-0.150$, and $-0.025$, while $\widehat\rho$ ranges only from $0.509$ to $0.530$.
The raw $\xi$ is not, however, the reduced-form effect corresponding to the lagged-neighbor coefficient in \textcite{fogli-veldkamp-2011}.
The dynamic
transition satisfies
\begin{equation*}
B(\rho)^{-1}(\phi I_N+\xi W)
=
\phi I_N+(\xi+\rho\phi)W
+\rho(\xi+\rho\phi)W^2+\cdots.
\end{equation*}
The leading neighbor-lag effect $\xi+\rho\phi$ is positive in every
specification, taking values $0.158$, $0.152$, and $0.173$.
Table~\ref{tab:application-projection-robustness} varies the SERL enrichment order
in the specification without wages or $WX$. Moving from $L=1$ to $L=2$ changes
the allocation between $\phi$ and $\xi$, but leaves $\rho$, the stability sum,
and the leading neighbor-lag effect relatively similar. By contrast, $L=0$
materially changes the dynamic coefficients.
\subsection{Second-Order QML}
The first-order estimates deliberately omit a spatial-error filter. To assess
whether this omission drives the covariate decomposition, our preferred
specification is the unrestricted second-order QML in
\eqref{eq:extension-second-order-model}. In the application it takes the form
\begin{equation*}
\begin{split}
B_*(\kappa)y_t
={}&(\phi I_N+\xi_1W+\xi_2W^2)y_{t-1}\\
&+X_t\beta_0+WX_t\beta_1+W^2X_t\beta_2
+\mathbf{1}_{N}\delta_t+\Lambda f_t+\varepsilon_t.
\end{split}
\end{equation*}
For the first two empirical specifications, the control polynomial includes
$X_t$, $WX_t$, and $W^2X_t$. For the specification that originally excludes
$WX_t$, it includes $X_t$ and $WX_t$: the latter is required because filtering
the original $X_t\beta$ term by $B(\psi)$ generates a first spatial lag. The
model contains no separately labeled spatial-error coefficient. Instead, it
allows the second-order filter and the transformed lag and control coefficients
implied by \eqref{eq:extension-spatial-error-restrictions} without imposing
those restrictions.
Let
\begin{equation*}
Q_2
=(1-\phi)I_N-(\rho_1+\xi_1)W-(\rho_2+\xi_2)W^2.
\end{equation*}
The contribution of control $k$ is computed as
\begin{equation*}
100\times\frac{1}{N}\mathbf{1}_{N}'Q_2^{-1}
(\beta_{0k}I_N+\beta_{1k}W+\beta_{2k}W^2)\mathbf{1}_{N}
\frac{\Delta x_k}{\Delta\mathrm{LFP}},
\end{equation*}
with an omitted coefficient set to zero. We report the symmetric filter
coefficients $(\rho_1,\rho_2)$, not a structural assignment of the two roots,
and use the corrected covariance construction in
Appendix~\ref{sec:extension-scores}.
Because $W\mathbf{1}_{N}=\mathbf{1}_{N}$,
the reduced-form transition maps a spatially uniform change into another
uniform change.
We summarize its one-period persistence by the dynamic coefficient for a
uniform change,
\begin{equation*}
d_{\mathrm{unif}}
:=
\frac{\phi+\xi_1+\xi_2}{(1-\kappa_1)(1-\kappa_2)}
=
\frac{\phi+\xi_1+\xi_2}{1-\rho_1-\rho_2}.
\end{equation*}
A value below one means that a spatially uniform change contracts from one
period to the next;
stability across all eigenvalues of $W$ is governed by
\eqref{eq:extension-second-order-stability}.
\begin{table}[htbp]
\centering
\caption{Preferred second-order QMLE estimates}
\label{tab:application-second-order}
\scriptsize
\setlength{\tabcolsep}{2.5pt}
\resizebox{\textwidth}{!}{
\begin{tabular}{lrrrrrr}
\toprule
& \multicolumn{2}{c}{Wages and $WX$} & \multicolumn{2}{c}{No wages; with $WX$} & \multicolumn{2}{c}{No wages or $WX$} \\
\cmidrule(lr){2-3}\cmidrule(lr){4-5}\cmidrule(lr){6-7}
Variable & Estimate (SE) & Contribution (SE) & Estimate (SE) & Contribution (SE) & Estimate (SE) & Contribution (SE) \\
\midrule
$\rho_1$ & 0.203 (0.032) & -- & 0.199 (0.020) & -- & 0.219 (0.020) & -- \\
$\rho_2$ & 0.468 (0.038) & -- & 0.487 (0.025) & -- & 0.453 (0.027) & -- \\
$y_{i,t-1}$ & 0.614 (0.025) & -- & 0.552 (0.029) & -- & 0.568 (0.032) & -- \\
$Wy_{t-1}$ & 0.037 (0.033) & -- & 0.006 (0.026) & -- & 0.033 (0.025) & -- \\
$W^2y_{t-1}$ & -0.356 (0.043) & -- & -0.266 (0.038) & -- & -0.303 (0.035) & -- \\
Urban population & 0.013 (0.004) & -17.00 (9.03) & 0.019 (0.005) & -21.58 (12.86) & 0.014 (0.005) & -5.80 (7.14) \\
Farm population & -0.024 (0.009) & 37.61 (27.89) & -0.049 (0.007) & 55.67 (29.64) & -0.057 (0.008) & 19.97 (18.74) \\
Education & 0.527 (0.084) & 82.39 (38.63) & 0.646 (0.072) & 73.70 (37.17) & 0.536 (0.065) & 71.79 (26.21) \\
Density/1000 & 0.060 (0.015) & -0.21 (0.09) & 0.098 (0.040) & -0.46 (0.27) & 0.071 (0.039) & -0.23 (0.16) \\
Wages/1000 & -0.060 (0.017) & 3.92 (28.81) & -- & -- & -- & -- \\
$W$ Urban population & -0.011 (0.007) & -- & -0.003 (0.006) & -- & -0.018 (0.004) & -- \\
$W$ Farm population & 0.090 (0.015) & -- & 0.081 (0.011) & -- & 0.052 (0.007) & -- \\
$W$ Education & 0.608 (0.177) & -- & 0.081 (0.130) & -- & -0.373 (0.074) & -- \\
$W$ Density/1000 & -0.089 (0.044) & -- & -0.060 (0.067) & -- & -0.130 (0.070) & -- \\
$W$ Wages/1000 & 0.030 (0.040) & -- & -- & -- & -- & -- \\
$W^2$ Urban population & -0.015 (0.008) & -- & -0.026 (0.007) & -- & -- & -- \\
$W^2$ Farm population & -0.078 (0.018) & -- & -0.042 (0.013) & -- & -- & -- \\
$W^2$ Education & -0.935 (0.210) & -- & -0.607 (0.152) & -- & -- & -- \\
$W^2$ Density/1000 & -0.031 (0.042) & -- & -0.124 (0.074) & -- & -- & -- \\
$W^2$ Wages/1000 & 0.034 (0.048) & -- & -- & -- & -- & -- \\
\midrule
Sum of contributions & -- & 106.70 & -- & 107.33 & -- & 85.73 \\
Dynamic coefficient for uniform change & 0.901 & -- & 0.930 & -- & 0.908 & -- \\
Observations & 9414 & -- & 18396 & -- & 18396 & -- \\
\bottomrule
\end{tabular}
}
\begin{minipage}{0.98\textwidth}
\scriptsize Notes: The filter is $(I_N-\kappa_1W)(I_N-\kappa_2W)=I_N-\rho_1W-\rho_2W^2$. All specifications use two common factors and $L=1$. Standard errors use the martingale-sample covariance estimator with the pair-only HC1 correction in \eqref{eq:Omega_martingale_pair_hc1}; contribution standard errors use the delta method.
\end{minipage}
\end{table}
Table~\ref{tab:application-second-order} reports the preferred estimates. The
two symmetric filter coefficients are stable across specifications:
$\widehat\rho_1$ ranges from $0.199$ to $0.219$, and
$\widehat\rho_2$ from $0.453$ to $0.487$.
The dynamic coefficient for a spatially uniform change ranges from $0.901$ to $0.930$, below one in every case.
The coefficient on education remains positive and precisely estimated
after the second-order filter and the additional spatial transformations are
included.
Its estimates are $0.527$, $0.646$, and $0.536$, with standard errors $0.084$, $0.072$, and $0.065$.
The corresponding contributions are $82.4\%$, $73.7\%$, and $71.8\%$, with delta-method standard errors of $38.6$, $37.2$, and $26.2$ percentage points.
Educational gains provide the largest positive contribution in all three specifications.
The decline in the farm population also produces a positive estimated contribution, but the contribution is not statistically significant.
\subsection{Discussion and Comparison with Earlier Estimates}
\textcite{fogli-veldkamp-2011} omit the contemporaneous spatial lag, corresponding to $\rho=0$ in our notation.
Their OLS estimates of the own and neighbor lags are $0.664$ and $0.195$, respectively;
their preferred GMM estimates are $0.916$ and $0.570$, whose sum exceeds one.
In our baseline model, contemporaneous feedback induces a first-order coefficient
$\rho\phi$ on $Wy_{t-1}$, ranging from $0.184$ to $0.224$. In the extension, the corresponding coefficient
$\xi+\rho\phi$ ranges from $0.152$ to $0.173$.
Thus, our reduced-form neighbor effect is close to their OLS estimate of $0.195$ on $Wy_{t-1}$.
\textcite{tziolas-elhorst-2023}, by contrast, allocate most spatial propagation
to $Wy_{t-1}$. Across their specifications, its coefficient is $0.504$,
$0.381$, and $0.454$, while the contemporaneous coefficient is $-0.075$,
$-0.340$, and $0.087$. They simultaneously estimate spatial-error coefficients
of $0.595$, $0.830$, and $0.432$. This is a specification
difference. Our model sets the spatial-error coefficient to zero.
This difference also helps explain why our contemporaneous coefficient remains
near $0.5$. Because the same $W$ enters the outcome and error filters in
\textcite{tziolas-elhorst-2023}, the two components can be difficult to
distinguish with six periods. Denoting their outcome coefficient by $\rho$ and
their spatial-error coefficient by $\psi$, the two filters satisfy
\begin{equation}\label{eq:application-two-spatial-filters}
(I_N-\psi W)(I_N-\rho W)
=I_N-(\rho+\psi)W+\rho\psi W^2.
\end{equation}
Thus, a model with one first-order spatial filter may load part of an omitted
error filter onto its outcome coefficient. Consistent with this interpretation,
subtracting the sum of their estimated outcome and error coefficients from our
estimated spatial coefficient gives $0.009$, $0.018$, and $-0.011$ in the
three specifications. These are not parameter equivalences: the product also
induces the $W^2$ term in
\eqref{eq:application-two-spatial-filters} and spatially transforms the dynamic
and covariate components.
Moreover, a direct two-filter diagnostic reaches distinct local likelihood
optima from different starting values in every specification at both factor ranks.
We therefore do not treat the two roots as separately identified
outcome and error coefficients.
The central substantive difference concerns education, a covariate for which
one would expect a meaningful relationship with female labor-force
participation. The OLS estimate reported by
\textcite{fogli-veldkamp-2011} is $0.643$ with a standard error of $0.036$,
close to the range from $0.527$ to $0.646$ across our preferred second-order
specifications in Table~\ref{tab:application-second-order}. Their three-lag difference-GMM
estimate, $-0.975$ with a standard error of $0.587$, is instead imprecise,
indicating that their education estimate is sensitive to the estimator. The
corresponding estimates of \textcite{tziolas-elhorst-2023} are $0.017$, $0.060$,
and $0.004$, all statistically insignificant. Our preferred estimates assign education
contributions of $82.4\%$, $73.7\%$, and $71.8\%$, compared with $-0.5\%$,
$1.9\%$, and $0.2\%$ in their corresponding results reported in
Table~\ref{tab:application-female-lfp-extension}.
Some of education's association could be mediated through lagged outcomes or
absorbed by the spatial-error process. Nevertheless, the preferred estimates
give education a positive own coefficient and a substantial long-run
contribution after admitting the polynomial terms generated by such a filter.
The education discrepancy therefore is not explained solely by omission of the
spatial-error filter. It may also reflect differences in the treatment of common
shocks, initial conditions, and unit heterogeneity. By projecting county-specific
factor loadings on observed histories rather than estimating each loading
separately, our fixed-$T$ procedure may allocate persistent variation between
education and heterogeneous exposure to common shocks differently.
Taken together, our estimates are more consistent with the literature documenting an
important positive relationship between women's education and
labor-force participation; see, e.g., \textcite{goldin-2006}.
\section{Conclusion}\label{sec:conclusion}
This paper develops a likelihood-based framework for dynamic spatial panels
with common shocks in the fixed-$T$, large-$N$ regime.
The SERL specification and variance-components representation avoid estimating unit-specific
incidental parameters, while the conditional likelihood retains the spatial
Jacobian and the dynamic recursion without a parameter-dependent temporal
Jacobian. We also propose an efficient block-coordinate profile
algorithm, which exploits the structure of the likelihood to replace a
high-dimensional joint optimization with closed-form and GLS updates, low-rank
covariance updates, and a scalar search over the spatial coefficient. We
establish large-$N$ distribution theory and implement spatially corrected
sandwich inference using the score decomposition of
\textcite{li-yang-2021}. A lagged spatial outcome fits into the existing GLS
block, and a spatial-error filter preserves the profiling structure while
requiring a two-dimensional outer optimization. The baseline simulations show
small estimation error and generally near-nominal coverage.
In the empirical application, the lagged-neighbor effect is close to the original
OLS estimate of \textcite{fogli-veldkamp-2011}. Under the preferred second-order
specification, education remains positive and economically important after
allowing the transformations generated by a spatial-error filter; its
coefficient is also close to their positive OLS estimate.
These results accord with the literature documenting an important positive
relationship between women's education and labor-force participation.
\printbibliography[heading=bibintoc,title={References}]
\clearpage