Finite-Sample Properties of Model Specification Tests for Multivariate Dynamic Regression Models

Finite-Sample Properties of Model Specification Tests for Multivariate Dynamic Regression Models

Introduction

Economic phenomena typically involve multiple variables that are simultaneously determined and interact with one another. To analyze these interdependencies, economists commonly employ multiple-equation regression frameworks, such as the seemingly unrelated regression (SUR) model introduced by zellner1962eme. When error terms are correlated across equations, SUR and related multiple-equation frameworks provide a natural way to account for such dependence. For example, in systems of household demand equations, unobservable household-specific attributes may jointly influence the demand for multiple commodities, generating cross-equation dependence in the error terms (e.g., fiebig2001cte,henley1994tue). Likewise, when a cost equation is estimated jointly with cost share equations, unobservable firm-specific conditions may simultaneously influence both input allocation and total costs, inducing cross-equation dependence in the error terms (e.g., kumbhakar1996tct,brown1995sse). Ignoring such cross-equation dependence can lead to size distortions and misleading inference; see laitinen1978wdh,bewley1983trl,sarafidis2012csd. Consequently, model specification tests that explicitly account for cross-equation dependence in the error terms are essential for rigorous evaluation.

This consideration is particularly relevant in the evaluation of asset pricing models. Asset returns are jointly affected by common economic and financial shocks, generating substantial co-movement and dependence across assets (e.g., cieslak2021css,asgharian2023eus,bali2017ieu). Such cross-sectional dependence is often particularly pronounced during episodes of severe market stress, including financial crises, pandemics, and periods of heightened geopolitical risk (e.g., longin2001eci,bae2003nam,diebold2009mfa,zhang2020fmg,umar2022iru,jalloul2023emc). When the common dependence across asset returns is not fully captured by the asset pricing model under evaluation, the unexplained common component remains in the residuals, causing the resulting error terms to exhibit substantial cross-equation dependence. Accordingly, model specification tests that explicitly account for the cross-equation error dependence are essential for evaluating asset pricing models.

To this end, gibbons1989teg developed a multivariate test (the GRS test) to assess the validity of asset pricing models. The GRS test remains widely used in empirical studies; see, for example, fama1993crf,fama2015ffa,fama2016daf,fama2017itf,fama2018cf,fama2020ccs and cakici2013svm. However, it is derived under the assumption of independently and identically distributed multivariate normal errors. This assumption may be violated in practice, because asset returns often exhibit fat tails, skewness, and other departures from normality (e.g., cont2001epa,harvey2002csa,huang2012edr). Accordingly, even after controlling for common risk factors, the error terms in asset pricing models may still exhibit departures from normality. Indeed, affleck1989nta and zhou1993apt show that violations of normality can lead to substantial size distortions in the GRS test by inflating its Type I error rate.

In response to the restrictive distributional assumptions underlying the GRS test, the literature has developed heteroskedasticity- and autocorrelation-robust (HAR) tests based on nonparametric estimators of the long-run covariance matrix, e.g., newey1987sps,newey1994als,andrews1991hac. While attractive in empirical asset pricing because they do not require normality, conventional HAR tests can still exhibit substantial finite-sample size distortions, especially when serial dependence is strong; see haan1997pgr,muller2014hcs. These distortions can become even more severe in multivariate asset-pricing applications when many intercept restrictions are tested jointly (ray2008pha,ray2009tcr). This has led to a large literature on improved HAR inference, including fixed-$b$ methods of kiefer2000srt,kiefer2005nat, exponentiated-kernel or fixed-$p$ approaches of phillips2006sde,phillips2007lrv, loss-function-based bandwidth selection of sun2008obs, and more recent practical and optimality-oriented refinements of lazarus2018hir,lazarus2021spt.

Despite these improvements, HAR tests still suffer from a fundamental limitation. Because they are typically built on OLS estimators, their validity ultimately depends on exogeneity conditions strong enough to ensure consistency of the underlying OLS estimator. This issue is well known in the simultaneous-equations literature, where simultaneity induces correlation between regressors and error terms, rendering OLS estimator inconsistent (haavelmo1943sis,koopmans1945ses). More generally, in time-series environments, error terms may exhibit both serial dependence and cross-equation dependence (wagner2020fmo), while regressors may also be dynamically dependent on the error terms, thereby violating exogeneity conditions (engle1983exo,boswijk1997lamt). When such dynamic regressor--error dependence is present, the consistency of the underlying OLS estimator is no longer guaranteed, and HAR tests may no longer deliver valid inference.

To address these limitations, several recent studies have developed alternative estimation frameworks for regressions with dynamic regressor--error dependence, including perron2026fgls,baillie2024rit. However, the related literature relies on different notions of exogeneity, and the relationships among them have not been fully clarified. In particular, it remains unclear how these alternative conditions are related to one another, which of them are stronger than others, and how they determine the validity of competing estimators in multiple-equation settings. Moreover, these studies are confined to single-equation regressions and exclude an intercept term. Consequently, these frameworks do not directly apply to multivariate model specification tests that require joint restrictions across equations to be assessed while allowing for cross-equation error dependence. This issue arises in a wide range of settings, including, for example, asset pricing tests of whether pricing errors are jointly zero across assets.

This paper makes five contributions to the literature on multivariate model specification testing. First, we clarify the relationships among the alternative exogeneity conditions considered in the related literature. Second, within the unified two-step estimation and inference framework of murphy1985eit,pagan1986sre, we clarify the asymptotic properties of competing estimators in multiple-equation settings under dynamic regressor--error dependence. Third, we generalize durbin1970tsc's durbin1970tsc regression to a multiple-equation system with an intercept, allowing for both cross-equation error dependence and dynamic regressor--error dependence. Within this generalized framework, we show that the generalized Durbin estimator remains consistent under the weakest exogeneity condition considered and derive its asymptotic distribution. Fourth, building on this asymptotic theory, we develop a Wald test for model specification testing. Last, to improve finite-sample performance, we propose a bootstrap Wald test based on the generalized Durbin estimator.

The remainder of the paper is organized as follows. Section (ref) introduces a multiple-equation regression model with joint time-series dynamics for the regressors and error terms. It clarifies the relationships among the relevant exogeneity conditions, studies the inconsistency of OLS and alternative GLS-type estimators under dynamic regressor--error dependence, and establishes the asymptotic properties of the generalized Durbin estimator and develops a bootstrap-corrected Wald test for joint specification testing. Section (ref) presents Monte Carlo evidence on the finite-sample performance of the competing procedures. Section (ref) provides an empirical application. Section (ref) concludes.

Related literature

The problem considered in this paper, namely how to conduct model specification testing for a multiple-equation system under dynamic regressor--error dependence, does not fall within a single strand of literature but instead lies at the intersection of at least three closely related strands. First, in the presence of dynamic regressor--error dependence, the consistency of an estimator depends on the exogeneity condition imposed. We therefore begin by reviewing the literature on exogeneity conditions in time-series regressions. Second, differences in exogeneity conditions are closely linked to the validity regions of competing estimators. We then review the literature on estimation under dynamic regressor--error dependence. Third, these estimators often take the form of two-step procedures involving auxiliary parameter estimation. When estimation takes this form, valid statistical inference requires that first-step estimation error be properly taken into account. We therefore also review the literature on variance correction that properly accounts for first-step estimation error.

The consistency of an estimator depends fundamentally on the exogeneity condition assumed. In cross-sectional environments, it is often sufficient to focus on the contemporaneous relationship between regressors and error terms. In time-series environments, by contrast, one must also account for the dynamic dependence between regressors and error terms. Accordingly, the time-series regression literature has developed alternative exogeneity conditions based on orthogonality restrictions between regressors and error terms, such as strict exogeneity and present-and-past exogeneity, e.g. engle1983exo,stock2019ite,mayer2020cts,mikusheva2025lrw. However, because the error term in a time-series regression can itself be represented as a linear combination of current and past innovations, imposing conditions on the error term alone does not necessarily make clear to which shocks the regressors are related.

In light of this observation, perron2026fgls define exogeneity with respect to the underlying innovations rather than the error term itself. By contrast, baillie2024rit analyze the joint process of regressors and error terms through a Wold/VAR representation and formulate exogeneity as restrictions on their joint dynamics, in particular on lagged feedback and the contemporaneous innovation covariance structure. Although a variety of exogeneity conditions have been proposed to determine estimator consistency under dynamic regressor--error dependence, their relative strength and hierarchical relationships remain unclear. Clarifying these relationships is therefore essential for characterizing the validity regions of existing estimators.

These differences in exogeneity conditions are directly relevant for the validity of competing estimators. In time-series regressions, whether the OLS estimator is consistent depends on the exogeneity condition imposed. Violations of the exogeneity condition have traditionally motivated the use of alternative estimators, most notably instrumental variables (IV) and generalized method of moments (GMM), to achieve consistent estimation under endogeneity; see, for example, hayashi1983nee,hansen1982lsp,hansen1982ivp,cumby1983tst. More recently, the literature has also explored alternative approaches that seek to restore consistency by explicitly exploiting the dynamic dependence between regressors and error terms. perron2026fgls reformulate exogeneity in terms of the innovations underlying the error process and show that their proposed FGLS estimator (hereafter, the FD estimator) remains consistent even when OLS is inconsistent. baillie2024rit extend the analysis to the joint dynamics of regressors and error terms, and show that the durbin1970tsc's durbin1970tsc estimator remains consistent under a broader range of dependence structures than the FD estimator. Taken together, these studies suggest that, when endogeneity arises primarily from structured dynamic regressor--error dependence, consistent estimation may in some cases be recovered by exploiting that dependence directly, rather than by relying on IV- or GMM-based procedures. However, these studies are confined to single-equation regressions and do not accommodate an intercept term. Consequently, they do not directly address joint hypothesis testing in multiple-equation systems with cross-equation error dependence.

Estimators designed to handle dynamic regressor--error dependence often involve auxiliary parameter estimation. In such cases, the procedure can be viewed as a two-step procedure in which the first step estimates nuisance quantities and the second step uses them to construct the estimator of interest. In such two-step procedures, first-step estimation error can affect the asymptotic distribution of the second-step estimator at first order. Consequently, naive plug-in variance formulas that ignore first-step estimation error generally do not deliver valid inference (newey1984mmi,cattaneo2019tse). For example, the FD estimator of perron2026fgls relies on auxiliary parameter estimation, yet its asymptotic variance does not explicitly account for first-step estimation error. If that error contributes at first order, the resulting variance formula is misspecified, and inference based on it may therefore be invalid.\footnote{In the \href{https://at-noda.com/appendix/durbin_appendix.pdf}{Online Appendix} A.1, we investigate the finite-sample size consequences of ignoring first-step estimation error in inference under single-equation environments considered in the existing literature, and confirm that the resulting naive variance formulas can produce non-negligible over-rejection.} In this regard, murphy1985eit,pagan1986sre provide a unified treatment of variance correction that properly accounts for first-step estimation error. Building on this literature, we develop estimators for regressions with dynamic regressor--error dependence in a unified two-step estimation and inference framework.

Model and Test

Setup

We consider the following linear regressions:

where $y_{i,t}$ is the scalar dependent variable, $\ve{x}_{i,t}=(x_{i,t,1},x_{i,t,2},\ldots,x_{i,t,k_i})^\prime$ is a $k_i\times 1$ vector of regressors, and $u_{i,t}$ is a scalar disturbance. The number of regressors $k_i$ is allowed to vary across equations.

Let $\ve{y}_t=(y_{1,t},\ldots,y_{N,t})^\prime$, $\ve{u}_t=(u_{1,t},\ldots,u_{N,t})^\prime$ with ${\mathbb{E}}[\ve{u}_t]=\ve{0}$, and define $\ve{\alpha}=(\alpha_{1},\ldots,\alpha_{N})^\prime$ and $\ve{\beta}=(\ve{\beta}_{1}^\prime,\ldots,\ve{\beta}_{N}^\prime)^\prime$, where $k:=\sum_{i=1}^N k_i$. Further, define the block-diagonal regressor matrix \[ \ve{X}_t:=\operatorname{blkdiag}(\ve{x}_{1,t},\ldots,\ve{x}_{N,t}) =

\in\mathbb{R}^{k\times N}. \] Stacking the $N$ equations, we can write the system compactly as \[

=

+

+

, \] equivalently

To describe the joint second-order dynamics of regressors and disturbances, let $\ve{x}_t=(\ve{x}_{1,t}^\prime,\ldots,\ve{x}_{N,t}^\prime)^\prime$ with ${\mathbb{E}}[\ve{x}_t]=\ve{\mu}_{x}$, and define \[ \ve{z}_t=

\in \mathbb{R}^{m}, \qquad m:= k + N, \] with mean $\ve{\mu}_z:=[\ve{\mu}_x^\prime,\ \ve{0}^\prime]^\prime = {\mathbb{E}}[\ve{z}_t]$. Let $\bar{\ve{z}}_t:=\ve{z}_t-\ve{\mu}_z$.

We work in the Hilbert space $L^2(\Omega,\mathscr{F},\mathbb{P};\mathbb{R}^{m})$ equipped with inner product $\langle U,V\rangle={\mathbb{E}}[U^\prime V]$. To define the one-step-ahead linear innovation of $\{\bar{\ve z}_t\}$, define \[ \mathscr{H}_{t-1} := \overline{ \mathopen{}\originalleft\{ \sum_{j=1}^{p} \ve{C}_j \bar{\ve z}_{t-j} \;:\; p<\infty,\; \ve{C}_j\in\mathbb{R}^{m\times m} \mathclose{}\originalright\} }, \] that is, the closed linear subspace generated by all finite matrix-coefficient linear combinations of past values of $\{\bar{\ve z}_t\}$. Let $\mathscr{P}_{t-1}:L^2(\Omega,\mathscr{F},\mathbb{P};\mathbb{R}^m)\to \mathscr{H}_{t-1}$ denote the $L^2$-orthogonal projection onto $\mathscr{H}_{t-1}$. The one-step-ahead linear innovation is then defined by \[ \ve{\varepsilon}_t := \bar{\ve z}_t-\mathscr{P}_{t-1}[\bar{\ve z}_t] =

. \] By construction, $\ve{\varepsilon}_t$ is orthogonal to every element of $\mathscr{H}_{t-1}$ and therefore represents the component of $\bar{\ve z}_t$ that cannot be linearly predicted from the past of the system. Our definition formalizes the idea in sims1980mr that an innovation is a new component of the system, namely, the part that is not predicted from its past.

In general, the closed linear subspace $\mathscr{H}_{t-1}$ does not exclude infinite-order linear combinations of past values of $\{\bar{\ve z}_t\}$. However, in a finite sample, an unrestricted infinite-order linear predictor is not directly operational. In our joint time-series representation, the innovation covariance matrix governs contemporaneous dependence between regressor-side and disturbance-side innovations, whereas lagged dependence in the system governs dynamic regressor--disturbance dependence. By treating $\ve{\varepsilon}_{x,t}$ and $\ve{\varepsilon}_{u,t}$ as distinct contemporaneous shock blocks, we separate instantaneous shock dependence from lagged transmission mechanisms. This interpretation is consistent with the literature in which distinct components of a joint system are modeled as separate shock blocks; see, for example, cushman1997imp,zha1999brs. The following assumption therefore provides a transparent benchmark that makes the dependence structure of the joint system easier to interpret and analyze.

Under Assumption (ref)-[1], the centered joint process $\{\bar{\ve z}_t\}$ is covariance-stationary and purely nondeterministic. Hence, by the Wold decomposition, $\bar{\ve z}_t$ admits a unique orthogonal VMA($\infty$) representation with no deterministic component. Moreover, finite-predictor exactness at lag $p_0$ implies that $\bar{\ve z}_t$ admits a VAR($p_0$) innovations representation and can also be interpreted as practical approximations to the Wold representation for a broad class of stationary processes (lutkepohl2005nim). Under Assumption (ref)-[2], the innovation covariance matrix is block diagonal: \[ \ve{\Sigma} = {\mathbb{E}}[\ve{\varepsilon}_t\ve{\varepsilon}_t^\prime] =

,\quad \ve{\Sigma}_{xx}\in\mathbb{R}^{k\times k},\quad \ve{\Sigma}_{uu}\in\mathbb{R}^{N\times N}. \] The following proposition formalizes these implications and establishes the stability of the resulting VAR($p_0$) representation.

Proposition (ref) provides the fundamental probabilistic structure used in the subsequent analysis. Part (i) shows that the centered joint process $\{\bar{\ve z}_t\}$ admits an orthogonal Wold representation, while part (ii) implies that finite-predictor exactness yields a finite-order VAR innovations representation. Part (iii) further establishes stability and absolute summability of the associated coefficients. Finally, part (iv) delivers the covariance-stationarity, ergodicity, and moment conditions needed for the asymptotic arguments below. Accordingly, the joint system can be studied through a stable VAR representation in which contemporaneous dependence is encoded in $\ve{\Sigma}$ and dynamic dependence is encoded in the lag structure.

We adopt a reduced-form VAR($p_0$) as the benchmark because it summarizes multivariate dynamics with minimal structural assumptions and serves as a standard workhorse in empirical macroeconomics and finance (sims1980mr,stock2001var). In asset-pricing applications, it is common for all equations to share the same regressor vector. Leading examples include the CAPM of sharpe1964cap and lintner1965vra, the arbitrage pricing theory (APT) of ross1976atc, and the multi-factor models of fama1993crf,fama2015ffa,fama2016daf. In these cases, $\ve{x}_{i,t}= \ve{x}_t$ for all $i$. This common-regressor case is nested in our general setup, and all results under Assumption (ref) remain valid after a notational simplification.

Hereafter, we establish conditions for the consistency and asymptotic normality of estimators of $(\ve{\alpha}^\prime, \ve{\beta}^\prime)^\prime$ across a variety of data-generating processes (DGPs), characterized by restrictions on $\{\ve{\Psi}_j\}_{j=1}^{p_0}$ and $\ve{\Sigma}$.

Relationships among exogeneity conditions under dynamic regressor--disturbance dependence

In linear regression analysis, exogeneity assumptions play a central role in determining how the stochastic relationship between regressors and error terms affects the consistency and efficiency of estimators, as well as the validity of associated test statistics. In cross-sectional settings, a contemporaneous orthogonality condition, ${\mathbb{E}}[\ve{x}_t \ve{u}_t^\prime] = \ve{0}$, is often sufficient to characterize the consistency of the OLS estimator. In time-series environments, by contrast, contemporaneous orthogonality alone is generally insufficient to describe the forms of regressor--error dependence that are relevant for estimation and inference.

To address this issue, the literature has developed a range of exogeneity conditions. A traditional approach, following stock2019ite, characterizes dynamic exogeneity through conditional mean restrictions linking the regressors and the error process across time. In particular, strict exogeneity requires

so that the error term at time $t$ is mean independent of the entire past, present, and future path of the regressors. A weaker notion is present-and-past exogeneity, which requires only

that is, the disturbance is mean independent of the contemporaneous and lagged regressors, but may still be related to future regressors. It follows immediately that strict exogeneity implies present-and-past exogeneity.

However, the error-based exogeneity conditions of stock2019ite are not sufficiently fine-grained to distinguish orthogonality with the error process itself from orthogonality with the innovations underlying it. To make this distinction explicit, perron2026fgls formulate exogeneity conditions in terms of the Wold innovations of the error process. In their single-equation framework, the regressors are assumed to satisfy the contemporaneous orthogonality condition ${\mathbb{E}}[\ve{x}_t \ve{\varepsilon}_{u,t}^\prime]=\ve{0}$ as a maintained assumption, and dynamic exogeneity is then distinguished according to whether $\ve{x}_t$ is orthogonal to future or past innovations of the error process. We extend these definitions to the multiple-equation setting by using the innovation of the error process $\{\ve{\varepsilon}_{u,t}\}$ from the joint innovation process in Proposition (ref). The regressors are pre-determined if

that is, $\ve{x}_t$ is uncorrelated with all future innovations of error terms. The regressors are exogenous if

so that $\ve{x}_t$ is uncorrelated with all past innovations of error terms.

While the innovation-based conditions of perron2026fgls are well suited to clarifying the conditions under which OLS and GLS are consistent, they do not by themselves fully capture the dynamic relationship between $\ve{x}_t$ and $\ve{u}_t$ in a joint regressor--error system. To make these forms of dynamic regressor--error dependence explicit, baillie2024rit introduce a set of exogeneity conditions as restrictions on a joint VAR representation for $\bar{\ve{z}}_t=(\bar{\ve{x}}_t^\prime,\ve{u}_t^\prime)^\prime$.

Under the block-diagonal ($BD$) condition, the joint dynamics of regressors and error terms are block-diagonal:

Hence, neither lagged feedback from past error terms to current regressors nor lagged feedback from past regressors to current error terms is allowed.

The GLS-exogeneity ($GEXOG$) condition relaxes the $BD$ condition by allowing feedback from past error terms to current regressors, while continuing to exclude feedback from past regressors to current error terms:

Finally, under the error-block-diagonal ($EBD$) condition, both directions of dynamic interaction are permitted through the VAR coefficients, while the innovation covariance matrix remains block-diagonal:

Thus, the $BD$--$GEXOG$--$EBD$ classification makes explicit which forms of lagged regressor--error dependence are ruled out and which are allowed in the joint system, with $BD \subsetneq GEXOG \subsetneq EBD$.

Because different strands of the literature adopt different notions of exogeneity, the relationships among these conditions and their relative strength are not always transparent. Under Assumption (ref), however, the joint process $((\ve{x}_t-\ve{\mu}_x)^\prime,\ve{u}_t^\prime)^\prime$ admits a joint VAR/VMA representation, which provides a convenient framework for comparing these exogeneity conditions through the VAR coefficients $\{\ve{\Psi}_j\}$ and the innovation covariance matrix $\ve{\Sigma}$. In particular, the $BD$, $GEXOG$, and $EBD$ conditions may be viewed as structural restrictions on $(\ve{\Psi}_j,\ve{\Sigma})$ that characterize the admissible forms of dynamic feedback between regressors and error terms. The following proposition makes these relationships explicit under Assumption (ref).

Proposition (ref) provides a unified characterization of the exogeneity concepts developed in different strands of the literature and clarifies how these conditions are related to one another. In particular, it identifies the hierarchy among the VAR-based classes and shows how the innovation-based exogeneity conditions are positioned relative to that hierarchy.

Proposition (ref) and Remarks (ref) and (ref) imply the relationships among the exogeneity conditions shown in the following figure.

While the exogenous condition and the present-and-past exogeneity condition remain useful as auxiliary concepts for relating our framework to alternative strands of the literature, they do not directly parameterize the lagged feedback structure that determines estimator consistency in the multivariate system considered below. By contrast, $BD$, $GEXOG$, and $EBD$ are defined as block restrictions on the joint VAR representation and therefore explicitly distinguish the direction and extent of dynamic regressor--disturbance dependence. For this reason, we use the $BD$--$GEXOG$--$EBD$ hierarchy as the organizing framework for assessing the validity of the multiple-equation estimators in the subsequent analysis.

The inconsistency of the OLS estimator in multi-equation models

In empirical applications, the ordinary least squares (OLS) estimator is routinely used to estimate the parameter vector in (ref). For notational convenience, rewrite (ref) as

where $\ve{Z}_t^\prime=[\ve{I}_N,\ve{X}_t^\prime]\in\mathbb{R}^{N\times(N+k)}$ and $\ve{\kappa}_0=(\ve{\alpha}_0^\prime,\ve{\beta}_0^\prime)^\prime\in\mathbb{R}^{N+k}$. The OLS estimator based on (ref) is defined by

Although the OLS estimator in (ref) is widely used in time-series regressions, its consistency requires not only a full-rank condition on the regressors but also an appropriate orthogonality restriction. We begin with a standard identification condition.

Assumption (ref) is a standard full-rank condition on the regressors. It rules out exact multicollinearity and guarantees that the population normal equations admit a unique solution.

baillie2024rit impose the $BD$ condition as a sufficient condition for OLS consistency. By contrast, perron2026fgls use innovation-based exogeneity conditions, defined from the marginal Wold decomposition of $\{\ve{u}_t\}$, in their consistency analysis. As shown in Proposition (ref)(i), under Assumption (ref), the $BD$ condition is equivalent to strict exogeneity in our framework. These conditions, however, are stronger than necessary. Relatedly, stock2019ite emphasize present-and-past exogeneity as a key condition for OLS consistency. Covariance-based present-and-past exogeneity implies the contemporaneous orthogonality condition ${\mathbb{E}}[\ve{X}_t\ve{u}_t]=\ve{0}$, but is not itself necessary. Thus, the exogeneity conditions used in these studies are nested within a broader class of data-generating processes under which OLS remains consistent. The next proposition characterizes the necessary and sufficient moment condition for OLS consistency in our setting and clarifies how it relates to the exogeneity classes discussed above.

Proposition (ref) shows that OLS consistency is equivalent to ${\mathbb{E}}[\ve Z_t\ve u_t]=\ve 0$, or, since ${\mathbb{E}}[\ve u_t]=\ve 0$, equivalently to the contemporaneous orthogonality condition ${\mathbb{E}}[\ve X_t\ve u_t]=\ve 0$. The key issue, therefore, is how this orthogonality can arise in a dynamic environment.

Under the joint VAR/VMA representation implied by Assumption (ref), blockwise orthogonality of contemporaneous innovations does not in general guarantee the contemporaneous orthogonality condition ${\mathbb{E}}[\ve X_t\ve u_t]=\ve 0$. Dynamic cross-dependence through the off-diagonal blocks of $\{\ve\Psi_j\}$, or equivalently through the impulse responses $\{\ve\Xi_i\}$, allows past innovations from one block to affect the current value of the other. As a result, ${\mathbb{E}}[\bar{\ve x}_t\ve u_t^\prime]$ can be nonzero, and hence ${\mathbb{E}}[\ve X_t\ve u_t]$ need not vanish.

The next proposition characterizes this moment condition within the joint VAR-based exogeneity classes. In particular, $BD$ is sufficient for contemporaneous orthogonality, whereas the condition need not hold under $GEXOG$. Hence, it need not hold more generally over the broader $EBD$ class either. This clarifies the precise sense in which OLS can be inconsistent in dynamic multi-equation environments.

Proposition (ref)(ii)--(iii) do not rule out the contemporaneous orthogonality condition for particular DGPs in $GEXOG$ or $EBD$; they only show that these classes do not guarantee it uniformly.

Two-Step Estimation and Asymptotic Variance Corrections

The previous subsection established that the exact condition for OLS consistency is the contemporaneous moment restriction ${\mathbb{E}}[\ve X_t\ve u_t]=\ve 0$. While the $BD$ condition guarantees this restriction, the broader exogeneity classes $GEXOG$ and $EBD$ do not do so in general. Hence, in time-series environments with regressor--error dependence, OLS-based estimation and inference are not generally reliable.

This observation motivates the use of estimators that remain consistent even outside $BD$ by explicitly accounting for dynamic regressor--error dependence. However, in multiple-equation systems, parameter consistency alone is not sufficient for valid specification testing. Because disturbances may remain cross-sectionally correlated across equations, valid and efficient inference must also account for the cross-equation covariance structure.

For this reason, we consider GLS-type estimators that account for dynamic regressor--error dependence while also exploiting the cross-equation covariance structure. Because the relevant error dynamics and covariance matrices are unknown, such estimators are naturally implemented in two steps. In particular, a first step estimates the nuisance parameters governing the feasible transformation, and a second step applies feasible GLS using those estimates.

This two-step setup introduces an additional layer of uncertainty for inference. The first-step estimation error propagates into the second-step estimator, so asymptotically valid inference requires a variance formula that explicitly incorporates this additional source of uncertainty. Accordingly, this subsection develops the asymptotic distribution of such GLS-type estimators within the two-step inference framework of murphy1985eit and pagan1986sre.

Throughout this subsection, let $\ve\psi_{jt}(\cdot)$ ($j=1,2$) denote moment functions based on the observed data $\{(\ve y_t,\ve X_t)\}$. Let $\ve\vartheta_0$ be the nuisance parameter identified in the first step and $\ve\kappa_0$ the parameter of interest identified in the second step. The first-step estimator $\widehat{\ve\vartheta}$ is defined as a solution to \[ \bar{\ve\psi}_1(\ve\vartheta) := \frac{1}{T}\sum_{t=1}^T \ve\psi_{1t}(\ve\vartheta) = \ve 0, \] whereas the second-step estimator $\widehat{\ve\kappa}$ is defined as a solution to \[ \bar{\ve\psi}_2(\ve\kappa,\widehat{\ve\vartheta}) := \frac{1}{T}\sum_{t=1}^T \ve\psi_{2t}(\ve\kappa,\widehat{\ve\vartheta}) = \ve 0. \]

Let $(\ve{\kappa}_0,\ve{\vartheta}_0)$ denote the true parameter values. Define the population Jacobian blocks, evaluated at $(\ve{\kappa}_0,\ve{\vartheta}_0)$, by \[ \ve{A}_{11} := {\mathbb{E}}\Bigg[ \frac{\partial \ve{\psi}_{1t}(\ve{\vartheta}_0)}{\partial \ve{\vartheta}^\prime} \Bigg], \qquad \ve{A}_{22} := {\mathbb{E}}\Bigg[ \frac{\partial \ve{\psi}_{2t}(\ve{\kappa}_0,\ve{\vartheta}_0)}{\partial \ve{\kappa}^\prime} \Bigg], \qquad \ve{A}_{21} := {\mathbb{E}}\Bigg[ \frac{\partial \ve{\psi}_{2t}(\ve{\kappa}_0,\ve{\vartheta}_0)}{\partial \ve{\vartheta}^\prime} \Bigg]. \] Since the first-step moment function $\ve{\psi}_{1t}(\ve{\vartheta})$ does not depend on $\ve{\kappa}$, there is no corresponding $\ve{A}_{12}$ block.

The two-step estimators admit the linear representations (ref)--(ref), where $\widetilde{\ve{\psi}}_{2t}$ is the adjusted second-step moment obtained by substituting the first-step linearization into the first-order expansion of the second-step moment around $(\ve\kappa_0,\ve\vartheta_0)$:

where

Let $\ve\psi_{1,t}:=\ve\psi_{1t}(\ve\vartheta_0)$ and $\ve\psi_{2,t}:=\ve\psi_{2t}(\ve\kappa_0,\ve\vartheta_0)$. Define the long-run covariance (LRV) blocks of the two moment processes by

Then

The naive plug-in covariance formula that treats $\widehat{\ve{\vartheta}}$ as known ignores the effect of first-step estimation error. In particular, it uses only the second-step moment $\ve{\psi}_{2t}$, together with its long-run covariance block $\ve S_{22}$, and omits the correction terms involving $\ve A_{21}\ve A_{11}^{-1}$ as well as the cross-covariance blocks $\ve S_{12}$ and $\ve S_{21}$. Remark (ref) shows that the second-step estimator is driven by the adjusted moment $\widetilde{\ve{\psi}}_{2t}$ rather than the original moment $\ve{\psi}_{2t}$. Accordingly, the relevant long-run covariance matrix is $\ve{S}_{\widetilde{2}\widetilde{2}}$, which generally differs from $\ve S_{22}$. The following remark is closely related to the special cases discussed in pagan1986sre.

The inconsistency of OLS-based GLS estimators in multiple-equation models

A common response to serial correlation in $\ve{u}_t$ is to adopt feasible GLS-type procedures that model and remove the error dynamics. However, although such procedures can address serial correlation under appropriate orthogonality conditions, they do not by themselves resolve the endogeneity problem highlighted above when ${\mathbb{E}}[\ve{X}_t\ve{u}_t]\neq\ve{0}$. In particular, when the first-step OLS estimator is inconsistent, the resulting residuals do not consistently recover the true disturbance process, so subsequent GLS transformations based on those residuals need not restore consistency. Nevertheless, such procedures remain a standard workhorse in empirical time-series regressions.

In practice, these methods are typically implemented in two steps. First, one estimates (ref) by OLS and obtains residuals $\{\widehat{\ve{u}}_t^{OLS}\}$. Second, one fits a finite-order model for the residual dynamics (for example, equation-by-equation AR($p$) models, or a vector autoregression for the residual process), applies the cochrane1949als's cochrane1949als quasi-differencing transformation based on the estimated lag polynomial, and then computes a feasible GLS estimator. We refer to such estimators as CO-type estimators. In the multiple-equation setting, nagakura2024cot develop a multivariate CO-type estimator in which the disturbance vector is modeled as a VAR($p$) process.

Suppose first that the disturbance admits the VAR($p_0$) representation

where the characteristic polynomial $\ve A_{uu}(z):=\ve I_N-\sum_{j=1}^{p_0}\ve\Psi_{uu,j}z^j$ has no zeros in or on the unit disk. Define the transformed variables $\ve{y}_{CO,t}:=\ve A_{uu}(L)\ve y_t$ and $\ve{Z}_{CO,t}^\prime:=\ve A_{uu}(L)\ve Z_t^\prime$. If $\ve A_{uu}(L)$ and $\ve\Sigma_{uu}$ were known, the corresponding infeasible GLS estimator would be

Assumption (ref) is the GLS analogue of the full-rank condition in Assumption (ref). It rules out degeneracy of the transformed regressors and ensures that the population GLS normal equations admit a unique solution.

Proposition (ref) shows that the FCO estimator is guaranteed to be consistent under $BD$, where the first-step OLS estimator consistently recovers the disturbance process and hence yields a valid feasible Cochrane--Orcutt transformation. Outside the $BD$ class, however, consistency is no longer guaranteed.

Since the FCO estimator is obtained in two steps, valid inference requires an asymptotic variance formula that accounts for first-step estimation error. Because consistency of the FCO estimator is guaranteed only under $BD$, we confine attention to the asymptotic variance under $BD$ condition.

Proposition (ref) shows that, under $BD$, the FCO estimator is asymptotically equivalent to the infeasible GLS estimator for the correctly specified transformed system. Consequently, no first-step asymptotic variance correction is required. In addition, under $BD$, the FCO estimator attains the GLS efficiency bound for that transformed system.

It follows immediately from Proposition (ref) that we can conduct model specification tests for linear restrictions on $\ve{\kappa}_0$ under the $BD$ condition.

The inconsistency of FD estimator in multiple-equation models

The previous subsection showed that the OLS-based CO procedure is guaranteed to be consistent only under $BD$. This naturally raises the question of whether one can construct a feasible GLS estimator that remains consistent over a broader class of dynamic regressor--disturbance dependence. A leading attempt in this direction is the FD procedure proposed by perron2026fgls for single-equation time-series regressions.

In the multiple-equation setting, the key maintained restriction underlying FD is that the disturbance depends only on its own lags. In the joint VAR for $\bar{\ve z}_t=((\ve x_t-\ve\mu_x)^\prime,\ve u_t^\prime)^\prime$, this amounts to imposing $\ve\Psi_{ux,j}=\ve 0$ for all $j$, so that lagged regressors are excluded from the $u$-equation, while lagged disturbances may still enter the regressor block through $\ve\Psi_{xu,j}\neq\ve 0$. Hence, in our notation, the maintained dynamic model of FD corresponds to the upper block-triangular region that we call $GEXOG$, together with innovation orthogonality $\ve\Sigma_{xu}=\ve 0$. Accordingly, we formulate the multivariate analogue of FD under $GEXOG$ and then characterize its consistency region.

Imposing these restrictions on $\ve{\Psi}_j$ yields

where the dynamic feedback from $\ve{x}_t$ to $\ve{u}_t$ is ruled out by imposing $\ve{\Psi}_{ux,j}=\ve{0}$ for all $j$, while feedback from past disturbances to the regressors is allowed through $\ve{\Psi}_{xu,j}$.

Substituting (ref) into (ref) and rearranging terms yields

Using $\ve X_t^\prime\ve\beta=\ve B\ve x_t$ with $\ve B:=\mathrm{blkdiag}(\ve\beta_1^\prime,\ldots,\ve\beta_N^\prime)\in\mathbb R^{N\times k}$, the term $-\sum_{j=1}^{p_0}\ve\Psi_{uu,j}\ve X_{t-j}^\prime\ve\beta$ in (ref) can be rewritten as $-\sum_{j=1}^{p_0}\ve\Psi_{uu,j}\ve X_{t-j}^\prime\ve\beta =-\sum_{j=1}^{p_0}\ve\Psi_{uu,j}\ve B\,\ve x_{t-j} =\sum_{j=1}^{p_0}\ve\Delta_j\,\ve x_{t-j}$ where, at the population level, $\ve\Delta_j := -\ve\Psi_{uu,j}\ve B$, $j=1,\ldots,p_0$.

For estimation, however, we treat $\{\ve\Delta_j\}_{j=1}^{p_0}$ as unrestricted nuisance parameters, which keeps the augmented regression linear in the unknown coefficients. Under $GEXOG$, the true parameter values satisfy $\ve\Delta_j=-\ve\Psi_{uu,j}\ve B$ for each $j$, so this reparameterization does not change the population model. Rather, it embeds the $GEXOG$ restriction in a larger linear regression system, thereby avoiding the bilinear cross-parameter restriction in estimation.

With this linear reparameterization, (ref) becomes a linear regression model in which $\ve c=(\ve{I}_N-\sum_{j=1}^{p_0}\ve{\Psi}_{uu,j})\ve{\alpha}$, $\ve\beta$, $\{\ve\Psi_{uu,j}\}_{j=1}^{p_0}$, and $\{\ve\Delta_j\}_{j=1}^{p_0}$ enter linearly. Since no cross-equation parameter restrictions are imposed in this augmented representation, the first step of the FD procedure can be implemented equation by equation with OLS. Specifically, for each $i=1,\ldots,N$, we regress $y_{i,t}$ on the contemporaneous regressors in equation $i$, the lagged dependent vectors $\{\ve y_{t-j}\}_{j=1}^{p_0}$, and the corresponding lagged regressor terms. The coefficients on $\{\ve y_{t-j}\}_{j=1}^{p_0}$ across equations then deliver the row blocks of $\widehat{\ve\Psi}_{uu,1},\ldots,\widehat{\ve\Psi}_{uu,p_0}$, while those on the lagged regressor terms deliver $\widehat{\ve\Delta}_1,\ldots,\widehat{\ve\Delta}_{p_0}$. The resulting residuals, denoted by $\widehat{\ve\varepsilon}_{u,t}$, are used to estimate $\ve\Sigma_{uu}$. This leads to the following two-step FD procedure.

Lemma (ref) establishes the basic stochastic properties of the augmented regression needed for equation-by-equation OLS estimation. We now impose a full-rank condition on the augmented regressors to obtain consistency of the first-step estimator.

In practice, the true lag order $p_0$ is unknown and must be estimated from the data. We can use the Bayesian information criterion (BIC) of schwarz1978edm, which consistently selects the true lag order under standard regularity conditions.

Let $\ve Z_{\mathrm{FD},t}$ denote the population filtered regressor matrix obtained by replacing $\widehat{\ve\Psi}_{uu,j}$ with $\ve\Psi_{uu,j}$ in Step 2. To ensure identification in the second-step filtered GLS regression, we impose the following full-rank condition.

Assumption (ref) is the filtered-system analogue of the full-rank condition in Assumption (ref). It rules out degeneracy of the filtered regressors and ensures that the population GLS normal equations admit a unique solution.

Proposition (ref) shows that FD has a larger consistency region than FCO estimators. While CO-type procedures are guaranteed to be consistent only under $BD$, FD remains consistent under the broader $GEXOG$ condition, and under $BD$ it is asymptotically equivalent to FCO. This reflects the fact that, under $GEXOG$, the disturbance process $\{\ve u_t\}$ admits a closed VAR($p_0$) representation, so the quasi-differencing step is asymptotically correctly specified.

Although the multivariate FD estimator is consistent under $GEXOG$, it is inherently a two-step estimator. The filtered regressors used in the second-step GLS regression, as well as the feasible weighting matrix, are constructed from first-step estimates of the disturbance dynamics and innovation covariance matrix. Consequently, the second-step estimating equation inherits first-step estimation error, and the naive plug-in covariance formula that treats these generated objects as fixed is generally asymptotically invalid.

To obtain valid inference, we analyze $\widehat{\ve\kappa}^{\mathrm{FD}}$ within the general two-step framework developed in Section (ref). Let $\ve\vartheta^{\mathrm{FD}} := \bigl(\operatorname*{vec}(\ve\Psi_{uu,1})^\prime,\ldots,\operatorname*{vec}(\ve\Psi_{uu,p_0})^\prime\bigr)^\prime$ collect the nuisance parameters governing the disturbance dynamics. The second-step moment is given by the GLS score associated with the filtered regression generated by $\ve\vartheta^{\mathrm{FD}}$, while the first-step augmented OLS estimator admits an asymptotically linear representation. As a result, the limit distribution of $\widehat{\ve\kappa}^{\mathrm{FD}}$ is governed not by the original second-step moment alone, but by the adjusted moment \[ \widetilde{\ve\psi}_{2t}^{\mathrm{FD}} := \ve\psi_{2t}^{\mathrm{FD}} - \ve A_{21}^{\mathrm{FD}} (\ve A_{11}^{\mathrm{FD}})^{-1} \ve\psi_{1t}^{\mathrm{FD}}, \] where the Jacobian blocks are \[ \ve A_{11}^{\mathrm{FD}} :={\mathbb{E}}\Big[\frac{\partial \ve\psi_{1t}^{\mathrm{FD}}(\ve\vartheta_0^{\mathrm{FD}})} {\partial \ve\vartheta^{\mathrm{FD}\prime}}\Big],\quad \ve A_{22}^{\mathrm{FD}} :={\mathbb{E}}\Big[\frac{\partial \ve\psi_{2t}^{\mathrm{FD}}(\ve\kappa_0,\ve\vartheta_0^{\mathrm{FD}})} {\partial \ve\kappa^\prime}\Big],\quad \ve A_{21}^{\mathrm{FD}} :={\mathbb{E}}\Big[\frac{\partial \ve\psi_{2t}^{\mathrm{FD}}(\ve\kappa_0,\ve\vartheta_0^{\mathrm{FD}})} {\partial \ve\vartheta^{\mathrm{FD}\prime}}\Big]. \] Hence the asymptotic variance takes the form $(\ve A_{22}^{\mathrm{FD}})^{-1}\ve S_{\widetilde2\widetilde2}^{\mathrm{FD}}(\ve A_{22}^{\mathrm{FD}})^{-1\prime}$, where $\ve S_{\widetilde2\widetilde2}^{\mathrm{FD}} := \sum_{h=-\infty}^{\infty}{\mathbb{C}\rm{ov}}\bigl(\widetilde{\ve\psi}_{2,t}^{\mathrm{FD}},\widetilde{\ve\psi}_{2,t-h}^{\mathrm{FD}}\bigr)$ is the long-run covariance matrix of the adjusted moment.

Explicit expressions for the Jacobian and covariance blocks entering $\ve V^{\mathrm{FD}}$ are reported in the \href{https://at-noda.com/appendix/durbin_appendix.pdf}{Online Appendix} A.2.1.

The FD estimator is a two-step estimator, since both the quasi-differenced regressors and the GLS weight are constructed from first-step estimates of the disturbance dynamics. Accordingly, whenever $\ve A_{21}^{\mathrm{FD}}\neq \ve 0$, first-step estimation error enters the first-order asymptotic variance, so a variance formula that treats these objects as known, as in perron2026fgls, is not valid in general. In that case, valid inference requires the Murphy--Topel correction in Theorem (ref). By contrast, if $\ve A_{21}^{\mathrm{FD}}=\ve 0$ as in Remark (ref)(i), the first-step effect is asymptotically irrelevant at the $T_{\mathrm{eff}}^{-1/2}$ order, the correction vanishes, and $\widehat{\ve\kappa}^{\mathrm{FD}}$ attains the infeasible GLS bound.

Theorem (ref) immediately yields a Wald test for linear restrictions on $\ve{\kappa}_0$ based on the FD estimator.

Proposition (ref) clarifies the scope of the FD estimator. Although FD can attain the infeasible GLS bound when the first-step effect is asymptotically negligible, its validity still hinges on correct specification of a restricted VAR for $\{\ve u_t\}$. This requirement fails once the DGP belongs to $EBD\setminus GEXOG$: in that case, $\{\ve u_t\}$ generally does not admit a closed finite-order VAR driven solely by its own lags, so the quasi-differencing filter underlying FD becomes misspecified. Hence, like OLS and FCO, FD is not robust to bidirectional dynamic feedback between regressors and disturbances. This limitation motivates the generalized Durbin-type estimators developed in the next subsection, which remain consistent under the weaker $EBD$ condition.

Asymptotic properties of generalized Durbin estimator in multi-equation models

We extend Durbin regression to a multiple-equation framework and derive an estimator that remains consistent even under the $EBD$ condition. Under Assumption (ref), the centered joint process $\bar{\ve z}_t$ admits the blockwise VAR($p_0$) representation in Proposition (ref), which can be written as

Substituting ((ref)) for $\ve{u}_t$ into ((ref)) yields the following regression model:

where $\ve{\gamma}=\bigl(\ve{I}_N- \sum_{j=1}^{p_0}\ve{\Psi}_{uu,j}\bigr)\ve{\alpha} - \sum_{j=1}^{p_0}\ve{\Psi}_{ux,j}\ve{\mu}_{x}$, $\ve{\Lambda}_j=\ve{\Psi}_{ux,j} -\ve{\Psi}_{uu,j}\ve{B}$ and $\ve{B}=\mathrm{blkdiag}\big(\ve{\beta}_1^\prime,\ve{\beta}_2^\prime,\ldots,\ve{\beta}_N^\prime \big)$.

For estimation, we treat $\{\ve\Lambda_j\}_{j=1}^{p_0}$ as unrestricted nuisance parameters so that the augmented regression remains linear in the unknown coefficients. This reparameterization does not alter the population regression implied by (ref); it simply avoids imposing the bilinear cross-parameter restriction $\ve\Lambda_j=\ve\Psi_{ux,j}-\ve\Psi_{uu,j}\ve B$ during estimation.

With this linear reparameterization, (ref) becomes a linear regression model in which $\ve\gamma$, $\ve\beta$, $\{\ve\Psi_{uu,j}\}_{j=1}^{p_0}$, and $\{\ve\Lambda_j\}_{j=1}^{p_0}$ enter linearly. Since no cross-equation parameter restrictions are imposed in this augmented representation, the first step of the generalized Durbin procedure can be implemented equation by equation by OLS. In particular, the coefficients on $\{\ve y_{t-j}\}_{j=1}^{p_0}$ deliver the row blocks of $\widehat{\ve\Psi}_{uu,1},\ldots,\widehat{\ve\Psi}_{uu,p_0}$, while those on $\{\ve x_{t-j}\}_{j=1}^{p_0}$ deliver $\widehat{\ve\Lambda}_1,\ldots,\widehat{\ve\Lambda}_{p_0}$. These first-step estimates are then used to recover $\{\ve\Psi_{ux,j}\}_{j=1}^{p_0}$, construct the generalized Durbin transformation, and estimate the innovation covariance matrix. This leads to the following two-step generalized Durbin procedure.

Lemmas (ref) and (ref) justify the first step of this procedure. In particular, for fixed lag order $p_0$, the equation-by-equation OLS estimators from the augmented regressions are consistent for the corresponding population coefficients. Hence the first-step estimates consistently recover the nuisance parameters that enter the generalized Durbin transformation. In practice, the lag order is unknown. Lemma (ref) shows that the Bayesian information criterion consistently selects the true lag order under the stated regularity conditions. Therefore, the feasible procedure based on the selected lag order remains asymptotically valid.

Let $\ve{Z}_{\mathrm{GD},t}$ denote the population transformed regressor matrix obtained by replacing the first-step estimates in Step 2 with their population counterparts. To ensure identification in the second-step GLS regression, we impose the following full-rank condition.

Theorem (ref) shows that the generalized Durbin estimator remains consistent on the broader $EBD$ region. However, it is still inherently a two-step estimator: both the transformed regressors used in the second-step GLS regression and the feasible weighting matrix are constructed from first-step estimates of the disturbance dynamics and innovation covariance matrix. As a result, the second-step estimating equation inherits first-step estimation error, so a naive plug-in covariance formula that treats these generated objects as fixed is generally asymptotically invalid.

To obtain valid inference, we analyze $\widehat{\ve\kappa}^{\mathrm{GD}}$ within the general two-step framework developed in Section (ref). In the generalized Durbin procedure, the first step estimates the primitive nuisance parameter $\ve\vartheta^{\mathrm{GD}}:=\bigl(\ve\theta_1^\prime,\ldots,\ve\theta_N^\prime,\ve\mu_x^\prime\bigr)^\prime$, where each $\ve\theta_i$ collects the coefficients in the $i$th augmented Durbin regression. The second-step transformed regression, however, depends on a different set of nuisance objects, \[ \ve\eta^{\mathrm{GD}} := \bigl( \operatorname*{vec}(\ve\Psi_{uu,1})^\prime,\ldots,\operatorname*{vec}(\ve\Psi_{uu,p_0})^\prime,\ \operatorname*{vec}(\ve\Psi_{ux,1})^\prime,\ldots,\operatorname*{vec}(\ve\Psi_{ux,p_0})^\prime,\ \ve\mu_x^\prime \bigr)^\prime =: h(\ve\vartheta^{\mathrm{GD}}), \] which enter the generalized Durbin transformation directly. In particular, the matrices $\{\ve\Psi_{ux,j}\}$ are recovered from the first-step estimates through the identity \[ \ve\Psi_{ux,j}=\ve\Lambda_j+\ve\Psi_{uu,j}\ve B, \qquad \ve B=\mathrm{blkdiag}(\ve\beta_1^\prime,\ldots,\ve\beta_N^\prime), \qquad j=1,\ldots,p_0. \]

Accordingly, the second-step estimating equation inherits first-step estimation error through the mapping $h(\cdot)$. Let $\ve\psi_{2t}^{\mathrm{GD}}(\ve\kappa,\ve\eta^{\mathrm{GD}})$ denote the GLS score associated with the generalized Durbin transformed regression. Since the first-step estimator of $\ve\vartheta^{\mathrm{GD}}$ admits an asymptotically linear representation, the limit distribution of $\widehat{\ve\kappa}^{\mathrm{GD}}$ is governed not by the original second-step moment alone, but by the adjusted moment \[ \widetilde{\ve\psi}_{2t}^{\mathrm{GD}} := \ve\psi_{2t}^{\mathrm{GD}} - \ve A_{21}^{\mathrm{GD}} (\ve A_{11}^{\mathrm{GD}})^{-1} \ve\psi_{1t}^{\mathrm{GD}}, \] where the Jacobian blocks are \[ \ve A_{11}^{\mathrm{GD}} := {\mathbb{E}}\mathopen{}\originalleft[\frac{\partial \ve\psi_{1t}^{\mathrm{GD}}(\ve\vartheta_0^{\mathrm{GD}})}{\partial \ve\vartheta^{\mathrm{GD}\prime}}\mathclose{}\originalright], \ \ve A_{22}^{\mathrm{GD}} :={\mathbb{E}}\Big[\frac{\partial \ve\psi_{2t}^{\mathrm{GD}}(\ve\kappa_0,\ve\vartheta_0^{\mathrm{GD}})} {\partial \ve\kappa^\prime}\Big],\ \ve A_{21}^{\mathrm{GD}} := {\mathbb{E}}\Big[\frac{\partial \ve\psi_{2t}^{\mathrm{GD}}(\ve\kappa_0,\ve\eta_0^{\mathrm{GD}})} {\partial \ve\eta^{\mathrm{GD}\prime}} \frac{\partial \ve\eta^{\mathrm{GD}}}{\partial \ve\vartheta^{\mathrm{GD}\prime}}\Big]. \] Hence the asymptotic variance takes the form $(\ve A_{22}^{\mathrm{GD}})^{-1}\ve S_{\widetilde2\widetilde2}^{\mathrm{GD}}(\ve A_{22}^{\mathrm{GD}})^{-1\prime}$, where $\ve S_{\widetilde2\widetilde2}^{\mathrm{GD}} := \sum_{h=-\infty}^{\infty}{\mathbb{C}\rm{ov}}\bigl(\widetilde{\ve\psi}_{2,t}^{\mathrm{GD}},\widetilde{\ve\psi}_{2,t-h}^{\mathrm{GD}}\bigr)$ is the long-run covariance matrix of the adjusted moment.

Explicit expressions for the Jacobian and covariance blocks entering $\ve V^{\mathrm{FD}}$ are reported in the \href{https://at-noda.com/appendix/durbin_appendix.pdf}{Online Appendix} A.2.2.

The generalized Durbin estimator is also a two-step estimator, since both the transformed regressors and the GLS weight used in the second step are constructed from first-step estimates of the nuisance parameters. In particular, the second-step transformation depends on the effective nuisance parameter $\ve\eta^{\mathrm{GD}}=h(\ve\vartheta^{\mathrm{GD}})$, which is itself generated from the primitive first-step estimator $\widehat{\ve\vartheta}^{\mathrm{GD}}$. Accordingly, whenever $\ve A_{21}^{\mathrm{GD}}\neq \ve 0$, first-step estimation error enters the first-order asymptotic variance, so a covariance formula that treats the transformed nuisance objects as known is not valid in general. In that case, valid inference requires the Murphy--Topel/Pagan correction in Theorem (ref). By contrast, if $\ve A_{21}^{\mathrm{GD}}=\ve 0$, the first-step effect is asymptotically irrelevant at the $T_{\mathrm{eff}}^{-1/2}$ order, the correction disappears, and $\widehat{\ve\kappa}^{\mathrm{GD}}$ attains the infeasible GLS bound.

As a direct consequence of Theorem (ref), we can conduct model specification tests for linear restrictions on $\ve{\kappa}_0$ under any of the exogeneity conditions $BD$, $GEXOG$, or $EBD$.

Bootstrap-based Wald test using the generalized Durbin estimator

Theorem (ref) establishes the $\sqrt{T_{\mathrm{eff}}}$-asymptotic normality of the generalized Durbin estimator $\widehat{\ve\kappa}^{\mathrm{GD}}$, which in turn justifies Wald inference as shown in Corollary (ref). In finite samples, however, it is well known that Wald tests based on first-order asymptotic approximations may exhibit non-negligible size distortions; see, for example, hall1996bcv,hwang2023fsc.

To mitigate these finite-sample distortions, we consider a bootstrap-based Wald test. Under suitable regularity conditions, bootstrap methods often improve finite-sample inference for asymptotically pivotal statistics by providing a more accurate approximation to their null distributions (hall1996bcv,davidson1999sdb,mackinnon2006bme). In a time-series setting, however, an ordinary i.i.d.\ bootstrap is inappropriate, because resampling the observations as if they were independent fails to preserve the dynamic dependence structure in the data. A bootstrap procedure for dependent data is therefore required. In our setting, bootstrap pseudo-samples must also be generated so as to be consistent with the null-imposed bootstrap data-generating process while preserving the relevant dynamic dependence. For this purpose, a sieve bootstrap is well suited, because it generates bootstrap pseudo-samples recursively from an estimated dynamic model. Accordingly, following buhlmann1997sbt, we approximate the relevant dependence structure by a VAR($p$), where $p$ may increase with $T$.

The bootstrap procedure is implemented as follows. We first impose the null hypothesis on the unrestricted GD estimator to obtain a restricted estimator and the associated null-imposed residuals. We then fit a VAR($\hat p$) approximation, generate bootstrap pseudo-samples by resampling the estimated innovations, and re-estimate the full feasible GD procedure in each bootstrap replication. The detailed algorithm is reported in the following table.

This bootstrap procedure is designed to capture the dynamic dependence structure relevant for the generalized Durbin setting. In particular, it is implemented under the null hypothesis, preserves the dependence structure of the null-imposed joint process, and incorporates the estimation uncertainty arising from the feasible generalized Durbin procedure. The bootstrap $p$-value is then obtained from the empirical distribution of the bootstrap Wald statistics.

The formal justification of the proposed bootstrap procedure is deferred to the Appendix. In particular, the \href{https://at-noda.com/appendix/durbin_appendix.pdf}{Online Appendix} A.3 establishes the consistency of the null-imposed bootstrap construction, the validity of the fitted VAR approximation, the conditional first-order expansion of the bootstrap feasible GD estimator, and the resulting asymptotic validity of the bootstrap GD-based Wald test.

Simulation Experiments

In this section, we examine the finite-sample properties of the estimators for multiple-equation systems and their associated test statistics under the $BD$, $GEXOG$, and $EBD$ conditions through Monte Carlo experiments. First, we assess the accuracy of the estimators under each exogeneity condition. Second, we investigate the size and size-adjusted power of the corresponding tests.

Simulation Design

In our simulation, we focus on the case in which the regressors are common across all equations, as discussed in Remark (ref). This specification is particularly relevant for asset pricing models, where each portfolio equation is typically assumed to depend on the same set of risk factors. Prominent examples include the multifactor models of fama1993crf's (fama1993crf,fama2015ffa) and the arbitrage pricing theory of ross1976atc's ross1976atc. Motivated by this setting, we consider a multiple-equation regression model in which the regressors are common across all equations:

where $\ve{y}_t = (y_{1,t},\ldots,y_{N,t})^\prime \in \mathbb{R}^{N}$, $\ve{x}_t$ is a $k \times 1$ vector of common regressors ($k=k_i$ for all $i$), and $\ve{\beta}=(\ve{\beta}_1^\prime,\ldots,\ve{\beta}_N^\prime)^\prime\in\mathbb{R}^{kN}$. Let $\ve{X}_t^\prime = \ve{I}_N \otimes \ve{x}_t^\prime \in \mathbb{R}^{N \times kN}$ and $\ve{Z}_t^\prime := [\ve{I}_N,\ve{X}_t^\prime] \in \mathbb{R}^{N \times (k+1)N}$, so that $\ve{\kappa} = (\ve{\alpha}^\prime,\ve{\beta}^\prime)^\prime \in \mathbb{R}^{(k+1)N}$ collects the intercept and slope parameters. As noted in Remark (ref), the specification in (ref) is a special case of (ref).

For the centered process $\bar{\ve{z}}_t=(\ve{x}_t^\prime-\ve{\mu}_x^\prime,\ve{u}_t^\prime)^\prime$, we consider the following three VAR(1) data-generating processes under the $BD$, $GEXOG$, and $EBD$ conditions, respectively.

In the simulation design of baillie2024rit, lagged dependent variables are explicitly included in the regression. In contrast, our framework does not require specifying a finite lag structure for $\ve{y}_t$ a priori in the estimating equation. Instead, under the maintained joint dynamic specification for $\ve{z}_t=(\ve{x}_t^\prime-\ve{\mu}_x^\prime,\ve{u}_t^\prime)^\prime$, together with the structural relation for $\ve{y}_t$, one obtains an equivalent Durbin-type representation. Hence, the dynamics of $\ve{y}_t$ are captured implicitly through the joint dynamics of $(\ve{x}_t,\ve{u}_t)$ rather than being specified directly in advance.

In all cases, we consider sample sizes $T \in \{100, 200, 400, 800\}$. To examine the effects of the number of regressors and the number of equations on the finite-sample accuracy of estimators and the performance of tests, we consider $k \in \{2, 4\}$ and $N \in \{5, 10\}$. In each replication, we generate $T+500$ observations and discard the first 500 observations as burn-in. For each DGP, we conduct 10,000 Monte Carlo replications.

For the VAR(1) coefficient matrix, the block restrictions depend on the DGP ($BD$/ $GEXOG$/ $EBD$) as in (ref)--(ref). For the diagonal blocks, we control persistence through the spectral norm. Specifically, we draw random orthogonal matrices $\ve Q_x \in \mathbb R^{k\times k}$ and $\ve Q_u \in \mathbb R^{N\times N}$ via QR decompositions of i.i.d.\ Gaussian matrices and set \[ \ve{\Psi}_{xx,1}=c_{xx}\,\ve Q_x,\qquad \ve{\Psi}_{uu,1}=c_{uu}\,\ve Q_u, \] so that $\|\ve{\Psi}_{xx,1}\|_2=c_{xx}$ and $\|\ve{\Psi}_{uu,1}\|_2=c_{uu}$ before any global rescaling for stationarity. For the cross blocks, we set $r=\min(k,N)$ and independently draw matrices with orthonormal columns, $\ve A_1,\ve A_2\in\mathbb R^{k\times r}$ and $\ve B_1,\ve B_2\in\mathbb R^{N\times r}$, such that \[ \ve A_1^\prime \ve A_1=\ve A_2^\prime \ve A_2=\ve I_r, \qquad \ve B_1^\prime \ve B_1=\ve B_2^\prime \ve B_2=\ve I_r. \] We then define \[ \ve{\Psi}_{xu,1}=\frac{c_{xu}}{\sqrt r}\,\ve A_1\ve B_1^\prime,\qquad \ve{\Psi}_{ux,1}=\frac{c_{ux}}{\sqrt r}\,\ve B_2\ve A_2^\prime. \] This scaling keeps the Frobenius norms of the cross blocks comparable across $r$. The baseline parameter values are set as $c_{xx}=0.4$, $c_{xu}=0.7$, $c_{ux}=0.3$, and $c_{uu}=0.5$. The full coefficient matrix is uniformly rescaled whenever necessary so that its spectral radius does not exceed $0.91$. Specifically, if the spectral radius of the initially generated VAR coefficient matrix, $\rho(\ve{\Psi})$, exceeds $0.91$, then all blocks of $\ve{\Psi}$ are multiplied by the common factor $0.91/\rho(\ve{\Psi})$.

For the innovation covariance matrices, we generate $\ve\Sigma_{xx}$ and $\ve\Sigma_{uu}$ separately as random positive definite matrices using a Cholesky-type construction. Specifically, for each block we draw a lower-triangular matrix $\ve L$ with unit diagonal and i.i.d.\ strictly lower-triangular entries from $\mathrm{Unif}(-\delta,\delta)$, and then set the covariance matrix equal to $\ve L\ve L^\prime$. In the implementation, we use $\delta=0.1$ for $\ve\Sigma_{xx}$ and $\delta=0.5$ for $\ve\Sigma_{uu}$. The two innovation blocks are generated independently, so that the contemporaneous innovation covariance between $\ve x_t$ and $\ve u_t$ is zero.

Given the coefficient matrix and the innovation covariance, we generate the centered vector $\ve{z}_t:=(\ve{x}_t^\prime-\ve{\mu}_x^\prime,\ve{u}_t^\prime)^\prime$ according to the VAR(1) recursion \[ \ve{z}_t=\ve{\Psi}\ve{z}_{t-1}+\ve{\varepsilon}_t, \qquad \ve{\varepsilon}_t\sim \mathcal N(\ve{0},\ve{\Sigma}), \quad \ve{\Sigma}=

. \] After generating $\ve x_t-\ve\mu_x$, we add back the mean vector $\ve\mu_x$ to obtain $\ve x_t$. Finally, the dependent variables are generated from the structural equation \[ \ve{y}_t=\ve{\alpha}+\big(\ve{I}_N \otimes \ve{x}_t^\prime\big)\ve{\beta}+\ve{u}_t,\quad t=1,\ldots,T. \] For the regressor mean, we set $\ve{\mu}_x=0.3\,\ve{1}_k$. All slope coefficients are set to one, i.e.\ $\ve{\beta}=\ve{1}_{Nk}$.

Estimation Accuracy

Following baillie2024rit, we examine the finite-sample accuracy of four estimators for multiple-equation systems: OLS, FCO, FD, and GD. Specifically, we compute and compare their bias and mean squared error (MSE). Because the parameters are vector-valued, we summarize bias and MSE using scalar measures based on the Euclidean norm for ease of interpretation.

Let $\ve\kappa_0=(\ve{\alpha}_0^\prime,\ve{\beta}_0^\prime)^\prime\in\mathbb R^{(k+1)N}$ denote the true parameter vector, and let $\widehat{\ve\kappa}_T$ be an estimator of $\ve\kappa_0$ based on a sample size $T$. As a scalar measure of bias, we report the Euclidean norm of the bias vector: \[ \mathrm{bias}(\widehat{\ve\kappa}_T) := \|{\mathbb{E}}[\widehat{\ve\kappa}_T]-\ve\kappa_0\|_2. \] We further define the scalar measure of MSE by \[ \mathrm{MSE}(\widehat{\ve\kappa}_T) :={\mathbb{E}}\mathopen{}\originalleft[\|\widehat{\ve\kappa}_T-\ve\kappa_0\|_2^2\mathclose{}\originalright] ={\mathbb{E}}\mathopen{}\originalleft[(\widehat{\ve\kappa}_T-\ve\kappa_0)^\prime(\widehat{\ve\kappa}_T-\ve\kappa_0)\mathclose{}\originalright]. \] Then the bias--variance decomposition is given by \[ \mathrm{MSE}(\widehat{\ve\kappa}_T) = \mathrm{bias}(\widehat{\ve\kappa}_T)^2 + \mathrm{tr}\mathopen{}\originalleft({\mathbb{V}\rm{ar}}(\widehat{\ve\kappa}_T)\mathclose{}\originalright), \] where ${\mathbb{V}\rm{ar}}(\widehat{\ve\kappa}_T)={\mathbb{E}}\mathopen{}\originalleft[(\widehat{\ve\kappa}_T-{\mathbb{E}}[\widehat{\ve\kappa}_T])(\widehat{\ve\kappa}_T-{\mathbb{E}}[\widehat{\ve\kappa}_T])^\prime\mathclose{}\originalright]$. This decomposition allows us to assess finite-sample estimation accuracy in terms of the relative contributions of squared bias and variance.

In Monte Carlo experiments with $M$ replications, let $\widehat{\ve\kappa}_T^{(m)}$ denote the estimator obtained in replication $m$. We define the Monte Carlo average by $\bar{\ve\kappa}_T := 1/M\sum_{m=1}^M \widehat{\ve\kappa}_T^{(m)}$ and the corresponding Monte Carlo covariance matrix by $\widehat{{\mathbb{V}\rm{ar}}}(\widehat{\ve\kappa}_T):=1/M\sum_{m=1}^M\bigl(\widehat{\ve\kappa}_T^{(m)}-\bar{\ve\kappa}_T\bigr)\bigl(\widehat{\ve\kappa}_T^{(m)}-\bar{\ve\kappa}_T\bigr)^\prime$. Then the Monte Carlo analogue of the scalar measure of bias and MSE are given by \[ \widehat{\mathrm{bias}}(\widehat{\ve\kappa}_T) := \mathopen{}\originalleft\|\bar{\ve\kappa}_T-\ve\kappa_0\mathclose{}\originalright\|_2, \quad \widehat{\mathrm{MSE}}(\widehat{\ve\kappa}_T) := \widehat{\mathrm{bias}}(\widehat{\ve\kappa}_T)^2 + \mathrm{tr}\bigl( \widehat{{\mathbb{V}\rm{ar}}}(\widehat{\ve\kappa}_T) \bigr). \]

In what follows, we compare the accuracy of the estimators for multiple-equation systems across the DGPs in (ref)--(ref). We set the true parameters to $\ve{\alpha}_0=\ve{0}_N$ and $\ve{\beta}_0=\ve{1}_{Nk}$, so that $\ve{\kappa}_0=(\ve{0}_N^\prime,\ve{1}_{Nk}^\prime)^\prime$.

Table (ref) reports the bias and MSE of the GD, FD, FCO, and OLS estimators under $BD$. For all estimators, the bias is small throughout, and the MSE declines as the sample size $T$ increases, indicating improved finite-sample accuracy in larger samples. Overall, the bias also tends to become smaller as $T$ grows. These findings are in line with the asymptotic theory, under which all four estimators are consistent under $BD$. The MSEs of FD and FCO are nearly identical across all designs and are uniformly smaller than those of GD and OLS. This pattern is consistent with Proposition (ref), which establishes the asymptotic equivalence of FD and FCO under $BD$. It is also in line with Proposition (ref), which shows that FCO attains the GLS efficiency bound under $BD$.

Table (ref) reports the bias and MSE of the GD, FD, FCO, and OLS estimators under $GEXOG$. The bias of the OLS estimator remains non-negligible across all combinations of $N$ and $k$, even as the sample size $T$ increases. This pattern is in line with Proposition (ref), which establishes consistency of OLS under $BD$, while the more general $GEXOG$ condition does not in itself guarantee OLS consistency. Although Proposition (ref) implies that FCO is not generally consistent under $GEXOG$, its finite-sample bias and MSE remain much smaller than those of OLS and do not deteriorate substantially in the reported designs. This finding is consistent with the broader literature showing that misspecified GLS procedures can still perform reasonably well in finite samples. For example, baillie2024rit report under their $GEXOG$ design that traditional FGLS can substantially outperform OLS and, in some small-sample parameterizations, even yield a smaller MSE than FD. Similarly, perron2026fgls show that even an incorrect GLS correction can achieve a lower MSE than OLS, and koreisha2001gls provide both theoretical efficiency bounds and simulation evidence indicating that incorrect GLS procedures can be quite robust in finite samples.

Table (ref) reports the bias and MSE of the GD, FD, FCO, and OLS estimators under $EBD$. The bias of GD tends to decline as the sample size $T$ increases across all combinations of $N$ and $k$, in line with Theorem (ref), which establishes consistency of GD under $EBD$. By contrast, FD, FCO, and OLS generally exhibit larger bias than GD and do not display the same clear tendency toward zero, consistent with the fact that these estimators are not generally consistent under $EBD$. GD also attains the smallest MSE in all reported cases.

Taken together, the results in Tables (ref)--(ref) show that the relative accuracy of the estimators depends crucially on the underlying exogeneity condition. Under $BD$, both FD and FCO perform well, while under $GEXOG$, FD continues to perform well and FCO can still perform reasonably well in finite samples despite the lack of general consistency. However, once the exogeneity conditions required for their consistency are violated, the finite-sample accuracy of FD and FCO deteriorates. By contrast, the GD estimator remains stable in terms of bias and MSE across all exogeneity conditions considered here. These findings are broadly consistent with the theoretical results established in the previous section.

Size and size-adjusted power of tests

In this subsection, we examine the size and size-adjusted power of test statistics constructed from the GD, FD, FCO, and OLS estimators. The size of a test is its probability of rejecting the null hypothesis when the null is true, and thus indicates whether the test controls the Type I error at the nominal significance level in finite samples. A test with severe size distortions may exhibit substantial over-rejection or under-rejection. By contrast, size-adjusted power measures the ability of a test to detect false null hypotheses after adjusting for differences in finite-sample size. Considering both measures allows for a fair comparison of the competing tests.

In linear asset-pricing models, the intercept vector $\ve{\alpha}$ captures pricing errors, because each $\alpha_i$ represents the component of the mean excess return on asset $i$ that is not explained by its factor exposure. Accordingly, the validity of a linear asset-pricing model is typically assessed through the joint null hypothesis $H_0:\ve{\alpha}=\ve{0}$. In empirical applications, this null hypothesis is often tested using the gibbons1989teg's gibbons1989teg test (the GRS test) or the HAR test of lazarus2018hir,lazarus2021spt.

\paragraph{GRS test} The GRS test of gibbons1989teg is one of the most widely used test statistics for assessing the validity of asset-pricing models; see, for example, fama1993crf,fama2015ffa,fama2016daf,fama2017itf,fama2018cf,fama2020ccs and cakici2013svm. Here, we suppose that the regressor vector is common across equations, i.e., $\ve{x}_{i,t}:=\ve{x}_t$ for all $i=1,\ldots,N$ as in Remark (ref). Define \[ \bar{\ve{x}}=\frac{1}{T}\sum_{t=1}^T\ve{x}_t,\quad \ve{S}_x=\frac{1}{T-1}\sum_{t=1}^T(\ve{x}_t-\bar{\ve{x}})(\ve{x}_t-\bar{\ve{x}})^\prime,\quad \widehat{\ve{\Sigma}}_{u}^{\mathrm{OLS}}=\frac{1}{T-k-1}\sum_{t=1}^T\widehat{\ve{u}}_t^{\mathrm{OLS}}\widehat{\ve{u}}_t^{\mathrm{OLS}\prime}, \] where $\widehat{\ve{u}}_t^{\mathrm{OLS}}$ denotes the OLS residual vector. Then the GRS statistic is given by \[ \mathrm{GRS} = \frac{T(T-N-k)}{N(T-k-1)} \mathopen{}\originalleft(1+\bar{\ve{x}}^\prime\ve{S}_{x}^{-1}\bar{\ve{x}}\mathclose{}\originalright)^{-1} \widehat{\ve{\alpha}}^{\mathrm{OLS}\prime}\widehat{\ve{\Sigma}}_{u}^{\mathrm{OLS},-1}\widehat{\ve{\alpha}}^{\mathrm{OLS}}. \] Under the assumption that $\{\ve{u}_t\}_{t=1}^T$ is a sequence of independently and identically distributed multivariate normal random vectors, the GRS statistic follows an $F(N,T-N-k)$ distribution under the null hypothesis $H_0:\ve{\alpha}_0=\ve{0}$. Since these conditions exclude the serial dependence and dynamic regressor--disturbance dependence studied in this paper, the GRS test is included as a conventional benchmark rather than as a procedure justified under our general setting.\footnote{kamstra2024tra propose a modified version of the GRS test to improve its small-sample properties. However, the performance of the original and modified GRS tests does not differ substantially once the sample size $T$ exceeds 200. For this reason, we adopt the original GRS test as the benchmark.}

\paragraph{HAR test} When regression errors exhibit heteroskedasticity and autocorrelation, tests that do not account for these features are generally invalid. One standard approach is the heteroskedasticity and autocorrelation robust (HAR) test, which is based on a HAC estimator of the long-run covariance matrix; see newey1987sps,andrews1991hac. Let \[ \ve{\Omega} = \ve{\Gamma}_0+\sum_{j=1}^{\infty}\bigl(\ve{\Gamma}_j+\ve{\Gamma}_j^\prime\bigr), \qquad \ve{\Gamma}_j={\mathbb{E}}[\ve{v}_t\ve{v}_{t-j}^\prime], \qquad \ve{v}_t=\ve{Z}_t\ve{u}_t, \] and estimate $\ve{\Omega}$ by the kernel-based HAC estimator \[ \widehat{\ve{\Omega}}(M) = \widehat{\ve{\Gamma}}_0 +\sum_{j=1}^{T-1} k\mathopen{}\originalleft(\frac{j}{M}\mathclose{}\originalright) \bigl(\widehat{\ve{\Gamma}}_j+\widehat{\ve{\Gamma}}_j^\prime\bigr), \qquad \widehat{\ve{\Gamma}}_j = \frac{1}{T}\sum_{t=j+1}^T \widehat{\ve{v}}_t\widehat{\ve{v}}_{t-j}^\prime, \] where $\widehat{\ve{v}}_t=\ve{Z}_t\widehat{\ve{u}}_t$. The corresponding HAR Wald statistic for testing $H_0:\ve{R}\ve{\kappa}=\ve{r}$ is \[ \mathcal{W}^{\mathrm{HAR}} = T(\ve{R}\widehat{\ve{\kappa}}^{\mathrm{OLS}}-\ve{r})^\prime \mathopen{}\originalleft[ \ve{R}\widehat{\ve{Q}}_{Z}^{-1} \widehat{\ve{\Omega}}(M) \widehat{\ve{Q}}_{Z}^{-1}\ve{R}^\prime \mathclose{}\originalright]^{-1} (\ve{R}\widehat{\ve{\kappa}}^{\mathrm{OLS}}-\ve{r}), \] where $\widehat{\ve{Q}}_{Z}=T^{-1}\sum_{t=1}^T\ve{Z}_t\ve{Z}_t^\prime$. Under the conventional small-$b$ asymptotics, if $M\to\infty$ and $M/T\to0$, then under standard regularity conditions, $\mathcal{W}^{\mathrm{HAR}} \xrightarrow{d} \chi_q^2$ with $q=\operatorname{rank}(\ve{R})$.

A key practical issue in HAC-based inference is the choice of the bandwidth parameter $M$. In small samples, conventional HAR tests can be highly sensitive to this choice and often exhibit substantial over-rejection. To address this difficulty, kiefer2000srt propose the fixed-$b$ approach, setting $M=bT$ with $b\in(0,1]$. In this case, the HAC estimator is no longer consistent, and the HAR statistic has a nonstandard null distribution that depends on the kernel and on $b$. This approach often yields more accurate finite-sample size control than conventional HAC-based HAR tests, i.e., kiefer2002har,jansson2004erp,kiefer2005nat,sun2014lfa.

Building on these insights, lazarus2018hir,lazarus2021spt provide a comprehensive comparison of fixed-$b$ methods, establish the size-power frontier, and propose loss-function-based rules for selecting $b$. They show that the frontier is achieved by the quadratic spectral (QS) kernel. Against this background, in our simulations we adopt the fixed-$b$ HAR test with the QS kernel recommended by lazarus2018hir,lazarus2021spt, so that inference is based on the corresponding nonstandard fixed-$b$ asymptotic distribution rather than on the conventional small-$b$ HAC approximation.

Rejection frequencies under the null hypothesis

We first examine whether the empirical rejection frequency of each test under the null hypothesis is close to the nominal significance level. We then turn to size-adjusted power, which measures how well each test detects false null hypotheses after correcting for size distortions. Taken together, these results help clarify the trade-off between size control and power across the three DGPs: $BD$, $GEXOG$, and $EBD$.

Throughout this subsection, $\widetilde{\mathcal{W}}^{\mathrm{GD}}$ denotes the bootstrap Wald test based on the GD estimator, while $\mathcal{W}^{\mathrm{GD}}$, $\mathcal{W}^{\mathrm{FD}}$, and $\mathcal{W}^{\mathrm{FGCO}}$ denote the asymptotic Wald tests based on the GD, FD, and FCO estimators, respectively. In addition, $\mathcal{W}^{\mathrm{HAR}}$ and $\mathrm{GRS}$ denote the HAR and GRS tests, respectively. Under the null hypothesis $H_0:\ve{\alpha}_0=\ve{0}_N$, we report, for each test, empirical rejection frequencies at the nominal 10%, 5%, and 1% significance levels.

Table (ref) reports empirical rejection frequencies under the null hypothesis for the $BD$ condition. Overall, the bootstrap Wald test based on the GD estimator, $\widetilde{\mathcal{W}}^{\mathrm{GD}}$, yields rejection frequencies close to the nominal levels across all reported designs, even when the sample size is as small as $T=100$, indicating satisfactory finite-sample size control. By contrast, the asymptotic Wald tests based on GD, FD, and FCO exhibit rejection frequencies well above the nominal levels in small samples, although these distortions gradually diminish as the sample size increases. The rejection frequencies of $\mathcal{W}^{\mathrm{FD}}$ and $\mathcal{W}^{\mathrm{FCO}}$ are also very similar throughout, consistent with their asymptotic equivalence under $BD$ established in Proposition (ref). The HAR test also performs reasonably well, especially when $N=5$, although for $N=10$ it tends to under-reject in small samples. By contrast, the GRS test exhibits non-negligible size distortions. In particular, in the reported designs, it tends to over-reject the null hypothesis when $N=5$ and to under-reject it when $N=10$. This pattern is consistent with the fact that the GRS test is derived under the assumption that the disturbance vector is independently and identically distributed as multivariate normal, whereas under the $BD$ condition the disturbances are serially correlated.

Table (ref) reports empirical rejection frequencies under the null hypothesis when the $GEXOG$ condition holds. As in the $BD$ case, the bootstrap Wald test based on the GD estimator, $\widetilde{\mathcal{W}}^{GD}$, continues to deliver rejection frequencies close to the nominal levels across all reported designs, regardless of the sample size $T$. In addition, the asymptotic Wald tests based on GD and FD tend to move closer to the nominal levels as the sample size increases. By contrast, the Wald test based on the FCO estimator performs less satisfactorily, especially when the number of equations is $N=10$. This is because the FCO estimator is not consistent under $GEXOG$, as shown in Proposition (ref). This tendency becomes more pronounced as the number of equations increases. The most severe size distortions arise for the HAR and GRS tests. This pattern is unsurprising, since both procedures rely on OLS-based estimation, whose consistency is not guaranteed under $GEXOG$. Taken together, the results indicate that under $GEXOG$, $\widetilde{\mathcal{W}}^{GD}$ provides the most reliable finite-sample size control among the competing procedures.

Table (ref) reports empirical rejection frequencies under the null hypothesis for the $EBD$ condition. The bootstrap Wald test based on the GD estimator, $\widetilde{\mathcal{W}}^{\mathrm{GD}}$, continues to yield rejection frequencies close to the nominal levels across all reported designs, indicating stable finite-sample size control even under the most general dependence structure considered in this paper. The asymptotic Wald test based on the GD estimator also improves as the sample size increases, although it exhibits noticeable over-rejection in smaller samples, especially when $N=10$. By contrast, the Wald tests based on the FD and FCO estimators exhibit substantial size distortions throughout, and these distortions do not vanish as the sample size increases. This pattern is consistent with the fact that these estimators are not consistent under $EBD$. The HAR and GRS tests also exhibit severe size distortions under $EBD$, with rejection frequencies far above the nominal levels in many designs. Overall, these results are consistent with the theoretical prediction that, among the competing procedures, only the GD-based procedures remain valid under $EBD$.

Taken together, the results in Tables (ref)--(ref) show that the relative performance of the tests depends crucially on the underlying dependence structure. Under the more restrictive $BD$ and $GEXOG$ conditions, FGLS-based tests perform well when the assumptions required for their validity are satisfied. By contrast, across all three dependence structures, the GD-based Wald test, and especially its bootstrap version, delivers the most reliable finite-sample size control.

Rejection frequencies under the alternative hypothesis

We now examine the size-adjusted power of each test for the $BD$, $GEXOG$, and $EBD$ designs. Because finite-sample size distortions can substantially affect raw rejection frequencies, unadjusted power is not directly comparable across procedures. We therefore focus on size-adjusted power, so that differences in rejection frequencies reflect each test's ability to detect false null hypotheses rather than differences in size distortions. Specifically, we consider the alternative hypothesis $H_1:\alpha_1=0.15$ and $\alpha_j=0$ for $j\neq 1$.

Figure (ref) reports the size-adjusted rejection frequencies under the alternative hypothesis for the $BD$ design. Among the competing procedures, the FGLS-based tests, especially $\mathcal{W}^{\mathrm{FD}}$ and $\mathcal{W}^{\mathrm{FCO}}$, generally attain the highest rejection frequencies across the reported designs. This pattern is consistent with the theoretical prediction that, under $BD$, FCO is asymptotically efficient (Proposition (ref)) and FD is asymptotically equivalent to FCO (Proposition (ref)). Accordingly, once size distortions are removed, the stronger performance of these tests can be attributed to genuine efficiency gains under correct specification rather than to spurious over-rejection under the null. In addition, $\widetilde{\mathcal{W}}^{\mathrm{GD}}$ and $\mathcal{W}^{\mathrm{GD}}$ track each other very closely in all panels, indicating that bootstrap size correction entails little, if any, loss of power under $BD$.

By contrast, $\mathcal{W}^{\mathrm{HAR}}$ is generally less powerful than the system-based procedures. This pattern is broadly consistent with the simulation evidence in baillie2024rit, who show in a single-equation time-series regression setting that HAR-based inference can be less powerful than procedures that exploit additional dynamic structure. In addition, the GRS test exhibits the lowest power across all reported designs. This finding is also consistent with earlier studies. In particular, gungor2013tlf show that the power of the GRS test does not necessarily increase monotonically with the number of test assets and may even decline when too many assets are included. More generally, because the exact finite-sample validity of the GRS test relies on independently and identically distributed multivariate normal disturbances, violations of these assumptions can further weaken its finite-sample performance; see, for example, affleck1989nta.

Figure (ref) reports the size-adjusted power under the $GEXOG$ condition. For this design, we report size-adjusted power only for procedures that remain theoretically valid. Accordingly, we exclude $\mathcal{W}^{\mathrm{FCO}}$ and the OLS-based benchmark tests, namely $\mathcal{W}^{\mathrm{HAR}}$ and the GRS test. This exclusion reflects not only their severe size distortions under the null, but also the more fundamental point that these procedures are not valid under $GEXOG$; see Propositions (ref) and (ref).

Among the remaining procedures, several patterns emerge. First, the size-adjusted power of all three tests increases steadily with the sample size $T$, as expected. Second, $\mathcal{W}^{\mathrm{FD}}$ attains the highest rejection frequencies throughout the reported designs. This finding is consistent with the fact that, under $GEXOG$, the FD estimator may be asymptotically efficient; see Corollary (ref). Third, $\widetilde{\mathcal{W}}^{\mathrm{GD}}$ and $\mathcal{W}^{\mathrm{GD}}$ remain very close in all panels, indicating that bootstrap size correction entails little, if any, loss of power under $GEXOG$. As $T$ increases, the differences across the three procedures narrow, and their rejection frequencies approach one.

Figure (ref) reports the size-adjusted power under the $EBD$ condition. For this design, we report size-adjusted power only for the GD-based procedures, since the remaining tests are not valid under $EBD$. Two patterns stand out. First, the size-adjusted power of both $\widetilde{\mathcal{W}}^{\mathrm{GD}}$ and $\mathcal{W}^{\mathrm{GD}}$ increases steadily with the sample size $T$ in all reported designs, and approaches one as $T$ becomes large. Power is somewhat lower in the designs with $k=4$ when $T$ is small, although this difference narrows rapidly as $T$ increases. Second, the two procedures remain very close in every panel, with only negligible differences throughout. This indicates that the bootstrap implementation improves finite-sample size control under $EBD$ without any meaningful loss of power relative to the asymptotic GD test.

The simulation results show that the finite-sample performance of the competing tests depends crucially on the underlying exogeneity condition. Under the more restrictive $BD$ and $GEXOG$ conditions, procedures that exploit stronger structural restrictions, such as $\mathcal{W}^{\mathrm{FD}}$ and, when valid, $\mathcal{W}^{\mathrm{FCO}}$, can attain higher size-adjusted power. By contrast, the OLS-based tests provide meaningful inference only under $BD$; under more general dependence structures, their size distortions become severe and power comparisons based on them are no longer informative. As the dependence structure becomes more general, robustness becomes more important than efficiency gains available under more restrictive specifications. In this respect, the GD-based procedures remain stable across the $BD$, $GEXOG$, and $EBD$ designs. In particular, the bootstrap GD test consistently delivers reliable size control while preserving power close to that of the corresponding asymptotic GD test. These findings highlight the practical usefulness of the bootstrap GD test when the true exogeneity condition is unknown or may be misspecified.

Empirical Application

We examine whether the Wald test based on the generalized Durbin estimator proposed in this paper can mitigate the over-rejection problem observed for Wald tests based on the FD and FCO estimators. In Section (ref), we show that, under the $EBD$ condition, the generalized Durbin estimator is consistent, whereas the FD and FCO estimators are not. At the same time, our simulation results in Section (ref) indicate that the Wald test based on the generalized Durbin estimator is not immune to finite-sample size distortions. We therefore employ a bootstrap bias correction to improve finite-sample size control.

To illustrate the empirical relevance of our approach, we revisit the Fama--French multifactor models studied by fama1993crf,fama2015ffa. As noted by moriya2024tif, there remains ongoing debate over the validity of the FF multifactor models. In empirical studies (e.g., fama1993crf,fama2015ffa,fama2016daf,fama2017itf,fama2018cf,fama2020ccs; cakici2013svm), it has become standard practice to apply the GRS test to assess the validity of FF multifactor models. However, our simulation evidence in Section (ref) shows that both the GRS test and the HAR test can suffer severe size distortions when the $BD$ condition fails. Because $BD$ is more restrictive, and thus less empirically plausible, than $EBD$, and because the GRS and HAR tests are consistent under $BD$ but not under $EBD$, we do not employ these procedures in our empirical analysis. Instead, we compare the bootstrap bias-corrected Wald test based on the generalized Durbin estimator with Wald tests based on the FD and FCO estimators.

Fama--French Multifactor Models

Following fama1993crf,fama2015ffa, we introduce the FF multifactor models. Equation ((ref)) presents fama1993crf's fama1993crf three-factor (hereafter referred to as FF3) model, which is the most widely known multifactor model:

where $R_{i,t}$ denotes the return on the portfolio $i$ at time $t$, $R_{f,t}$ is the risk-free rate at time $t$, $R_{m,t}$ is the returns on the market portfolio at time $t$, and $\varepsilon_{i,t}$ is the error term for portfolio $i$ at time $t$. The FF3 model expands the CAPM by adding size and value risk factors to capture market anomalies. The size risk factor ($SMB_t$) reflects the empirical regularity that stocks with smaller (or larger) market capitalizations tend to earn higher returns, while the value risk factor ($HML_t$) accounts for the superior performance of stocks with low (or high) price-to-book ratios. Building on the FF3 model, fama2015ffa propose a five-factor model (hereafter referred to as FF5), as shown in Equation ((ref)):

This model adds two risk factors to FF3: the profitability risk factor ($RMW_t$) and the investment risk factor ($CMA_t$). The models are said to correctly capture the behavior of stock returns if all intercept terms $\alpha_i$ for each portfolio in Equations ((ref)) and ((ref)) are not significantly different from zero. Therefore, to assess the validity of the FF multifactor models, we test the following null hypothesis: \[ H_0: \ve{\alpha}=\ve{0},\ \ H_1: {\rm{not}} \ H_0 \] where $\ve{\alpha}=\mathopen{}\originalleft\{\alpha_i\mathclose{}\originalright\}_{i=1}^N$ in each of Equations ((ref)) and ((ref)).

Dataset

fama1993crf,fama2015ffa provide the returns on the benchmark portfolios and risk factors for examining the FF multifactor models. In these models, the benchmark portfolios are formed by classifying all stocks in the market into either $2\times 3=6 \ {\rm{or}} \ 5\times 5=25$ categories using two criteria (e.g., “market capitalization (firm size) and book-to-market ratio (MB)” as proposed by fama1993crf). The average return of each portfolio is regarded as the representative portfolio in the market. Other sorting criteria have also been proposed, such as “market capitalization and profitability (MO)” and “market capitalization and investment growth rate (MI)” (fama2015ffa).

lewellen2010asa point out that many studies test the FF multifactor models using portfolios sorted by “market capitalization (firm size)” and “book-to-market ratio,” but the models may not perform as well when using portfolios sorted by other criteria. To consider the possibility that different sorting methods may affect the results, we use three types of benchmark portfolios to compare the performance of the Wald test based on the generalized Durbin estimator proposed in this paper with the GRS test.

Furthermore, the FF multifactor models are often rejected when U.S. data are used, as documented in fama2012svm,fama2017itf. In this study, following much of the previous literature, we use U.S. monthly data to examine whether such rejections may be driven by the use of unsuitable test statistics, such as the GRS statistic. All datasets are available from \href{http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/data_library.html}{Professor Kenneth French's website}. For the U.S. monthly data used in this study, the available sample spans July 1990 to December 2025. Our estimation period is restricted to the post-global financial crisis (GFC) period, from October 2008 (the month following the Lehman Brothers bankruptcy) to December 2025. We adopt this window because the estimators we employ are time-invariant and do not explicitly accommodate structural breaks.

A large literature documents that the GFC coincided with shifts in equity-market comovements and factor sensitivities. In particular, bekaert2014gce report crisis-period changes in CAPM betas and increases in residual correlations around the Lehman episode. In addition, lehkonen2015smi document regime-dependent changes in international stock-market integration during the GFC, suggesting that coefficients in FF multi-factor models may have shifted in that period. While the COVID-19 pandemic is widely viewed as a potential source of structural change in global equity markets, ndako2025sbg report that, using bai1998etl's (bai1998etl,bai2003cam,bai2003cvm) multiple structural break tests, no additional break aligned with the onset of COVID-19 is detected across the countries and regions they analyze. Instead, their evidence points to the GFC as the dominant common event driving structural change in global stock market, and they conclude that the GFC generated the most pronounced breaks in those markets. Accordingly, our baseline sample spans October 2008 to December 2025.

Table (ref) shows the descriptive statistics and the results of the unit root test for the returns on risk factors in the FF multifactor models. The descriptive statistics show that for all datasets, there are no inconsistent data in the perspective of the traditional mean-variance approach. In the estimations, each variable that appeared in the moment conditions should be stationary. We apply the augmented Dickey--Fuller (ADF) test of said1984tur to check whether the variables satisfy the stationarity condition. The ADF test rejects the null hypothesis that each variable contains a unit root at the 5% significance level.\footnote{The ADF test also rejects the null hypothesis that the returns on the 6 and 25 benchmark portfolios contain a unit root at the 5% significance level, regardless of portfolio sorting method.}

Empirical Results

In this subsection, we evaluate the FF multifactor models using the Wald test based on the generalized Durbin estimator and compare its performance with competing Wald tests based on the FD and FCO estimators. To improve finite-sample size control, we also report a bootstrap bias-corrected Wald test based on the generalized Durbin estimator. Overall, Table (ref) shows that inference on the null hypothesis $H_0:\ve{\alpha}=\ve{0}$ depends materially on the estimator underlying the Wald statistic. In particular, the competing Wald tests based on the FD and FCO estimators, together with the conventional Wald test based on the generalized Durbin estimator, tend to produce smaller $p$-values and hence more frequent rejections than the bootstrap bias-corrected Wald test based on the generalized Durbin estimator, especially when the number of benchmark portfolios is large.

Table (ref) reports the test results for the FF3 and FF5 models. For the FF3 model, the evidence against the null hypothesis is noticeably stronger for the MB25 portfolios than for the MB6 portfolios. In the MB6 case, the bootstrap bias-corrected Wald test based on the generalized Durbin estimator yields $p=0.03$, whereas the other Wald statistics yield $p=0.01$. Thus, the non-bootstrap procedures imply stronger evidence against the null. For MB25, all reported statistics reject the null hypothesis at conventional significance levels; however, the bootstrap bias-corrected Wald test remains distinctly less aggressive, with $p=0.01$, whereas the other procedures yield $p$-values that are essentially zero. These results indicate that, as the cross-sectional dimension increases, the strength of evidence against the FF3 model becomes increasingly sensitive to the construction of the Wald statistic.

For the FF5 model, the results reveal a similar, though more nuanced, pattern across alternative portfolio sorts. For portfolios sorted by MB, the bootstrap bias-corrected Wald test does not reject the null at the 5% level for MB6, although the other Wald statistics yield $p=0.02$. For MB25, all procedures reject the null, but once again the bootstrap bias-corrected Wald test is less aggressive than its competitors. For portfolios sorted by MI, none of the tests rejects the null for MI6, with all reported $p$-values lying comfortably above conventional significance thresholds. By contrast, for MI25, the bootstrap bias-corrected Wald test yields $p=0.06$ and therefore does not reject at the 5% level, whereas the other Wald statistics yield $p$-values that are essentially zero and reject decisively. A similar pattern emerges for portfolios sorted by MO. For MO6, none of the statistics rejects the null, whereas for MO25 the bootstrap bias-corrected Wald test continues not to reject the null, with $p=0.19$, even though the competing Wald tests yield much smaller $p$-values and reject at conventional levels.

Taken together, these empirical findings suggest that the bootstrap bias-corrected Wald test based on the generalized Durbin estimator provides more conservative and comparatively stable inference across portfolio sorts and is less prone to frequent rejection in higher-dimensional systems. By contrast, the competing Wald tests based on the FD and FCO estimators, as well as the conventional Wald test based on the generalized Durbin estimator, often generate substantially more aggressive rejections in the same settings. This contrast is consistent with our simulation evidence that misspecification in the dynamic dependence between regressors and errors can induce size distortions in competing procedures, whereas the generalized Durbin approach, when combined with bootstrap bias correction, provides more reliable inference for multi-equation specification testing.

Conclusion

This paper proposes a new model specification test for multiple-equation systems in the presence of cross-equation error dependence and dynamic regressor--error dependence. Conventional specification tests for such systems often rely on exogeneity conditions strong enough to ensure consistency of the OLS estimator. In the presence of dynamic regressor--error dependence, these exogeneity conditions may fail, in which case the OLS estimator becomes inconsistent and OLS-based model specification tests no longer deliver valid inference.

To address these limitations, we clarify the relationships among alternative exogeneity conditions and characterize how the validity of competing estimators depends on the underlying dependence structure in multiple-equation systems. Building on this analysis, we propose a generalized Durbin (GD) estimator for multiple-equation systems with an intercept and cross-equation error dependence, and show that the resulting estimator remains consistent under the weakest exogeneity condition considered in the paper. We then derive its asymptotic distribution, construct Wald tests for model specification testing, and develop a bootstrap version to improve finite-sample inference.

The simulation results show that the finite-sample performance of the competing multivariate tests depends crucially on the underlying exogeneity condition. Under the more restrictive conditions such as $BD$ and $GEXOG$, tests based on stronger structural restrictions can perform well and may achieve higher size-adjusted power when correctly specified. However, once those restrictions fail, their finite-sample performance can deteriorate substantially. By contrast, the GD-based procedures remain stable across the exogeneity conditions considered in the paper. In particular, the bootstrap GD test delivers the most reliable finite-sample size control, even when the sample size is as small as $T=100$, while preserving sufficient size-adjusted power.

Our empirical application further illustrates the practical relevance of these findings. We examine the validity of the fama1993crf's (fama1993crf,fama2015ffa) multifactor models and find that inference on their validity depends crucially on the choice of test statistic, especially when the number of portfolios is large. In particular, the bootstrap Wald test based on the generalized Durbin estimator often leaves the null hypothesis unrejected in cases where the competing Wald tests reject it more aggressively. These results suggest that empirical rejections based on more restrictive procedures may partly reflect finite-sample distortions or misspecification of the underlying dependence structure, rather than decisive evidence against the asset pricing model itself. This contrast is consistent with our simulation evidence that misspecification in the dynamic dependence between regressors and errors can induce size distortions in competing tests.

Overall, the proposed bootstrap GD test provides a reliable and practically useful tool for multivariate model specification testing in the presence of cross-equation dependence and dynamic regressor--error dependence. When the true exogeneity structure is unknown, it offers a robust alternative to conventional procedures whose validity depends on stronger maintained assumptions.

Acknowledgments

The authors thank Tirthatanmoy Das, Katsuhito Iwai, Akitada Kasahara, Genya Kobayashi, Daisuke Nagakura, Yohei Yamamoto, Tatsuma Wada, and participants at the 100th Annual Conference of the Western Economic Association International and the Japan Society of Monetary Economics 2025 Autumn Meeting for helpful comments and suggestions. The author also gratefully acknowledges financial support from the Japan Society for the Promotion of Science, Grant-in-Aid for Scientific Research (Grant Nos. 23H00838 and 23K25535), and from the Japan Science and Technology Agency, Moonshot Research and Development Program (Grant No. JPMJMS2215). All data and programs used in this paper are available upon request.

\setcounter{table}{0}