Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
92,622 characters · 20 sections · 62 citation commands
Identification and Estimation of Simultaneous Equation Models Using Higher-Order Cumulant Restrictions
\smallKeywords: simultaneous equations, factor model, VAR, higher-order cumulants, independent component analysis
\smallJEL: C10, C30, C32, C38
Linear simultaneous-equation models are among the most commonly used tools in economics. Their appeal lies in their ability to capture equilibrium relationships and other scenarios where variables are determined jointly. However, as illustrated by the canonical “supply and demand” system, the structural equations that encode these relationships are typically not identified.
As HAUSMAN1983 discusses, the predominant solution to this identification challenge has been the instrumental variables (IV) approach. Yet, finding valid instruments can be difficult, as it often requires extensive structural modeling, deep subject-matter knowledge, and considerable ingenuity. This difficulty motivates a natural question: Is there a pathway forward when instruments are unavailable?
Such concerns have appeared frequently in econometric applications. A classic non-identification result arises in the vector autoregression (VAR) literature: if no further restrictions are imposed and one assumes a simple setting in which structural errors are uncorrelated and normally distributed, then the distribution of reduced-form errors---completely characterized by their mean and covariance---is insufficient to identify the structural parameters (see LANNE2017288 for a detailed argument). A common remedy involves exploiting higher-order information, such as identification from changes in second moments (heteroskedasticity), which has spawned a rich body of literature proposing influential identification strategies (e.g., 10.1162/003465303772815727,LANNE2010121; see annurev:/content/journals/10.1146/annurev-economics-070124-051419 for a comprehensive review). Here, we pursue higher-order moments in a more direct manner---namely, by moving away from the Gaussian assumption and focusing on distributions that encode additional information in their higher moments, such as skewness and kurtosis. Pioneering work in this direction can be traced back to bbfc08fa-3d75-37e2-b3fe-02d4101a7538 and 22ff115f-e119-3e81-8033-1597bfc7003b, who showed that, in errors-in-variables settings, non-Gaussianity can deliver identification without instruments.
Beyond econometrics, the assumption of mutually independent, non-Gaussian errors has also been widely investigated in signal processing. Under this framework, COMON1994287 developed one of the most influential blind source separation methods, later termed independent component analysis (ICA). In the simultaneous-equations context, ICA theory implies that if the structural errors are mutually independent and non-Gaussian (with at most one Gaussian component), then the structural parameter matrix (which defines the linear relationships across equations) is identified up to signed permutation and scale. Although this does not provide full point identification in every case, it is still a powerful starting point that, coupled with basic and relatively uncontroversial economic theory, can deliver exact identification.
In principle, ICA identification results translate into viable estimators. One prominent example is the FastICA algorithm introduced by 761722, which has been widely applied in fields ranging from telecommunications to medical imaging hyvarinen2013independent. Nevertheless, within econometrics---especially in microeconometrics---ICA-based estimators have seen limited use. One reason is that the statistical properties of these estimators are often not the central focus of ICA research, with 10.1214/009053606000000939 standing out as a notable exception. A more fundamental modeling concern arises from how structural errors are typically viewed in econometrics: rather than “signals” broadcasting independently, they are usually conceived as latent factors, which makes the strong assumption of mutual independence seem less realistic.
bonhomme2009consistent are among the earliest econometric attempts to move beyond the strict independence setup in a linear factor model with additive error, allowing error terms that may be correlated to some degree. Working in this setting, their identification argument shows that second- to fourth-order information can identify part or all of the loading/structural matrix, with the identifiable portion shrinking as error correlation strengthens; they also propose a JADE-type joint-diagonalization estimator that is consistent and asymptotically normal. Subsequently, in the SVAR literature, lanne2021gmm and guay2021identification show that, while maintaining uncorrelated structural errors, identification follows from imposing diagonality of either the third- or fourth-order cumulant tensor.\footnote{Formal definitions appear in Section 2.} mesters2024non extend this line of work, still assuming uncorrelated structural errors, by allowing richer dependence structures beyond diagonal fourth cumulants.
A recurring theme in this literature is the assumption that structural errors are uncorrelated. The chief reason is that it enables whitening, a preprocessing step that transforms the variables so that the mixing/structural‑coefficient matrix becomes orthogonal. This whitening step underpins most existing identification proofs. In structural VAR (SVAR) settings the assumption is usually innocuous, because uncorrelated shocks are both standard and, to some extent, desirable. In many other simultaneous-equation settings, however, the assumption is overly restrictive: it rules out linear dependence among structural errors, including dependence induced by measurement error or omitted common shocks. For instance, if an omitted causal variable appears in two equations, it induces correlation between their structural errors, violating the assumption.
This paper's first contribution is to show that the zero-covariance assumption is unnecessary for structural-parameter identification when a diagonal condition on a single higher-order moment (cumulant) is available. In Section 5 we prove that, under a diagonal higher-order cumulant assumption, whitening is unnecessary: the identification problem reduces to a simple eigenvector problem. This argument is both simple and constructive: it leads directly to a straightforward sample-analogue estimator. Moreover, even in settings such as SVAR models—where relaxing the uncorrelatedness assumption is arguably less crucial—our framework still offers a transparent test of that uncorrelatedness within the same higher-order cumulant framework used in earlier work.
Our second contribution is to propose a conceptually and computationally simple estimator with desirable statistical properties. In Section 6, we establish that the sample analogue estimator derived from our identification argument is consistent and asymptotically normal (under standard regularity conditions). We also show how the identification framework simplifies both the asymptotic analysis and the implementation. Section 7 then evaluates the estimator’s finite-sample performance in Monte Carlo experiments, and we conclude by illustrating its practical usefulness through two empirical applications.
Let $X \in \mathbb{R}^{d}$ be a random vector. Its characteristic function is
The cumulant generating function is
Here each index $i_{1},\dots,i_{r}$ runs independently over $\{1,\dots,d\}$, so the inner sum contains $d^{r}$ terms.
which defines the $r$th-order cumulants. Collecting all such coefficients yields the $r$th-order tensor $C_{r}(X)\in(\mathbb{R}^{d})^{\otimes r}$.
For expositional clarity, the main text focuses on the third- and fourth order cases; extensions to any fixed order $r\ge 3$ are stated and proved in the appendix. Write $\mu_i:=\mathbb{E}[X_i]$ and $\mathrm{Cov}(X_a,X_b):=\mathbb{E}\!\big[(X_a-\mu_a)(X_b-\mu_b)\big]$. Then, assuming the relevant moments exist,
In the special case $\mathbb{E}[X]=0$, these reduce to
By convention, the entries of the third cumulant tensor $C_3(X)_{ijk}$ with $i=j=k$ are referred to as diagonal entries. Denoting the $i$th diagonal entry by $\kappa_3(X_i)$, we have \[ \kappa_3(X_i)=C_3(X)_{iii}=\mathbb{E}\!\big[(X_i-\mu_i)^3\big], \] the skewness of $X_i$. Similarly, the $i$th diagonal entry of the fourth cumulant tensor is \[ \kappa_4(X_i)=C_4(X)_{iiii}=\mathbb{E}\!\big[(X_i-\mu_i)^4\big]-3\Big(\mathbb{E}\!\big[(X_i-\mu_i)^2\big]\Big)^2, \] commonly called the excess kurtosis.\footnote{ In this paper, “skewness” and “excess kurtosis” denote the raw third and fourth cumulants; the variance-normalized versions are not used.}
Cumulants have many useful properties. The following list records those used in this paper.
This section introduces the setup for our baseline model. Although the assumptions used here are already weaker than those in much of the current literature, they can be relaxed further. For clarity, we postpone the discussion of possible relaxations until Section 5.
Let $S$ be an unobserved $d$-dimensional random vector of structural errors, where $S_q$ denotes its $q$-th component. We observe a $d$-dimensional random vector $X$ generated by the linear structural equation system \[ \Lambda X = S. \]
Our primary goal is to recover the structural parameter matrix $\Lambda$ using $n$ independent and identically distributed (i.i.d.) observations of $X$.
The model imposes the following restrictions on $\Lambda$ and a fixed $h$th \((h > 2)\) cumulant tensor of $S$:
Before delving into these assumptions, two observations are in order. First, write $\mu_X:=\mathbb{E}[X]$ and $\mu_S:=\mathbb{E}[S]$, so $\Lambda \mu_X=\mu_S$, and define the centered variables \[ X_c := X-\mu_X,\qquad S_c := S-\mu_S,\qquad \text{so that}\quad \Lambda X_c = S_c. \] In the upcoming sections, for estimation, we preprocess by subtracting the sample mean, and in the population, we work with $X_c$ and $S_c$. Because cumulants of order $h\ge 2$ are translation invariant, we have \[ C_h(X)=C_h(X_c),\qquad C_h(S)=C_h(S_c),\qquad \kappa_h(S)=\kappa_h(S_c), \] so working with the mean zero variables $X_c$ and $S_c$ yields simpler and clearer arguments without loss of generality for both identification and estimation. Second, although the current model representation lacks exogenous observable covariates, these can be incorporated through a simple partialling-out argument.
Assumption A1) is standard: the observable variables $X$ can be written as a linear combination of the structural errors $S$ via the matrix inverse $\Lambda^{-1} = A$. The matrix $A$ is commonly referred to as the mixing matrix. While this assumption is natural in simultaneous-equation models, it is less convenient in factor-model settings where the loading/mixing matrix is typically rectangular. In Section 5, we show how to relax A1 to require only that $A$ have full column rank.
Assumption A3) implies that the $h$th cumulant tensor $C_h(S)$ is diagonal; A2) additionally rules out zero diagonal entries. Note that the commonly used non‑Gaussian, mutually independent assumption for structural errors is stronger than what is imposed here: if all components of $S$ are mutually independent, then (by the standard additivity property of cumulants) all order-$r\ge 3$ cumulant tensors of $S$ are diagonal. Furthermore, if these components deviate from Gaussianity due to asymmetry (nonzero skewness) or heavy or light tails (nonzero excess kurtosis), the third or fourth cumulant will have nonzero diagonal entries. Many classical ICA algorithms and identification protocols based on higher-order information, either implicitly or explicitly, exploit these diagonal structures. In particular, the diagonality of the fourth cumulant, alongside zero covariance (i.e., a diagonal second cumulant), is often the central assumption guaranteeing identification.
As is evident from our assumptions, we depart from the conventional approach—which typically combines a diagonal covariance (second cumulant) with a diagonal order-$h$ cumulant—and show that diagonality of a single higher-order cumulant by itself is sufficient for both identification and estimation, without imposing any second-moment restriction. This relaxation is important not only because it grants greater modeling flexibility, but also because it implies that previously used specifications are overidentified and hence testable. For concreteness, we focus on the case $h=3$, with occasional references to the kurtosis case ($h=4$). There are two main reasons for this choice. First, the structure of theorems and proofs remains the same for any fixed $h>2$, but skewness ($h=3$) yields the cleanest and most interpretable representation. Second, the $h=3$ case benefits most from relaxing full independence, thereby enabling model setups that were previously unattainable.
Under $h=3$, the model assumptions specialize to:
Below, we list a few examples for which these assumptions hold, alongside their corresponding econometric applications.
Consistent with most studies that leverage higher‑order cumulant structures for identification, Assumptions A1--A3 at best ensure the identification of \(A\) (and hence \(\Lambda\)) only up to a permutation and scaling. A familiar illustration of this concept appears in the Gaussian VAR literature. Specifically, if one assumes independence across structural errors and imposes a triangular (or acyclic) system, then the Cholesky decomposition of the variance matrix can identify the contemporaneous interaction matrix only up to a permutation and scaling. \footnote{ Formally, conditional on a chosen ordering, the standard Cholesky factor with positive diagonal fixes a particular scale. In applied work it is common to re‑normalize (e.g., setting \(\operatorname{diag}(\Lambda)=\mathbf{1}\) or shock variances to unity) after the Cholesky step. Our discussion treats such conventional re‑scalings as part of the diagonal‑scaling indeterminacy.} Concretely:
Scale indeterminacy is a routine and generally minor issue in linear models; it is typically addressed by normalizing the diagonal of \(\Lambda\) to unity. While this normalization resolves scale indeterminacy, the situation becomes more intricate when permutation indeterminacy is present. Overcoming permutation indeterminacy is often the principal hurdle between the aforementioned identification result and point identification. Fortunately, economic theory can be instrumental in pinpointing the correct permutation. Once an economically meaningful permutation or labeling is established, a diagonal normalization can then be used to eliminate scale indeterminacy.
Although resolving permutation indeterminacy is not the central focus of this paper, we outline two practical examples of how one might obtain point identification in structural models, guided by economic reasoning:
A crucial advantage of these procedures is that, once a consistent estimator of \(\widetilde{A}\) is available, one can choose the correct permutation of its columns with probability tending to 1 as the sample size goes to infinity after estimation. In lewis2021identifying, rules like this are referred to as consistent labeling criteria. Theorem 4 in lewis2021identifying shows that the error introduced by the post‑estimation labeling step is asymptotically negligible: the permutation choice does not alter the first‑order asymptotic distribution of the estimator. Consequently, once an asymptotically normal estimator of the permuted matrix is available, standard delta‑method arguments yield valid inference for normalized structural parameters.
Throughout this paper, we focus on settings analogous to P1 and P2, in which economic theory supplies a consistent‑labeling criterion that eliminates both scale and permutation indeterminacies. In this context, once a consistent and asymptotically normal estimator of the permuted structural‑parameter matrix \(\widetilde{A}\) is available, inference on the structural parameters follows directly from Theorem 4 in lewis2021identifying. Consequently, our objective is to show that, under Assumptions A1–A3, the mixing matrix \(A\) is identified up to a column permutation and scaling, and that consistent, asymptotically normal estimators of \(\widetilde{\Lambda}\) can be constructed.
In this section, we state the main identification theorem under the baseline model for $h=3$. Although the argument extends to all $h>2$, the general proof differs from the $h=3$ case only by minor technical details. For clarity and brevity, we relegate the general proof to the Appendix. We then discuss some useful extensions of the baseline model.
Let $X_c := X - \mathbb{E}[X]$ denote the centered observables. By translation invariance of cumulants (orders $\ge 2$), $\kappa_3(w^\top X)=\kappa_3(w^\top X_c)$. Let $w$ be a $d$-dimensional column vector and consider the objective
We begin with a lemma concerning the Hessian of this objective function.
We now state and prove our main identification theorem:
Working in the coordinates of the latent variable \(S\) yields straightforward proofs. However, the crucial role of Lemma (ref) can be obscured by this simplicity. Its intuition is more transparent when viewed in the coordinates of the observed \(X\) through tensor contraction.
For centered data, the third-order cumulant tensor \(C_3(X_c)\) has entries \[ \big[C_3(X_c)\big]_{ijk}=\mathbb{E}\!\big[(X_c)_i (X_c)_j (X_c)_k\big], \] so for any \(w\in\mathbb{R}^d\), the expression \(C_3(X_c)[w,\cdot,\cdot]\) denotes its mode-1 contraction with \(w\), i.e., the \(d\times d\) matrix \[ \big[C_3(X_c)[w,\cdot,\cdot]\big]_{jk} =\sum_{i=1}^d w_i\,\big[C_3(X_c)\big]_{ijk} =\mathbb{E}\!\big[(w^{\top}X_c)\,(X_c)_j\,(X_c)_k\big]. \] Intuitively, we can view the third-order cumulant tensor as slices of \(d\times d\) matrices indexed by \(i\), and the contraction takes a linear combination of these slices to form a weighted sum, with \(w_i\) the weight on the \(i\)th slice.
From a straightforward application of the chain rule in the coordinates of \(X\), the Hessian \[ \nabla_w^2 Q(w)=6\,\mathbb{E}\!\big[(w^{\top}X_c)\,X_cX_c^{\top}\big] =6\,\big[C_3(X_c)[w,\cdot,\cdot]\big], \] is exactly this tensor-contraction operation with weight \(w\).
In the latent coordinates, \(X_c=AS_c\) and \(w^{\top}X_c=v^{\top}S_c\) with \(v=A^{\top}w\), so \[ \mathbb{E}\!\big[(w^{\top}X_c)\,X_cX_c^{\top}\big] = A\,\mathbb{E}\!\big[(v^{\top}S_c)\,S_cS_c^{\top}\big]A^{\top} = A\Big(\sum_{i=1}^d v_i\,\mathbb{E}[\,S_{ci}\,S_cS_c^{\top}\,]\Big)A^{\top}. \] Thus, contracting the third-order cumulant tensor for \(X\) is, up to a common congruence by \(A\), equivalent to contracting \(C_3(S_c)\) along the vector \(v\).
Because cumulants are multilinear and translation-invariant, under the diagonality of \(C_3(S)\) (Assumption A3$^\prime$), any contracted matrix \(C_3(S_c)[v,\cdot,\cdot]\) is diagonal, with diagonal entries proportional to \(\kappa_3(S_i)\,v_i\). Lifting back to the coordinates of the observables applies a congruence transformation by \(A\) to these diagonal matrices:
\[ \nabla_w^2 Q(w)=A\,\mathrm{diag}\!\big(6\,\kappa_3(S_i)\,v_i\big)\,A^{\top}\equiv A\,D_wA^{\top}, \] where the weights depend on the contraction vector \(w\) (equivalently, \(v=A^{\top}w\)). This Hessian device gives fine control over the collapse: the diagonal weights inherit the variation in \(w\), so different directions produce different diagonal scalings.
The main theorem then observes that two such contractions—at \(w_1\) and \(w_2\)—yield \(G(w_\ell)\coloneqq A D_{w_\ell}A^{\top}\) with a common congruence. Forming \[ H=G(w_2)^{-1}G(w_1)=\Lambda^{\top}D_{w_2}^{-1}D_{w_1}(\Lambda^{\top})^{-1} \] removes the congruence and produces a matrix similar to a diagonal one; its eigenvectors are exactly the rows of \(\Lambda=A^{-1}\), up to scaling and permutation. Choosing \(w_1\) from a distribution with a density ensures distinct eigenvalues almost surely, guaranteeing uniqueness of the eigendirections.
The above argument is inspired by the joint diagonalization used by bonhomme2009consistent and ICA joint-diagonalization variants, but it follows a different and simpler path. Identification here requires only two matrices, and—by Lemma (ref)—both arise solely from the third cumulant (or any fixed order \(h>2\)) via tensor contraction, so no fourth- or second‑order information is required, in contrast to bonhomme2009consistent. Consequently, there is no need to assemble large families of cumulant matrices, introduce off‑diagonal minimization objectives, or impose second‑moment diagonality to enable pre‑whitening cardoso1993blind,hyvarinen2013independent. Identification is also robust to non‑diagonality in other cumulants, and when researchers do wish to maintain assumptions on those cumulants (e.g., in SVAR applications), the restriction becomes a testable implication in our framework rather than a maintained assumption. Relative to non‑orthogonal JD methods that also avoid whitening, our construction does not rely on fixed‑point iterations 1011195, tuning‑intensive off‑diagonality criteria 4671095, or positive‑definiteness constraints tied to the covariance 10.5555/1005332.1016784. As we will show in Section 6, this simplicity leads to transparent large‑sample analysis for the sample‑analogue estimator.
The baseline identification result relies on two conditions: (i) every structural error has a nonzero third cumulant, and (ii) the structural parameter matrix $A$ is square and invertible. Some popular empirical settings in economics and statistics violate one or both of these assumptions. Throughout this subsection we continue to work with the centered observables $X_c:=X-\mathbb{E}[X]$ and the cumulant objective $Q(w)=\kappa_3(w^\top X)=\mathbb{E}[(w^\top X_c)^3]$, so that \[ G(w):=\nabla_w^{2}Q(w)=A D_w A^\top, \qquad D_w:=\operatorname{diag}\!\big(6\,\kappa_3(S_i)\,(A^\top w)_i\big), \] and we let $X\in\mathbb{R}^{d_1}$, $S\in\mathbb{R}^{d_2}$, $A\in\mathbb{R}^{d_1\times d_2}$ (so $d_1=d_2=d$ in the baseline case). As shown before, by translation invariance of cumulants (orders $\ge 2$), working with $X_c$ is without loss for identification.
For example, in independent component analysis (ICA) and other higher‑order‑cumulant methods it is common to allow one (or several) structural errors to exhibit zero skewness, in which case $G(w)=A D_w A^{\top}$ is singular for all $w$ and ordinary inversion is unavailable. Moreover, in many economic applications, even when only a subset of structural errors exhibits nonzero skewness, it is still desirable for the model to deliver partial information about the structural parameters (e.g., for the columns of $A$ associated with the skewed shocks). On the other hand, in factor models the loading matrix $A$ is typically tall ($d_{1}>d_{2}$), creating a similar noninvertibility problem because $\operatorname{rank}(G(w))\le d_2<d_1$ for all $w$.
By contrast, non‑Gaussian VAR applications often satisfy the baseline invertibility and nonzero skewness conditions, yet practitioners also impose the additional assumption that structural shocks are uncorrelated. This second‑moment information is valuable and can be exploited whenever available via $\mathrm{Var}(X)=A D_2 A^\top$ with diagonal $D_2$.
The remainder of this subsection develops three extensions that accommodate these practical scenarios:
\paragraph{Extension 1: Nonsquare $A$ with skewed factors}
A common setting for factor models is when the dimension of the observed variables exceeds the dimension of the latent factors. Concretely, let \(X=A\,S\) with \(A\in\mathbb{R}^{d_1\times d_2}\), \(d_1>d_2\), and \(A\) full column rank. Under A2\(^{\prime}\)–A3\(^{\prime}\) and for choices of \(w_2\) such that \((A^\top w_2)_i\neq 0\) for all \(i\) (e.g., \(w_2=\mathbf{1}\) with \(\sum_q A_{qi}\neq 0\) for all \(i\)), we have for any \(w\): \[ G(w):=\nabla_w^2 Q(w)=A D_w A^\top, \qquad D_w:=\operatorname{diag}\!\big(6\,\kappa_3(S_i)\,(A^\top w)_i\big). \] In the tall case (full column rank \(A\) and invertible \(D_{w_2}\)), \[ G(w_2)^{+}=(A D_{w_2} A^\top)^{+}=(A^+)^{\!\top}\,D_{w_2}^{-1}\,A^{+}, \] where $A^+$ denotes the Moore–Penrose inverse. Therefore \[ H \;:=\; G(w_1)\,G(w_2)^{+} = A\,\big(D_{w_1}D_{w_2}^{-1}\big)\,A^{+} = A\,\Delta\,A^{+}, \qquad \Delta=\operatorname{diag}\!\Big(\tfrac{(A^\top w_1)_i}{(A^\top w_2)_i}\Big)_{i=1}^{d_2}. \] It follows that \[ H\,A_{\cdot i}=\delta_i\,A_{\cdot i},\qquad \delta_i=\Delta_{ii}, \] so the \(d_2\) right eigenvectors of \(H\) associated with its nonzero eigenvalues are precisely the columns of \(A\), identified up to column scaling and permutation. The remaining \(d_1-d_2\) eigenvalues equal zero, with eigenvectors spanning \(\operatorname{col}(A)^\perp\). As in the main theorem, if \(w_1\) is drawn from any absolutely continuous distribution on \(\mathbb{R}^{d_1}\) (with \(w_2\) fixed as above), the \(\{\delta_i\}\) are almost surely pairwise distinct.
\paragraph{Extension 2: Identification when there are non‑skewed structural errors}
Suppose some components of the structural error have zero third cumulant. Let \(J:=\{i:\kappa_3(S_i)\neq 0\}\) and assume \(J\neq\varnothing\). Then, for every \(w\), the diagonal matrix \(D_w=\operatorname{diag}\!\big(6\,\kappa_3(S_i)\,(A^\top w)_i\big)\) has zeros on the coordinates \(J^c\), so \(G(w)=A D_w A^\top\) is singular and the ordinary inverse used in Theorem (ref) is unavailable.
In this case, what is identifiable depends on the target parameter. If our parameter of interest is the mixing matrix \(A\), we lose very little. For \(\ell=1,2\) we can write \[ G(w_\ell):=A_J D_{\ell,J} A_J^\top, \qquad D_{\ell,J}:=\operatorname{diag}(d_{\ell i})_{i\in J}, \quad d_{\ell i}:=6\,\kappa_3(S_i)\,(A^\top w_\ell)_i, \] where \(A_J\) collects the columns of \(A\) indexed by \(J\). Assuming \((A^\top w_2)_i\neq 0\) for all \(i\in J\) (so \(D_{w_2,J}\) is invertible), we are back in the tall case of Extension 1 with \(A_J\in\mathbb{R}^{d_1\times |J|}\) full column rank. Hence \[ H_{\text{mix}} \;:=\; G(w_1)\,G(w_2)^{+} = A_J\,\big(D_{w_1,J}D_{w_2,J}^{-1}\big)\,A_J^{+}, \] and the \( |J| \) right eigenvectors corresponding to its nonzero eigenvalues are exactly the columns \(\{A_{\cdot i}: i\in J\}\), identified up to column scaling and permutation (with almost‑sure eigenvalue distinctness for random \(w_1\)).
By contrast, identification of the structural parameter matrix \(\Lambda=A^{-1}\) is generally not possible from third‑cumulant information alone when \(|J|<d_2\). The third‑order cumulant family \(\{G(w):w\}\) depends only on \(A_J\) and is invariant to how the complementary columns \(A_{J^c}\) are chosen (as long as \(A=[A_J\;A_{J^c}]\) has full column rank). Different completions \(A_{J^c}\) lead to different inverses \(\Lambda\) but the same \(\{G(w)\}\). Thus, without additional structure (e.g. second‑moment restrictions such as uncorrelated structural errors), the rows of \(\Lambda\) cannot be identified when some components are non‑skewed.
\paragraph{Extension 3: Uncorrelated structural errors and an overidentification test}
In some econometric applications one may wish to impose that the structural errors are uncorrelated. Writing $X_c:=X-\mathbb{E}[X]$, this means the second cumulant (covariance) satisfies \[ \mathrm{Var}(X)\;=\;\mathbb{E}[X_c X_c^\top]\;=\;A\,D_2\,A^\top, \] where $D_2$ is diagonal with the structural variances on its diagonal.
\noindentIdentification remark (optional). If one were to exploit the second‑moment structure in estimation, one could replace $\bigl(\nabla_w^2 Q(w_2)\bigr)^{-1}$ by $\mathrm{Var}(X)^{-1}$ in Theorem (ref), yielding \[ H_\Sigma \;:=\; \mathrm{Var}(X)^{-1}\,\nabla_w^2 Q(w_1) \;=\;\Lambda^\top\bigl(D_2^{-1}D_{w_1}\bigr)(\Lambda^\top)^{-1}, \] so the eigenvectors are again the rows of $\Lambda$ (up to scaling/permutation).Note that with a diagonal second cumulant, the non-inverse problem does not arise, and hence we can recover all rows of $\Lambda$ corresponding to the distinct eigenvalues. For the test below, however, we deliberately estimate the demixing matrix using third‑cumulant information only, so that the second‑moment implication can be tested rather than imposed.
\noindentTestable implication (Testing joint diagonality of the second and the third cumulants). Let $\widetilde{\Lambda}$ denote the (scaled‑and‑permuted) demixing matrix recovered solely from third‑cumulant diagonalization, i.e. \[ \widetilde{\Lambda} \equiv \text{(rows of eigenvectors of }H:=(\nabla_w^2 Q(w_2))^{-1}\nabla_w^2 Q(w_1)\text{, oriented and normalized)}. \] Under A1$^\prime$--A3$^\prime$ plus uncorrelated structural errors, \[ \widetilde{\Lambda}\,\mathrm{Var}(X)\,\widetilde{\Lambda}^\top \;=\; P^* D^*\,\Lambda\,\mathrm{Var}(X)\,\Lambda^\top D^* P^{*\,\top} \;=\; P^* D^*\,D_2\,D^* P^{*\,\top} \] is diagonal. Equivalently, its unique off‑diagonal entries are zero in population, regardless of the unknown row scaling and permutation. Let $\operatorname{vech}_{\mathrm{off}}(\cdot)$ denote the operator that stacks the strict upper‑triangular entries of a symmetric matrix into a vector in $\mathbb{R}^{\binom{d}{2}}$. The overidentifying restrictions are \[ \operatorname{vech}_{\mathrm{off}}\!\Bigl(\,\widetilde{\Lambda}\,\mathrm{Var}(X)\,\widetilde{\Lambda}^\top\,\Bigr)\;=\;0\in\mathbb{R}^{\binom{d}{2}}. \]
\noindentSample statistic and asymptotics (outline). Estimate $\widetilde{\Lambda}$ from third‑order information only by \[ \hat H \;=\; \bigl(\nabla_w^2 \hat Q(w_2)\bigr)^{-1}\,\nabla_w^2 \hat Q(w_1), \quad \hat{\widetilde{\Lambda}}=\text{(oriented, normalized rows of eigenvectors of }\hat H\text{)}, \] and estimate the covariance by the centered sample covariance \[ \widehat{\Sigma}_X \;=\; \frac{1}{n}\sum_{i=1}^n (X_i-\bar X)(X_i-\bar X)^\top. \] Define the $\binom{d}{2}\times 1$ vector of unique off‑diagonals \[ \hat r \;=\; \operatorname{vech}_{\mathrm{off}}\!\Bigl(\,\hat{\widetilde{\Lambda}}\,\widehat{\Sigma}_X\,\hat{\widetilde{\Lambda}}^\top\,\Bigr). \] If an asymptotically normal estimator $\hat{\widetilde{\Lambda}}$ is available (as established in the Section 6; see Theorem (ref)), then, under mild moment conditions, a joint Delta‑method argument implies \[ \sqrt{n}\,\hat r \;\xrightarrow{d}\; \mathcal{N}(0,\Omega), \] for some positive definite $\Omega$ that can be consistently estimated (e.g., via a plug‑in Delta method based on raw moments up to order 6, or via jackknife/bootstrap). Consequently, the quadratic form \[ T_n \;=\; n\,\hat r^\top\,\hat\Omega^{-1}\,\hat r \] is asymptotically $\chi^2$ with $\binom{d}{2}$ degrees of freedom under the null. The construction is invariant to row scaling and permutation of $\hat{\widetilde{\Lambda}}$, so no additional normalization or row matching is required. Full details—covariance construction and the proof that $T_n\Rightarrow\chi^2_{\binom{d}{2}}$—are provided in the Appendix. Moreover, we show that under the usual regularity conditions for VAR models (e.g., those adopted in DAVIS2023180), the test statistic is asymptotically unaffected by the first-stage OLS estimation used to partial out the lags.
Given the constructive nature of the identification proof, a plug-in estimator is a natural choice. In this section, we establish the asymptotic properties of the simple plug-in estimator under the baseline model by showing that it is a smooth function of sample averages. Estimation procedures for the model extensions discussed in Section 5 are provided in the appendix.
Recall the population objective function (defined via the third cumulant) \[ Q(w) \;=\; \kappa_3\!\bigl(w^\top X\bigr) \;=\; \mathbb{E}\Bigl[\bigl(w^\top X_c\bigr)^3\Bigr], \qquad X_c := X - \mathbb{E}[X]. \] It is straightforward to see that the entries of its Hessian matrix are a linear combination of the entries of the third-order cumulant tensor; for order three, cumulants coincide with centered moments: \[ \nabla_w^2 Q(w) \;=\; 6\,\mathbb{E}\!\big[(w^\top X_c)\,X_c X_c^\top\big] \;=\; 6\sum_{r=1}^d w_r\,\mathbb{E}\!\big[(X_c)_r\,X_c\,X_c^\top\big]. \]
For asymptotics it is convenient to parameterize everything by raw moments and apply the cumulant map inside the moment-to-Hessian construction. Let \[ M(X)\in\mathbb{R}^{D_3(d)-1} \quad\text{stack all raw monomials up to total degree 3 (excluding the constant),} \] so that \(D_3(d)-1=\binom{d+3}{3}-1\). For example, when \(d=2\), \[ M(X)^\top = \bigl(X_1,\,X_2,\,X_1^2,\,X_1X_2,\,X_2^2,\,X_1^3,\,X_1^2X_2,\,X_1X_2^2,\,X_2^3\bigr). \]
We then define the cumulant map (still denoted \(C\)) \[ C:\mathbb{R}^{D_3(d)-1}\to\mathbb{R}^{d_3},\qquad c_3 := C(M), \] where \(d_3=\binom{d+2}{3}\) is the number of distinct entries in the (symmetric) third-order cumulant tensor. Intuitively, \(C\) takes raw moments and, via the standard moment–cumulant identities, subtracts off lower-order products to return third-order cumulants (which, at order three, equal centered third moments).\footnote{For any fixed order \(r\), the map from raw moments up to order \(r\) to cumulants up to order \(r\) is a multivariate polynomial}
The Hessian \(\nabla_w^2 Q(w)\) is an affine function of \(c_3\), so the population matrix \[ H(M) \;=\; \bigl(\nabla_w^2 Q(w_2)\bigr)^{-1}\,\nabla_w^2 Q(w_1) \;=\; H\!\big(C(M)\big) \] is an analytic function of \(M\) on any neighborhood where \(\nabla_w^2 Q(w_2)\) is nonsingular. Order the eigenvalues of \(H(M)\) in decreasing order and let \(u_k(M)\) denote the associated normalized eigenvector, oriented by one of the following conventions.
\noindentOrientation A (row-sum rule). Assume each row of \(\Lambda\) has nonzero sum, i.e.\ \(\mathbf{1}^\top \Lambda_{i\cdot}\neq 0\) for all \(i\). Fix the sign by requiring \[ \mathbf{1}^\top u_k(M) \;>\; 0 . \]
\noindentOrientation B (largest-entry rule). Alternatively, assume that for each \(k\) the largest absolute coordinate of the population eigenvector \(u_k(M)\) is unique. Let \(j_k=\arg\max_j |[u_k(M)]_j|\) and fix the sign by requiring \([u_k(M)]_{j_k}>0\).
Either convention yields a well-defined, smooth map \[ g:\mathbb{R}^{D_3(d)-1}\to\mathbb{R}^d,\qquad g(M)=u_k(M), \] in a neighborhood of the population \(M=\mathbb{E}[M(X)]\).
For the point estimator we use the natural centered sample analogue of the Hessian: \[ \nabla_w^2 \hat{Q}(w) \;=\; 6\,\frac{1}{n}\sum_{i=1}^n \bigl(w^\top X_{c,i}\bigr)\,X_{c,i} X_{c,i}^\top, \qquad X_{c,i}:=X_i-\bar X,\ \bar X:=\tfrac1n\sum_{i=1}^n X_i . \] Equivalently, one may compute \(\nabla_w^2 \hat{Q}(w)\) by first forming the vector of sample raw moments \(\hat M=\tfrac{1}{n}\sum_{i=1}^n M(X_i)\), then applying the cumulant map \(C\) to obtain the sample third cumulants \(C(\hat M)\), and finally plugging these into the affine formula for \(\nabla_w^2 Q(w)\); the two procedures coincide algebraically because third-order cumulants equal centered third moments.
Define the sample analogue of \(H\) by \[ \hat H \;=\; \bigl(\nabla_w^2 \hat{Q}(w_2)\bigr)^{-1}\,\nabla_w^2 \hat{Q}(w_1). \] Let \(\tilde u_k\) denote a (possibly complex-valued) normalized eigenvector of \(\hat H\) associated with its \(k\)th largest eigenvalue, oriented by the same rule used in population (Orientation A or B, applied after normalization). We then define our estimator as the real part \[ \hat{u}_k \;=\; \Re\bigl(\tilde u_k\bigr). \] This “real-part” safeguard is asymptotically inactive: at the population \(M\), the relevant eigenvalue is real and simple (Theorem (ref)), so there exists a neighborhood of \(M\) on which the eigenvector map is real-analytic and real-valued; within that neighborhood \(\Re(\tilde u_k)=\tilde u_k\). Moreover, \(\Re(\cdot)\) is a real-linear projection \(\mathbb{C}^d\to\mathbb{R}^d\), so composing with \(\Re\) preserves differentiability at \(M\).
For asymptotics we view \(u_k\) as a function of the raw moment parameter \(M\) and write \(g(M)=u_k(M)\). Our estimator is the plug-in \[ \hat{u}_k \;=\; g(\hat{M}), \qquad \hat M=\frac{1}{n}\sum_{i=1}^n M(X_i). \]
From the Continuous Mapping Theorem and the Delta Method, we obtain the main estimation result:
To determine the finite sample performance of our method, we design our simulation experiment following the composite structural error setup described in Section 3. This framework captures all three classical sources of endogeneity: measurement error, simultaneity, and omitted variable bias. The data-generating process is given by:
Here, \(X_1\) and \(X_2\) denote the observed counterparts of the latent regressors \(X_1^{*}\) and \(X_2^{*}\), each contaminated by additive measurement errors \(\epsilon_1\) and \(\epsilon_2\). We treat these errors as classical—mean zero and independent of all latent economic variables—while allowing them to be strongly correlated with one another. Specifically, \((\epsilon_1,\epsilon_2)\) is drawn from a bivariate normal distribution with unequal variances (marginal variances \(1\) and \(0.25\)) and a high negative covariance \(-0.45\) (implying correlation \(-0.9\)). The factor \(\sqrt{k}\) in Eq. (ref) scales the measurement‑ error variance, so larger \(k\) values monotonically worsen the signal‑to‑noise ratio and let us trace how estimator performance deteriorates as data quality declines. Importantly, our estimator does not require cross‑equation independence between the two measurement errors; it remains valid even when \(\operatorname{Cov}(\epsilon_1,\epsilon_2)\neq0\). We therefore impose a large negative covariance to stress‑test this robustness. A leading empirical case where such correlation occurs—and where the relaxation is essential—is the household‑budget data studied by 2a001323-07d6-31c1-910a-be0b0a4eecf7, in which expenditure and quantity are recorded while unit price is derived as their ratio. A positive recording error in quantity mechanically induces an equal‑and‑opposite error in the derived price, producing the negative correlation we mirror here (with \(X_1\) representing \(\log\)‑quantity and \(X_2\) representing \(\log\)‑price).
We generate the latent variables \(\{X_1^*, X_2^*\}\) from a simultaneous‑equations model with structural parameter \[ \Lambda=
, \] where the sign pattern replicates a supply‑and‑demand environment (downward‑sloping demand and upward‑sloping supply). The associated structural errors consist of two parts:
We assume mutual independence across the three blocks \((s_1,s_2)\), \((e_1,e_2,e_3)\), and \((\epsilon_1,\epsilon_2)\). The scaling \(\sqrt{k/3}\) in Eq. (ref) assigns total variance contribution of order \(k\) to the omitted‑shock component across its three independent sources (each contributes \(\approx k/3\)).
Both the measurement errors \(\epsilon\) and the omitted common effects \(\{e_1, e_2, e_3\}\) act as “noise,” inducing correlation in second and higher‑order moments. By varying \(k\), we directly control the magnitude of this noise and investigate its impact on estimation.
Our simulation focuses on recovering the demand elasticity \(-\Lambda_{12}=-1.5\). For reporting, we equivalently work with the positive slope parameter \(b_1\equiv \Lambda_{12}=1.5\) and compute MSE for \(\hat b_1\) relative to \(b_1\) (IV slopes are negated post‑estimation and signs are aligned so that \(\Lambda_{11},\Lambda_{22}>0\), \(\Lambda_{12}>0\), \(\Lambda_{21}<0\)).
Alongside our proposed method, we implement three other estimators for comparative benchmarking:
We consider sample sizes \(n\in\{500,3000,5000\}\) and conduct \(10{,}000\) Monte Carlo replicates for each \((n,k)\) and each estimator. Within each replicate we use a common random numbers design: a single “big” dataset is generated and then reused across \((n,k)\); the Fast‑ICA initialization is re‑seeded identically across \((n,k)\) to isolate design effects rather than algorithmic randomness. Point identification for both our proposed method and Fast‑ICA is achieved via sign restrictions.
We measure estimation accuracy using the mean squared error (MSE) between the true parameter \(b_1=1.5\) and its estimate \(\hat b_1\). Finally, we vary the noise‑to‑signal parameter \(k\in\{0,0.1,\dots,0.5\}\). Increasing \(k\) inflates the measurement‑error component \((\sqrt{k}\,\epsilon)\) and the omitted‑shock component \((\sqrt{k/3}\,e)\), thereby increasing the variance of the composite error. Because both variances and covariances scale proportionally in \(k\), the cross‑equation correlation of the composite error remains unchanged for any \(k>0\) (it is undefined at \(k=0\)).
From Table 1, we see that Fast‑ICA (F‑ICA) performs very well at \(k=0\), the independence case: its MSE is comparable to (and for larger \(n\), slightly below) that of the proposed eigenvector estimator (M1). Once \(k>0\), however, F‑ICA’s MSE grows sharply with \(k\) and, at higher noise levels, even fails to improve with larger \(n\) (e.g., at \(k=0.5\), MSE rises from \(0.916\) at \(n=500\) to \(1.249\) at \(n=5000\)), indicating asymptotic inconsistency under dependent composite shocks. The mechanism is subtle: although the fixed‑point iteration itself optimizes a higher‑moment criterion and does not use second moments directly, identification in Fast‑ICA relies on a whitening step that premultiplies the data by the inverse square root of the covariance, yielding \(z=\Sigma_x^{-1/2}x=B s\). This step forces the subsequent search over an orthogonal demixing matrix, which is valid only when the latent signals (here, the composite structural errors) are uncorrelated —a condition implied by independence in ICA, but violated in our design for \(k>0\). Consequently, no orthogonal demixer can recover the sources, and the algorithm converges to a biased limit even though its higher‑order objective is correctly specified within the (incorrect) orthogonal constraint set. This failure mode is representative: many higher‑cumulant‑based procedures impose whitening‑plus‑orthogonality as a precondition, and thus inherit a similar bias in this setup.
By contrast, M1’s MSE increases only modestly with \(k\) but declines sharply as \(n\) grows (e.g., from \(0.0113\) to \(0.00123\) at \(k=0\), and from \(0.0524\) to \(0.00439\) at \(k=0.5\)), consistent with the root‑\(n\) rate predicted by our asymptotic theory. Quantitatively, for \(k>0\) F‑ICA is about \(4\times\) to \(284\times\) less accurate than M1 across our \((n,k)\) grid.
Turning to the IV benchmarks, IV‑1 (oracle) uniformly dominates, as expected from using the latent shifter \(s_2\): its MSE is about half of M1’s at \(k=0\) and about one‑fifth to one‑quarter at \(k=0.5\) across \(n\). Importantly, although M1 does not use any instrument, it achieves the same order of magnitude MSE as IV‑2, which benefits from partial access to the latent shifter through \(\tilde z=\sqrt{0.3}\,s_2+\sqrt{0.7}\,z\). Across all \((n,k)\), the MSE ratio \(\text{M1}/\text{IV-2}\) lies between \(\approx 0.53\) and \(\approx 1.43\) (median \(\approx 0.83\)), with M1 typically outperforming IV‑2 for low‑to‑moderate noise (\(k\le 0.3\); e.g., \(n=5000,k=0.2\): \(0.00201\) vs \(0.00264\)) and IV‑2 modestly ahead at higher noise (\(k\ge 0.4\); e.g., \(n=5000,k=0.5\): \(0.00439\) vs \(0.00350\)).
Overall, this experiment shows that neglecting dependence in the composite errors is highly consequential for whitening‑based, higher‑order methods, whereas M1 remains reliable across finite samples and a wide range of noise levels—delivering near‑IV performance without any instrument.
Table 2 reports empirical coverage of two-sided 95% confidence intervals for the demand elasticity \(b_1=\Lambda_{12}=1.5\) (equivalently, \(-\Lambda_{12}\)). We use the same DGP as in the simulation study with the noise level fixed at \(k=0.5\), across sample sizes \(n\in\{500,3000,5000\}\) and 5{,}000 Monte Carlo replications.
As discussed in Section 6, the variances are estimated using two methods: the leave‑one‑out jackknife and the delta method. Confidence intervals are then constructed by asymptotic normal approximation using the standard normal critical value \(z_{0.975}\).
As shown in Table 2, both procedures achieve coverage close to the nominal 95% level. The jackknife is slightly closer to nominal at \(n=500\), while for \(n\ge 3000\) the two methods are essentially indistinguishable. With 5{,}000 replications, the Monte Carlo standard error of a 95% coverage estimate is about 0.003 (i.e., roughly \(\pm 0.6\) percentage points for a 95% simulation band), so small differences of that magnitude are not substantively meaningful. These findings are consistent with our asymptotic normality derivation and support the use of both plug‑in (delta) and resampling (jackknife) inference in this setting.
We study the finite-sample behavior of the over-identification test from Section 5 (Extension 3) using the same data-generating process and simulation protocol as in Section 7.1. In each replication we estimate the (scaled/permuted) structural parameter matrix \(\widehat{\widetilde{\Lambda}}\) solely from third-cumulant information via the eigenvector method, fixing \(w_2=\mathbf{1}\) and drawing a single \(w_1\) at random (uniformly from the unit cube, as in Theorem 5.1) once and for all. We then test the null of the joint validity of uncorrelated structural errors and A1$^\prime$–A3$^\prime$ by examining whether \(\widehat{\Theta}=\widehat{\widetilde{\Lambda}}\;\widehat{\operatorname{Var}}(X)\;\widehat{\widetilde{\Lambda}}^{\!\top}\) is diagonal. With \(d=2\) this yields a single over-identifying restriction; we use a Wald statistic with variance computed by delta‑method linearization and compare to a \(\chi^2_1\) critical value at \(\alpha=5\%\).
Table (ref) reports rejection rates for sample sizes \(n\in\{500,750,1000,5000\}\) and noise scales \(k\in\{0,0.1,\dots,0.5\}\) (larger \(k\) produces stronger departures from joint diagonality). Under the null (\(k=0\)) the test modestly over-rejects at small samples (10.4% at \(n=500\)) but improves with \(n\) (6.0% by \(n=5000\)), indicating that the size distortion recedes as the sample grows. Power rises with both the severity of the violation and with \(n\): for instance, at \(k=0.2\) the rejection rate increases from 0.442 at \(n=500\) to 0.813 at \(n=1000\); by \(n=5000\) power is essentially one for all \(k\ge 0.2\) (and 0.985 even at \(k=0.1\)). At very low sample sizes the test still struggles to detect mild violations—performance at \(n=500\) remains moderate—so some VAR applications with monthly data may face limited power against small departures from joint diagonality. That said, the variance regime used here is deliberately challenging; in many empirical settings (e.g., when working with residuals after partialling out lags) the effective noise level is lower, which should improve performance.
Our first application estimates returns to schooling with the proposed method. Data come from the Wooldridge textbook empirical exercise; details are provided in the wooldridge package manual on CRAN.\footnote{\url{https://cran.r-project.org/web/packages/wooldridge/wooldridge.pdf}}
Following card1993using, we adopt the linear system
Here, \(lwage_i\) denotes the natural log of wages, \(educ_i^{obs}\) is reported schooling, and \(X_i\) denotes a vector of control variables—such as race, potential experience, and location dummies. Following card1993using, there are two main sources of inconsistency in the OLS estimate of \(\beta\). The first is measurement error in reported schooling. Consistent with the original study, we model this error as classical: mean‑zero, independent of the true education level and other regressors, and approximately normal (specifically, symmetric). Let \(educ_i^*\) denote true schooling and write the measurement equation \[ educ_i^{obs} = educ_i^* + \Delta_{educ,i}. \]
The second—and more challenging—issue is omitted variable bias. As emphasized by card1993using, much of the literature attributes this bias to unobserved “ability,” a latent individual trait that is inherently difficult to measure. In this paper, we assume that ability is the primary source of structural error correlation. That is, if ability were observed, all parameters in the model could be estimated consistently. Moreover, we assume that ability, conditional on \(X_i\), is symmetrically distributed. This assumption is consistent with precedents in applied econometrics (e.g., ee98fe21-acdf-39c2-ada6-7518abe3ad0b,doi:10.1086/504455 ) and is supported by empirical findings from psychometric studies (e.g., jensen1998g).
In summary, we use the following triangular linear structural model, which mirrors the composite structural error framework discussed in Section 3 and verified by simulation in Section 7: \[
\] where \(u_i\) and \(v_i\) are the structural errors after accounting for \(ability_i\). These errors capture factors with no cross‑equation effects independent of ability and measurement error. The term $v_i$ collects idiosyncratic shocks that affect schooling but not wages except through schooling (given $X_i$); in IV setups these shocks are the source of exogenous variation typically targeted by instruments. We assume \(u_i\) and \(v_i\) are mutually independent and have non-zero skewness. \(\Delta_{educ,i}\) denotes measurement error in education, assumed to be normally distributed and independent of the true level of education. We also assume \(ability_i\) is the sole source of omitted‑variable bias and is independent of \(u_i\) and \(v_i\) (conditional on \(X_i\)). Finally, once the effects of \(X_i\) are accounted for, the distribution of \(ability_i\) is assumed symmetric.
As shown in Table 2 of card2001estimating, previous estimates of the return to schooling range from $0.0245$ to $0.36$. card1993using IV estimate based on distance to college is $0.132$. Using our proposed estimator, we obtain a point estimate of $0.0987$ with a 95% jackknife confidence interval of \([0.0358,\,0.1500]\). This interval contains the majority of the estimates reviewed by card2001estimating in this setup, and our point estimate is close to the distance‑to‑college IV benchmark, lending credibility to our approach.
Identification through higher moments or cumulant restrictions has recently become more popular in the macroeconometrics literature, particularly within the vector autoregression (VAR) framework. As highlighted in the introduction, in this literature the assumption of uncorrelated structural shocks is not only relatively uncontroversial but also often desirable for interpretable impulse responses. Our results show that the commonly imposed model—characterized by a diagonal second cumulant and diagonal third or fourth cumulants of the structural shocks—is overidentified and therefore testable. One reason structural shocks can appear dependent in practice is that the specified VAR system is simply too small, so the included variables do not support a linear causal model. The empirical exercise below illustrates how the test proposed in Section 5 can detect such misspecification and how the resulting evidence can answer a substantive economic question.
Uncertainty typically rises in economic downturns. ludvigson2021uncertainty use a VAR to ask whether uncertainty helps cause recessions or is instead an endogenous response. A central conclusion of the paper is that the type of uncertainty matters: innovations to financial uncertainty (\(UF\)) behave more like an exogenous driver, raising macro uncertainty (\(MU\)) and lowering industrial production (\(IP\)), whereas downturns in \(IP\) have limited feedback to \(UF\).\footnote{See ludvigson2021uncertainty for details on the uncertainty measures and VAR specification.} In their analysis, different variables were considered to capture real (macro) uncertainty. The main specification used \(MU\) to obtain the aforementioned result, while the economic policy uncertainty (EPU) index constructed by 10.1093/qje/qjw024 was employed as a robustness check. Interestingly, in that specification \(UF\) exhibited a statistically significant contemporaneous response to \(IP\). However, the EPU sample there was much shorter (358 observations).
Importantly, the view that \(UF\) acts as a driver of the business cycle is not isolated. https://doi.org/10.1002/jae.2672 use a heteroskedasticity‑based identification method and conclude that financial uncertainty does not respond to shocks in real activity nor to shocks in macro uncertainty. DAVIS2023180 revisit the same question using an ICA‑based approach: they posit that the reduced‑form residuals are linear combinations of three independent non‑Gaussian shocks and develop a permutation‑based independence test under a proposed causal ordering. Within the \((MU,UF,IP)\) system, the lower‑triangular orderings that place \(UF\) on top (most exogenous) are not rejected, providing further support for ludvigson2021uncertainty. When \(MU\) is replaced by EPU, however, independence is rejected for all triangular orderings. This may reflect data limitations (the shorter EPU sample used in those papers), or it may indicate that a strict triangular contemporaneous structure does not hold (EPU also responds contemporaneously to \(IP\)). Building on the statistical framework of DAVIS2023180, we therefore re‑estimate the EPU specification using updated EPU data that cover the full ludvigson2021uncertainty sample, 1960:07–2015:04. Our goal is to assess the claim that \(UF\) is approximately exogenous by evaluating two conditions: (1) shocks to \(UF\) act as the exogenous driver; and (2) the shocks are uncorrelated and satisfy A2$^\prime$–A3$^\prime$.
We adopt DAVIS2023180's VAR set‑up \[ Y_t \;=\; A_1 Y_{t-1} + \cdots + A_p Y_{t-p} + e_t, \qquad e_t \;=\; B u_t, \] with \(Y_t'=(UF_t,EPU_t,IP_t)\) and \(u_t\) a three‑dimensional vector of shocks. We assume \(u_t\) are i.i.d.\ across \(t\), and we set the lag length to \(p=6\) as in ludvigson2021uncertainty. The structural object of interest is \(B\). Our test, however, differs from DAVIS2023180 (and, e.g., 10.1257/pandp.20221047): we test only the minimal restriction set that suffices for identification here—diagonality of a single higher‑order cumulant (third order) together with uncorrelated structural shocks. We implement the joint‑diagonality over‑identification test on the estimated reduced‑form residuals. We assume the lag order is correctly specified so that the Delta‑method large‑sample approximation for our Wald statistic (Section 5) is valid.
Ideally, one would implement a joint three‑equation test and, if not rejected, recover the full \(B\) from third‑order information and then inspect any implied lower‑triangular pattern. Given that the sample size in this exercise is still relatively moderate to estimate higher cumulants accurately, we instead focus on three testable implications that follow if \(UF\) is effectively “on top” of the contemporaneous causal order. The idea of the test mirrors DAVIS2023180, which effectively tests the joint plausibility of a given rotation combined with independent shocks.
Concretely, if the first row of \(B\) loads only on \(u_{UF}\) (we do not require full triangularity beyond this) and the shocks \(u_t\) satisfy A2$^\prime$–A3$^\prime$ and uncorrelatedness, then the \(2\times2\) subsystems \((e_{UF},e_{EPU})\) and \((e_{UF},e_{IP})\) can be written as linear combinations of shocks that satisfy joint diagonality of the third and second cumulants. To see this, \[
=
, \] so \(e_{UF}=u_{UF}\), \(e_{EPU}=b_{21}u_{UF}+u_{EPU}+b_{23}u_{IP}\), and \(e_{IP}=b_{31}u_{UF}+b_{32}u_{EPU}+u_{IP}\). Let \(u_{\mathrm{comb}1}:=u_{EPU}+b_{23}u_{IP}\) and \(u_{\mathrm{comb}2}:=b_{32}u_{EPU}+u_{IP}\). Then \[ T_{UF,EPU}
=
, \qquad T_{UF,EPU} :=
, \] and analogously \( T_{UF,IP}
=
\) with \( T_{UF,IP}=
. \) If \(u_t\) satisfies A2$^\prime$–A3$^\prime$ and uncorrelatedness, the pairs \((u_{UF},u_{\mathrm{comb}1})\) and \((u_{UF},u_{\mathrm{comb}2})\) retain diagonal third cumulants and diagonal covariance, so each \(2\times2\) subsystem meets the joint diagonality restriction under the null.
By contrast, the subsystem \((e_{EPU},e_{IP})\) generically mixes three shocks in a two‑dimensional space. Absent knife‑edge cancellations (e.g., \(b_{21}=b_{31}=0\) or vanishing third cumulants in just the right combination), \((e_{EPU},e_{IP})\) will violate joint diagonality of the second and third cumulants. This is why we can test the joint validity of the assumptions by focusing on these three two‑equation tests.
Applying this procedure, the over‑identification test rejects joint diagonality for \((e_{EPU},e_{IP})\) at the 5% level (Delta‑method \(p<0.01\)), while it does not reject for the pairs \((e_{UF},e_{EPU})\) and \((e_{UF},e_{IP})\) at the 5% level. On the surface, this pattern is what the “\(UF\) on top’’ hypothesis predicts: pairs that include \(UF\) behave like two‑source systems and pass, whereas \((EPU,IP)\) fails the minimal over‑identifying restriction. Looking more closely, however, \((e_{UF},e_{IP})\) is borderline at the 5% level (Delta‑method \(p=0.055\)), which is suggestive of modest tension with exact joint diagonality—consistent with a small contemporaneous channel between uncertainty and policy uncertainty or with specification noise. By contrast, the large \(p\)‑value for \((e_{UF},e_{EPU})\) (\(p=0.65\)) indicates no detectable deviation from the two‑source benchmark in that pair. Taken together, these results are broadly in line with https://doi.org/10.1002/jae.2672, in the sense that uncertainty shocks appear to help drive the cycle, while also indicating that \(UF\) is not perfectly exogenous in the EPU specification. Moreover, it shows that a strict triangular contemporaneous ordering may be too rigid for this model, as ludvigson2021uncertainty suggested.
I would like to express my sincere gratitude to my supervisors, Aureo de Paula and Andrei Zeleneev, for their continuous support, guidance, and encouragement throughout this work. I am also indebted to Tim Christensen, Ben Deaner, Raffaella Giacomini, Dennis Kristensen, and Daniel Lewis for their invaluable feedback and advice. I am grateful to Yanziyi Zhang and Chen-Wei Hsiang for many helpful discussions. Any remaining errors are solely my own responsibility.
The author declares that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.