Identifying and Testing Recursive vs. Interdependent Links in Simultaneous Equation Models via the SIRE Package

1 Introduction

As is well known, each equation in a simultaneous equation model (SEM) represents a specific link between a dependent (endogenous) variable and a set of other variables which play an explicative role for the former. These links can reflect either one-way relations between the dependent and their explicative variables or two-ways relations, ascribable to the presence of feedback loops operating either at a systematic or a stochastic level. SEMs are of recursive type as long as the equations represent unidirectional links. Otherwise, if the equations are bidirectional, the SEM (or part of it) is interdependent. Interdependence is, both structurally connected to the presence of current endogenous variables playing an explicative role, and can result as a by-product of error-term dependencies.

Investigating the nature, causal rather than interdependent, of a SEM is important in several respects. First the analysis, unfolding the dynamics among variables, sheds more light on the rationale behind the theoretical assumptions of the model. For instance, in an economic framework, the distinction between interdependent and causal SEMs leads to models which can be traced back to two main streams of economic theory: Neoclassical and Keynesian (Bellino et al. (2018), (2018)). Furthermore, the implication of interdependence vs. causality is crucial for undertaking parameter estimation, given that a set of causal equations can be estimated equation by equation by ordinary least squares (OLS), while simultaneous estimation methods, like three stage least squares (3SLS) are required when interdependence occurs. Given that large SEMs have become increasingly popular, the need for an analytical set-up, able to effectively detect and test causality versus interdependence, has of course become more urgent.

Starting from this premise and following Strotz and Wold (1960), (1960); Wold (1964), (1964); and more recently Faliva (1992), (1992); Faliva and Zoia (1994), (1994)); in this paper we have devised an operational method to distinguish the causal from the interdependent equations of a SEM.

Other approaches for detecting feedback-loops arising in deterministic (error free) models are based on either graph or system theory (see e.g., Gilli (1992,)). Our methodological proposal goes beyond the aforementioned methods, as besides covering both the cases of deterministic and error-driven feedback effects, it provides a way for testing the feedback effectiveness. In addition, it differs in principle from other approaches, as the one proposed by Granger (see Granger (1980), (1980)) and the Covariance Structural Analysis (CSA; Jöreskog (1978)). The former essentially rests on a predictability criterion for defining causality regardless of the theory behind the model. The latter, which is meant to find the best parametric approximation of the sample covariance matrix in terms of a given theoretical SEM structure; as such, it does not lead to a causal/interdependent interpretation of the model links as the one developed in our paper.

The feedbacks identified by the method proposed here demand statistical confirmation on certain empirical evidence arguments. Lack of significance of (one or more of) the estimated feedbacks can overturn the nature of the connections among model variables. To this end, a Wald type test is devised to check whether a given equation is significantly affected by feedback or not. The statistic of this test hinges on the parameter matrices of the model: the matrix associated to the endogenous variables playing an explicative role and the dispersion matrix of the error terms. If an equation is affected by feedback loops, the testing procedure allows to diagnose which endogenous variables are significantly connected in the loop of interest. Indeed, testing the significance of feedbacks means also checking if the links among variables, suggested by the theory at the basis of the model, are confirmed according to an empirical evidence argument.

The methodological approach put forth in this paper is implemented in R with the SIRE package. Besides integrating functions usually employed for the estimation of SEM’s, the package provides new functions meant to duly split a system of equations into its unidirectional and bidirectional links, and test their significance. To our knowledge, extant alternative approaches to causality do not offer a similar test.

The paper is structured as follows. The first section provides the methodological set-up devised to single out causal and interdependent relations in a SEM. In the second section, a Wald-type test is worked out to check whether a given equation is affected by feedbacks or not. The third section shows how the method and the R code work for detecting and testing feedback-loops in a macroeconomic model. An Appendix, with proofs of the main theoretical results, completes the paper.

2 Detecting Loops in an Equation System

An equation system is a set of structural equations representing economic theory-driven relations linking the variables relevant to the study at hand.

It is customary to specify an equation system as follows $$\begin{array}{}\text{(1)}& {\mathbf{y}}_t=\mathbf{\Gamma }{\mathbf{y}}_t+\mathbf{A}{\mathbf{z}}_t+{\mathit{ϵ}}_{t}t=1,\dots ,T\end{array}$$ where ${\mathbf{y}}_t$ is a $L\times 1$ vector of current dependent or endogenous variables, ${\mathbf{z}}_t$ is a $J\times 1$ vector of explicative variables and ${\mathit{ϵ}}_t$ is a $L\times 1$ vector of error terms. $T$ is the sample period. $\mathbf{\Gamma }$ and $\mathbf{A}$ are, respectively, $L\times L$ and $L\times J$ sparse parameter matrices. In particular $\mathbf{\Gamma }$, expressing the relations among current endogenous variables, is a hollow matrix to prevent any endogenous variable from explaining itself. Furthermore, it is assumed that ($\mathbf{I}-\mathbf{\Gamma }$) is of full rank, meaning that the equations are linearly independent.

Error terms are assumed to be non-systematic, stationary in a wide sense, and uncorrelated over time, that is $$\begin{array}{}\text{(2)}& \begin{array}{rl}E({\mathit{ϵ}}_t)& ={\mathbf{0}}_L\\ E({\mathit{ϵ}}_t{\mathit{ϵ}}_{\tau }^{\prime })& =\{\begin{array}{ll}{\mathbf{\Sigma }}_{(L\times L)}& \text{if}t=\tau \\ {\mathbf{0}}_{(L\times L)}& \text{if}t\ne \tau \end{array}\end{array}\end{array}$$ Actually, the pattern of relations recognizable in an econometric model can be interpreted either in terms of causal or interdependent schemes. A causal relation among variables is an asymmetric, theoretically-grounded and predictive relations which can be ideally meant as a stimulus-response mechanism (see Wold (1964), (1964) and Strotz and Wold (1960,)). The equations of a model form a causal chain when, once they are properly ordered, each current endogenous variable turns out to be, on the one hand, resultant of the joint effect of the endogenous which precede it in the chain and, on the other hand, cause of the current endogenous which follow the same endogenous in the chain. A model with equations that form a causal chain is defined recursive. The following simple equation system provides an example of a recursive model (see Figure 1, left panel) $$\begin{array}{}\text{(3)}& \begin{array}{rl}{y}_{1,t}=& {\mathit{a}}_1^{\prime }{\mathit{z}}_t+{ϵ}_{1,t}\\ y_{2,t}=& {\gamma }_{2,1}{y}_{1,t}+{\mathit{a}}_2^{\prime }{\mathit{z}}_t+{ϵ}_{2,t}\\ y_{3,t}=& {\gamma }_{3,2}{y}_{2,t}+{\gamma }_{3,1}{y}_{1,t}+{\mathit{a}}_3^{\prime }{\mathit{z}}_t+{ϵ}_{3,t}\\ y_{4,t}=& {\gamma }_{4,3}{y}_{3,t}+{\gamma }_{4,1}{y}_{1,t}+{\mathit{a}}_4^{\prime }{\mathit{z}}_t+{ϵ}_{4,t}\end{array}\end{array}$$ Recursive systems can be easily estimated, equation by equation, using OLS, starting from the top of the chain.

When a causal chain exists among blocks of current endogenous variables, a causal order can be established among those blocks of equations. In this case, the current endogenous variables of a block are effects of the variables belonging to the blocks which come before them in the chain, as well as the causes of the variables belonging to blocks which follow the block at stake in the chain. In this case, the model is of block-recursive type. The following simple equation system provides an example of a recursive model (see Figure 1, middle panel) $$\begin{array}{}\text{(4)}& \begin{array}{rl}{y}_{1,t}=& {\gamma }_{1,2}{y}_{2,t}+{\mathit{a}}_1^{\prime }{\mathit{z}}_t+{ϵ}_{1,t}\\ y_{2,t}=& {\gamma }_{2,1}{y}_{1,t}+{\mathit{a}}_2^{\prime }{\mathit{z}}_t+{ϵ}_{2,t}\\ y_{3,t}=& {\gamma }_{3,2}{y}_{2,t}+{\gamma }_{3,4}{y}_{4,t}+{\mathit{a}}_3^{\prime }{\mathit{z}}_t+{ϵ}_{3,t}\\ y_{4,t}=& {\gamma }_{4,3}{y}_{3,t}+{\gamma }_{4,1}{y}_{1,t}+{\mathit{a}}_4^{\prime }{\mathit{z}}_t+{ϵ}_{4,t}\end{array}\end{array}$$ Here, the chain is formed by two blocks of variables ($y_1$, $y_2$) and ($y_3$ and $y_4$) with the variables of the first block explaining those of the second.

Sometimes the composite nature of the connections among variables leads to a closed sequence of dependencies among variables to be ascribed to feedback loops. This type of interaction among endogenous variables is usually called interdependence. Interdependence is structurally connected to the presence of both current endogenous variables on the right-hand side of the model and the correlation between contemporaneous error terms.See the system below as an example in this regard (see Figure 1, right panel) $$\begin{array}{}\text{(5)}& \begin{array}{rl}{y}_{1,t}=& {\gamma }_{1,2}{y}_{2,t}+{\mathit{a}}_1^{\prime }{\mathit{z}}_t+{ϵ}_{1,t}\\ y_{2,t}=& {\gamma }_{2,1}{y}_{1,t}+{\gamma }_{2,3}{y}_{3,t}+{\mathit{a}}_2^{\prime }{\mathit{z}}_t+{ϵ}_{2,t}\\ y_{3,t}=& {\gamma }_{3,2}{y}_{2,t}+{\gamma }_{3,4}{y}_{4,t}+{\mathit{a}}_3^{\prime }{\mathit{z}}_t+{ϵ}_{3,t}\\ y_{4,t}=& {\gamma }_{4,3}{y}_{3,t}+{\gamma }_{4,1}{y}_{1,t}+{\mathit{a}}_4^{\prime }{\mathit{z}}_t+{ϵ}_{4,t}\end{array}\end{array}$$

Recursive model .	block-recursive model .	interdependent model .

Figure 1: The three patterns of relations in a simultaneous equation model.

Based on this premise, it is clear that the causal or interdependent features of a model’s equations depend on the pair of matrices $\mathbf{\Gamma }$ and $\mathbf{\Sigma }$. The former matrix highlights the possible (circular) dependencies or feedbacks among endogenous variables, while the latter features those induced by the stochastic components. In fact, the correlation of error terms associated to an equation-pair may transform the link between the endogenous, explained by these equations, into a relation with feedback.

Moreover, the essential information concerning the causal structure of a model can be obtained from the topological propertiesThe term topological properties refers to those properties of a matrix which depend exclusively on the number and the relative position of its null and non-null elements (Marimont (1969), (1969)). of the pair of the mentioned matrices and, at the very end, from the topological properties of the associated binary matrices ${\mathbf{\Gamma }}^b$and ${\mathbf{\Sigma }}^b$. A binary matrix associated to a matrix $\mathbf{G}$ is a matrix whose entries are equal to 1 if the corresponding entries of $\mathbf{G}$ are non-null, or 0 otherwise. Binary matrices preserve the topological properties of the parent matrices.

Following Faliva (Faliva (1992), (1992)) matrix $\mathbf{\Gamma }$ can be split as follows $$\begin{array}{}\text{(6)}& \mathbf{\Gamma }=\stackrel{\mathbf{~}}{\mathbf{C}}+{\mathbf{\Psi }}_0\end{array}$$ where $\stackrel{\mathbf{~}}{\mathbf{C}}$ includes the coefficients associated to current endogenous variables involved in feedback loops, and ${\mathbf{\Psi }}_0$ those associated to endogenous variables involved in causal relations.

Matrix $\stackrel{\mathbf{~}}{\mathbf{C}}$ is specified as follows $$\begin{array}{}\text{(7)}& \stackrel{\mathbf{~}}{\mathbf{C}}=\mathbf{C}+{\mathbf{\Psi }}_1\end{array}$$ where $\mathbf{C}$ includes the feedbacks arising in the systematic part of the model and matrix ${\mathbf{\Psi }}_1$ those induced by the correlation of the error terms. Matrices $\mathbf{C}$ and ${\mathbf{\Psi }}_1$ are defined as follows $$\begin{array}{}\text{(8)}& \begin{array}{r}\mathbf{C}=\mathbf{\Gamma }\ast \mathbf{R}\mathbf{R}={\left\{{\left[\sum _{r=1}^{L-1}{\left({\mathbf{\Gamma }}^b\right)}^r\right]}^b\right\}}^{\prime }\end{array}\end{array}$$

$$\begin{array}{}\text{(9)}& {\mathbf{\Psi }}_1=(\mathbf{\Gamma }-\mathbf{C})\ast {\left[{\mathbf{\Sigma }}^b(\mathbf{I}+\mathbf{R})\right]}^b,\end{array}$$ where the symbol "$\ast$" denotes the Hadamard product.The Hadamard product of two matrices, $\mathbf{A}$ and $\mathbf{B}$ of the same order, is defined as the matrix of the term-to-term products of the elements of these matrices, that is $(\mathbf{A}\ast \mathbf{B}{)}_{(i,j)}=a_{(i,j)}{b}_{(i,j)}$.An alternative approach for determining the feedbacks operating at a systematic level in a model is based on graph theory (see Jöreskog and Wold (1982), (1982) and Ponstein (1966), (1966)).

The rationale of (8) hinges on the fact that a direct feedback between variables $y_i$ and $y_j$ corresponds to the simultaneous non-nullity of ${\gamma }_{i,j}$ and ${\gamma }_{j,i}$ of coefficient matrix $\mathbf{\Gamma }$. This entails that a direct feedback between these two variables exists if the $(i,j)$-th element of the matrix The element ${\gamma }_{j,i}$ of $\mathbf{\Gamma }$ corresponds to the element ${\gamma }_{i,j}$ of ${\mathbf{\Gamma }}^{\prime }$ $$\begin{array}{}\text{(10)}& \mathbf{\Gamma }\ast ({\mathbf{\Gamma }}^b{)}^{\prime }\end{array}$$ is non null. An indirect feedback between the same variables is instead associated to a bidirectional connection between $y_i$ and $y_j$ established

through other variables and equations. In algebraic terms this corresponds to the simultaneous non-nullity of the $(i,j)$-th element of $\mathbf{\Gamma }$ and of the $(i,j)$-th element of a positive power of ${\mathbf{\Gamma }}^{\prime }$ (Fiedler (2013), (2013)). This entails that an indirect feedback exists between the mentioned variables if the $(i,j)$-th element of the following matrix $$\begin{array}{}\text{(11)}& {\mathbf{\Gamma }}^{\prime }\ast {\left\{{\left[\sum _{r=2}^{L-1}{\left({\mathbf{\Gamma }}^b\right)}^r\right]}^b\right\}}^{\prime }\end{array}$$ is non-null.

Accordingly, matrix $$\begin{array}{}\text{(12)}& \mathbf{\Psi }=\mathbf{\Gamma }-\mathbf{C}\end{array}$$ includes the coefficients associated to endogenous variables which, as far as the systematic aspects of the model are concerned, have a causal role.It is worth mentioning that $\mathbf{\Psi }$ is Hadamard-orthogonal to $\mathbf{C}$ (two matrices $\mathbf{A}$ and $\mathbf{B}$ are said to be Hadamard-orthogonal if $\mathbf{A}\ast \mathbf{B}$=$\mathbf{0}$). Furthermore, while matrix $\mathbf{C}$ is co-spectral to $\mathbf{\Gamma }$ (i.e., they have the same eigenvalues), matrix $\mathbf{\Psi }$ is a hollow-nilpotent matrix, like $\mathbf{\Gamma }$ (a square matrix $\mathit{N}$ is nilpotent if ${\mathit{N}}^k=\mathbf{0}$ for some $k<M$, where M is the matrix dimension). A hollow, nilpotent matrix can always be expressed in triangular form.

In order to show how feedbacks operating in the systematic part of a model can be detected, let us consider as an example the following deterministic model $$\begin{array}{}\text{(13)}& \begin{array}{rl}{y}_{1,t}=& {\gamma }_{1,5}{y}_{5,t}+{\gamma }_{1,7}{y}_{7,t}+{\mathit{a}}_1^{\prime }{\mathit{z}}_t\\ y_{2,t}=& {\mathit{a}}_2^{\prime }{\mathit{z}}_{2,t}\\ y_{3,t}=& {\gamma }_{3,11}{y}_{11,t}+{\mathit{a}}_3^{\prime }{\mathit{z}}_t\\ y_{4,t}=& {\gamma }_{4,3}{y}_{3,t}+{\mathit{a}}_4^{\prime }{\mathit{z}}_t\\ y_{5,t}=& {\gamma }_{5,10}{y}_{10,t}+{\mathit{a}}_5^{\prime }{\mathit{z}}_t\\ y_{6,t}=& {\gamma }_{6,5}{y}_{5,t}+{\gamma }_{6,9}{y}_{9,t}\mathbf{+}{a}_6^{\prime }{\mathit{z}}_t\\ y_{7,t}=& {\gamma }_{7,6}{y}_{6,t}+{\mathit{a}}_7^{\prime }{\mathit{z}}_t\\ y_{8,t}=& {\gamma }_{8,12}{y}_{12,t}+{\mathit{a}}_8^{\prime }{\mathit{z}}_t\\ y_{9,t}=& {\gamma }_{9,7}{y}_{7,t}+{\mathit{a}}_9^{\prime }{\mathit{z}}_t\\ y_{10,t}=& {\gamma }_{10,5}{y}_{5,t}+{\mathit{a}}_{10}^{\prime }{\mathit{z}}_{2,t}\\ y_{11,t}=& {\gamma }_{11,12}{y}_{12,t}+{\mathit{a}}_{11}^{\prime }{\mathit{z}}_t\\ y_{12,t}=& {\gamma }_{12,4}{y}_{4,t}+{\gamma }_{12,11}{y}_{11,t}+{\mathit{a}}_{12}^{\prime }{\mathit{z}}_t\\ y_{13,t}=& {\gamma }_{13,2}{y}_{2,t}+{\gamma }_{13,6}{y}_{6,t}+{\mathit{a}}_{13}^{\prime }{\mathit{z}}_t\end{array}\end{array}$$ Matrix ${\mathbf{\Gamma }}^b$ is given by $${\mathbf{\Gamma }}^b=\left[\begin{array}{ccccccccccccc}\cdot & \cdot & \cdot & \cdot & 1& \cdot & 1& \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 1& \cdot & \cdot \\ \cdot & \cdot & 1& \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 1& \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & 1& \cdot & \cdot & \cdot & 1& \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & 1& \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 1& \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 1& \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & 1& \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 1& \cdot \\ \cdot & \cdot & \cdot & 1& \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 1& \cdot & \cdot \\ \cdot & 1& \cdot & \cdot & \cdot & 1& \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \end{array}\right]$$ Using (8) and (12), ${\mathbf{\Gamma }}^b$ is split in the following two submatrices $$\begin{array}{}\text{(14)}& {\mathbf{C}}^b=\left[\begin{array}{ccccccccccccc}\cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 1& \cdot & \cdot \\ \cdot & \cdot & 1& \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 1& \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 1& \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & 1& \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 1& \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & 1& \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 1& \cdot \\ \cdot & \cdot & \cdot & 1& \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 1& \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \end{array}\right],{\mathbf{\Psi }}^b=\left[\begin{array}{ccccccccccccc}\cdot & \cdot & \cdot & \cdot & 1& \cdot & 1& \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & 1& \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 1& \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & 1& \cdot & \cdot & \cdot & 1& \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \end{array}\right]\end{array}$$
Looking at matrix ${\mathbf{C}}^b$, we see that the simultaneous non-nullity of the $c_{5,10}$, $c_{10,5}$, $c_{11,12}$, and $c_{12,11}$ elements imply the existence of two direct feedbacks: one between the variable-pair $y_5$ and $y_{10}$, and the other between $y_{11}$ and $y_{12}$. The non-nullity of the $c_{3,11}$, $c_{4,3}$, and $c_{12,4}$ elements denotes the existence of indirect feedbacks between the four variables $y_3$, $y_4$, $y_{11}$, and $y_{12}$. Similarly, variables $y_6$, $y_7$, and $y_9$ are connected by an (indirect) feeback as a consequence of the non-nullity of the $c_{6,9}$, $c_{7,6}$, and $c_{9,7}$ elements. Looking at matrix $\mathbf{\Psi }$ we conclude that variables $y_5$ and $y_7$ have a causal role in the first equation. Variables $y_5$ and $y_{12}$ have the same role in the equations six and eight, while variables $y_2$ and $y_6$ play a causal role in the last equation. The results ensuing from the decomposition of ${\mathbf{\Gamma }}^b$ are depicted in Figure 2.

Figure 2: Interdependent links (in red) and causal links (in black) operating in the model (13).

If the error terms are correlated, the causal structure of a model could no longer match that of its systematic counterpart since part of the relations that are recursive at systematic level, namely $\mathbf{\Psi }{\mathbf{y}}_t$, may become interdependent as a consequence of the feedback mechanisms induced by the stochastic terms in $\mathbf{\Sigma }$. In this case, matrix $\mathbf{\Psi }$ turns out to be the sum of two Hadamard-orthogonal matrices, ${\mathbf{\Psi }}_0$ and ${\mathbf{\Psi }}_1$, that is $$\begin{array}{}\text{(15)}& \mathbf{\Psi }={\mathbf{\Psi }}_0+{\mathbf{\Psi }}_1{\mathbf{\Psi }}_0\ast {\mathbf{\Psi }}_1={\mathbf{0}}_{(L\times L)}\end{array}$$ where $$\begin{array}{}\text{(16)}& {\mathbf{\Psi }}_1^b=\mathbf{\Psi }\ast \mathbf{F}\mathbf{F}={\left[{\mathbf{\Sigma }}^b(\mathbf{I}+\mathbf{R})\right]}^b\end{array}$$ Here, matrix ${\mathbf{\Psi }}_1$ includes the coefficients associated to the endogenous variables involved in loops induced by disturbances. In fact, it can be proved (see 1. in Appendix) that the matrix ${\left[{\mathbf{\Sigma }}^b(\mathbf{I}+\mathbf{R})\right]}^b$ is the binary counterpart of the covariance matrix between the error terms and the endogenous variables given by $$\begin{array}{}\text{(17)}& \mathbb{E}({\mathit{ϵ}}_t^{\prime }{\mathbf{y}}_t)=[\mathbf{\Sigma }{(\mathbf{I}-\mathbf{\Gamma })}^{-1}]\end{array}$$ The non-null elements of the above matrix express the effect of the model’s left-hand side (LHS) endogenous variables on the right-hand side (RHS) ones, which are induced by the error term correlation.

Equations (15) and (16) rest on the logical relations between the concepts of causality and predictability, where the notion of optimal predictor (in mean-square sense) tallies with that of conditional expectation. In fact, given that causal relations are also predictive, but not vice-versa, we can define as causal those relations that are both causal in the deterministic model and predictive in a stochastic context. This means that if the conditional expectations of the relations, which are causal in the deterministic model, namely $\mathbf{\Psi }{\mathbf{y}}_t$, are not affected by the error terms, then $\mathbf{\Psi }{\mathbf{y}}_t$ turns out to also have a causal role in a stochastic context.

Accordingly, we can say that the stochastic specification is neutral with respect to the underlying systematic causal structure if the following holds (Faliva (1992), (1992)) $$\begin{array}{}\text{(18)}& \mathbb{E}(\mathbf{\Psi }{\mathbf{y}}_t+{\mathit{ϵ}}_t|\mathbf{\Psi }{\mathbf{y}}_t)=\mathbf{\Psi }{\mathbf{y}}_t+\mathbb{E}({\mathit{ϵ}}_t|\mathbf{\Psi }{\mathbf{y}}_t)=\mathbf{\Psi }{\mathbf{y}}_t\end{array}$$

meaning that $$\begin{array}{}\text{(19)}& \mathbb{E}({\mathit{ϵ}}_t|\mathbf{\Psi }{\mathbf{y}}_t)=\mathbf{0}\end{array}$$

Otherwise, the correlation between the error terms and the endogenous variables may affect the conditional expectation of the error term as follows (see Faliva (1992), (1992)) $$\begin{array}{}\text{(20)}& \mathbb{E}({\mathit{ϵ}}_t|\mathbf{\Psi }{\mathbf{y}}_t)=-{\mathbf{\Psi }}_1{\mathbf{y}}_t\end{array}$$ which, in turn, implies that $$\begin{array}{}\text{(21)}& \mathbb{E}(\mathbf{\Psi }{\mathbf{y}}_t+{\mathit{ϵ}}_t|\mathbf{\Psi }{\mathbf{y}}_t)=\mathbf{\Psi }{\mathbf{y}}_t-{\mathbf{\Psi }}_1{\mathbf{y}}_t={\mathbf{\Psi }}_0{\mathbf{y}}_t\end{array}$$ In this case, only the subset ${\mathbf{\Psi }}_0{\mathbf{y}}_t$ of the original set of causal relations, playing a predictive role, is causal. This, in turn, implies that the overall feedback operating in the system is included in matrix $\stackrel{\mathbf{~}}{\mathbf{C}}=\mathbf{C}+{\mathbf{\Psi }}_1$.

To highlight the role played by the stochastic specification on the model causal structure, let us consider as an example the following specification for matrix ${\mathbf{\Sigma }}^b$ $$\begin{array}{}\text{(22)}& {\mathbf{\Sigma }}^b=\left[\begin{array}{cc}1& \\ \cdot & 1& \\ \cdot & \cdot & 1& \\ 1& 1& \cdot & 1& \\ 1& \cdot & \cdot & 1& 1& \\ \cdot & 1& 1& 1& \cdot & 1& \\ \cdot & \cdot & 1& \cdot & \cdot & 1& 1& \\ 1& 1& \cdot & 1& 1& 1& \cdot & 1& \\ \cdot & 1& \cdot & 1& \cdot & 1& \cdot & 1& 1& \\ 1& \cdot & \cdot & 1& 1& \cdot & \cdot & 1& \cdot & 1& \\ \cdot & \cdot & 1& \cdot & \cdot & 1& 1& \cdot & \cdot & \cdot & 1& \\ 1& \cdot & \cdot & 1& 1& \cdot & \cdot & 1& \cdot & 1& \cdot & 1& \\ \cdot & \cdot & 1& \cdot & \cdot & 1& 1& \cdot & \cdot & \cdot & 1& \cdot & 1& \end{array}\right]\end{array}$$ Then, matrices $\widetilde{{\mathbf{C}}^b}$ and ${\mathbf{\Psi }}_0^b$ are $$\begin{array}{}\text{(23)}& \widetilde{{\mathbf{C}}^b}=\left[\begin{array}{ccccccccccccc}\cdot & \cdot & \cdot & \cdot & 1& \cdot & 1& \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 1& \cdot & \cdot \\ \cdot & \cdot & 1& \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 1& \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 1& \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & 1& \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 1& \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 1& \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & 1& \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 1& \cdot \\ \cdot & \cdot & \cdot & 1& \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 1& \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & 1& \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \end{array}\right]={\mathbf{C}}^b+{\mathbf{\Psi }}_1^b=\left[\begin{array}{ccccccccccccc}\cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 1& \cdot & \cdot \\ \cdot & \cdot & 1& \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 1& \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 1& \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & 1& \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 1& \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & 1& \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 1& \cdot \\ \cdot & \cdot & \cdot & 1& \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 1& \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \end{array}\right]+\left[\begin{array}{ccccccccccccc}\cdot & \cdot & \cdot & \cdot & 1& \cdot & 1& \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & 1& \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & 1& \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \end{array}\right]\end{array}$$

$$\begin{array}{}\text{(24)}& {\mathbf{\Psi }}_0^b=\left[\begin{array}{ccccccccccccc}\cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & 1& \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \\ \cdot & 1& \cdot & \cdot & \cdot & 1& \cdot & \cdot & \cdot & \cdot & \cdot & \cdot & \cdot \end{array}\right]\end{array}$$ The non-null correlation between the pairs of error terms $\{{ϵ}_5,{ϵ}_1\}$, $\{{ϵ}_6,{ϵ}_{13}\}$, $\{{ϵ}_7,{ϵ}_1\}$ and $\{{ϵ}_{12},{ϵ}_8\}$ (see Equation (22)) has transformed the relations among the pairs of variables $\{y_5,y_1\}$, $\{y_6,y_{13}\}$, $\{y_7,y_1\}$, and $\{y_{12},y_8\}$, which were causal in the deterministic model (13), into interdependent links. Figure 3 shows the effect of the stochastic specification (22) on the feedbacks originally detected in the deterministic model (13).

Figure 3: Interdependent (in red) and causal (in black) links operating in the model (13) when the stochastic specification is as in (22). Dashed red lines with double-headed arrows denote interdependent links induced by the correlation of the error terms.

The flow-chart in Figure 4 shows the different cases, according to the structure of matrices $\mathbf{\Gamma }$ and $\mathbf{\Sigma }$.

Figure 4: Flow-chart showing the possible outcome of the system decomposition in terms of $\mathbf{\Gamma }$ and $\mathbf{\Sigma }$.

3 Testing the Significance of Feedback Loops

In the previous section an analytic framework was set up to describe the potential feedbacks operating in a model. In fact, the analysis developed, relying on binary matrices, was meant to be qualitative since it only highlights the feedback set that potentially operates in a model, given the characteristics of its relations and its stochastic specification. Only once the model has been duly estimated, can the coefficients of matrix $\stackrel{\mathbf{~}}{\mathit{C}}$ be properly evaluated. At this point, it proves useful to devise a procedure for testing the significance of the estimated loops (see Faliva and Zoia (1994), (1994)). To this end, let us observe that, once the matrix including all the feedbacks operating in the model $$\begin{array}{}\text{(25)}& \begin{array}{r}\mathbf{C}+{\mathbf{\Psi }}_1=\mathbf{C}+{\mathbf{\Psi }}_1=\mathbf{C}+(\mathbf{\Psi }\ast \mathbf{F})\mathbf{F}=\mathbf{\Sigma }(\mathbf{I}-\mathbf{\Gamma }{)}^{-1}\end{array}\end{array}$$ have been properly estimated, a test for the effective functioning of feedback loops can be established, based on the significance of its non-null entries. Any given equation, say the $j$-th one, turns out to be involved in feedback loops with other equations of the model whenever the $j-$th row of the above matrix is not a null vector. Should the $(j,i)$-th entry of this matrix be non-null, then a feedback between the $j$-th and the $i$-th equation would be expected to exist (see (55) in the Appendix). Actually, it can be proved (see 2. in Appendix) that, in light of the identity

$$\begin{array}{}\text{(26)}& \begin{array}{r}\mathbf{C}+(\mathbf{\Psi }\ast \mathbf{F})=(\mathbf{C}+\mathbf{\Psi })\ast \mathbf{F}=\mathbf{\Gamma }\ast \mathbf{F}\end{array}\end{array}$$ a test for the significance of the loops can be based on the exam of the statistical non-nullity of the elements of matrix $\mathbf{\Gamma }\ast \mathbf{F}$ which, unlike $\stackrel{\mathbf{~}}{\mathbf{C}}$, does not require the preliminary split of $\mathbf{\Gamma }$ into its components, given the feedback loops $\mathbf{C}+{\mathbf{\Psi }}_1$ and causal links ${\mathbf{\Psi }}_0$.

In this context (following Faliva and Zoia (1994), (1994)), it can be proved that the $j$-th row of matrix $\mathbf{\Gamma }\ast \mathbf{F}$ measures both the direct effect of the RHS endogenous variables on the $j$-th one and the feedback effect of the latter on the former variables. In fact, the direct effects of the RHS endogenous variables, collected in vector ${\mathbf{y}}_o$, on variable $y_j$ are included in the $j$-th row of matrix $\mathbf{\Gamma }$ (excluding its $j$-th element), that is $$\begin{array}{}\text{(27)}& \frac{\partial E(y_j|{\mathbf{y}}_o)}{\partial {\mathbf{y}}_o}={\mathbf{e}}_j^{\prime }\mathbf{\Gamma }{\mathbf{M}}_j\end{array}$$ Here, ${\mathbf{e}}_j$ is the $L$-dimensional $j$-th elementary vector and ${\mathbf{M}}_j$ is the $(L\times (L-1))$ selection matrix obtained from the identity matrix by deleting its $j$-th column, that is $$\begin{array}{}\text{(28)}& {\mathbf{M}}_j=\left[\begin{array}{ccccc}\underset{(L,1)}{{\mathbf{e}}_1},& \dots & \underset{(L,1)}{{\mathbf{e}}_{j-1}},& \underset{(L,1)}{{\mathbf{e}}_{j+1}},& \dots \underset{(L,1)}{{\mathbf{e}}_{L-1}}\end{array}\right]\end{array}$$

The feedback effects of the $y_j$ variable on its explicative endogenous variables, ${\mathbf{y}}_o$, are included in the $j$-th row of matrix $\mathbf{F}$ (excluding its $j$-th element), that is $$\begin{array}{}\text{(29)}& \frac{\partial E({\mathbf{y}}_o^{\prime }|y_j)}{\partial y_j}=({\mathbf{M}}_j^{\prime }{\mathbf{F}}^{\prime }{\mathbf{e}}_j{)}^{\prime }\end{array}$$ To prove (29), let us focus on the $j$-th equation and consider this equation as the first of the system, with the others in sequence, that is $$\begin{array}{}\text{(30)}& \begin{array}{rl}\underset{(1,1)}{y_j}& =\underset{(1,L-1)}{{\mathit{\gamma }}_j^{\prime }}\underset{(L-1,1)}{{\mathbf{y}}_o}+\underset{(1,J)}{{\mathbf{a}}_j^{\prime }}\underset{(J,1)}{\mathbf{z}}+\underset{(1,1)}{{ϵ}_j}\end{array}\end{array}$$

$$\begin{array}{}\text{(31)}& \begin{array}{rl}\underset{(L-1,1)}{{\mathbf{y}}_o}& =\underset{(L-1,1)}{\mathit{\eta }}\underset{(1,1)}{y_j}+\underset{(L-1,L-1)}{{\mathbf{\Gamma }}_o}\underset{(L-1,1)}{{\mathbf{y}}_o}+\underset{(L-1,J)}{{\mathbf{A}}_o}\underset{(J-1,1)}{\mathbf{z}}+\underset{(L-1,1)}{{\mathit{ϵ}}_o}\end{array}\end{array}$$

$$\begin{array}{}\text{(32)}& \begin{array}{rl}\left(\begin{array}{c}{ϵ}_j\\ {\mathit{ϵ}}_o\end{array}\right)& \sim {\mathcal{N}}_L(\mathbf{0},\mathbf{\Sigma })\text{where}\mathbf{\Sigma }=\left(\begin{array}{cc}{\sigma }_{jj}& {\mathit{\sigma }}_{jo}\\ {\mathit{\sigma }}_{oj}& {\mathbf{\Sigma }}_o\end{array}\right)\end{array}\end{array}$$ Looking at the $j$-th equation, it is clear that vector ${\mathit{\gamma }}_j^{\prime }={\mathbf{e}}_j^{\prime }\mathbf{\Gamma }{\mathbf{M}}_j$ measures the direct effect of the (RHS) endogenous variables on $y_j$.

In order to determine the feedback effect of $y_j$ on ${\mathbf{y}}_o$, let us rewrite (31) as follows $$\begin{array}{}\text{(33)}& \begin{array}{rl}{\mathbf{y}}_o& =\mathit{\eta }({\mathit{\gamma }}_j^{\prime }{\mathbf{y}}_o+{\mathbf{a}}_j^{\prime }\mathbf{z}+{ϵ}_j)+{\mathbf{\Gamma }}_o{\mathbf{y}}_o+{\mathbf{A}}_o\mathbf{z}+{\mathit{ϵ}}_o\end{array}\end{array}$$ Next, given that, under normality, the following holds $$\begin{array}{}\text{(34)}& \begin{array}{r}{\mathit{ϵ}}_o=\frac{{\mathit{\sigma }}_{oj}}{{\sigma }_{jj}}{ϵ}_j+{\mathit{\zeta }}_o;{\mathit{\zeta }}_o\mathrm{\perp }{ϵ}_j\end{array}\end{array}$$ the set of equations ((33)) can be conveniently rewritten in the form $$\begin{array}{}\text{(35)}& \begin{array}{rl}(\mathbf{I}-\mathbf{G}){\mathbf{y}}_o& =\mathbf{D}\mathbf{z}+\mathbf{d}{\mathit{ϵ}}_j+{\mathit{\zeta }}_o\end{array}\end{array}$$ where $$\begin{array}{}\text{(36)}& \begin{array}{r}\mathbf{G}=\mathit{\eta }{\mathit{\gamma }}_j^{\prime }+{\mathbf{\Gamma }}_o;\mathbf{D}=\mathit{\eta }{\mathbf{a}}_j^{\prime }+{\mathbf{A}}_o;\mathbf{d}=\mathit{\eta }+\frac{{\mathit{\sigma }}_{oj}}{{\sigma }_{jj}}\end{array}\end{array}$$ This, in turn (see 3. in Appendix) entails $$\begin{array}{}\text{(37)}& \frac{\partial E({\mathbf{y}}_o^{\prime }|{ϵ}_j)}{\partial y_j}=\frac{\partial E({\mathbf{y}}_o^{\prime }|{ϵ}_j)}{\partial {ϵ}_j}\frac{\partial {ϵ}_j}{\partial y_j}=[(\mathbf{I}-\mathbf{G}{)}^{-1}\mathbf{d}{]}^{\prime }={\mathit{\phi }}_j^{\prime }=\frac{1}{{\sigma }_{jj}}{\mathbf{e}}_j^{\prime }\mathbf{F}{\mathbf{M}}_j\end{array}$$ Thus, we can conclude that the presence of non-zero elements in the vector $$\begin{array}{}\text{(38)}& \begin{array}{r}{\mathit{\rho }}_j^{\prime }={\mathit{\gamma }}_j^{\prime }\ast {\mathit{\phi }}_j^{\prime }=({\mathbf{e}}_j^{\prime }\mathbf{\Gamma }{\mathbf{M}}_j)\ast (\frac{1}{{\sigma }_{jj}}{\mathbf{e}}_j^{\prime }\mathbf{F}{\mathbf{M}}_j)={\mathbf{e}}_j^{\prime }\frac{1}{{\sigma }_{jj}}(\mathbf{\Gamma }\ast \mathbf{F}){\mathbf{M}}_j\end{array}\end{array}$$ reveals the simultaneous action of both the direct effects of ${\mathbf{y}}_o$ on $y_j$ and the feedback effects of $y_j$ on ${\mathbf{y}}_o$.

Accordingly, testing the significance of ${\mathit{\rho }}_j$ means checking whether the $j$-th endogenous is involved in feedback loops with other endogenous variables.

Actually, the statistic of the test can be derived from (38), by deleting from ${\mathit{\gamma }}_j^{\prime }$ the elements that, according to the exclusion constraints postulated by the economic theory, are null. This leads to move from the former ${\mathit{\rho }}_j$ vector to the following compressed vector $$\begin{array}{}\text{(39)}& \begin{array}{r}{\stackrel{\mathbf{~}}{\mathit{\rho }}}_j={\stackrel{\mathbf{~}}{\mathit{\gamma }}}_j^{\prime }\ast {\widetilde{{\mathit{\phi }}_j}}^{\prime }=({\mathbf{S}}_j{\mathit{\gamma }}_j{)}^{\prime }\ast ({\mathbf{S}}_j{\mathit{\phi }}_j{)}^{\prime }\end{array}\end{array}$$ which has no zero entries. Here ${\mathbf{S}}_j$ is a selection matrix selecting from ${\mathit{\gamma }}_j$ and ${\mathit{\phi }}_j$ the non-null entries. Accordingly, the reference model (30)-(32) can be restated as $$\begin{array}{}\text{(40)}& \begin{array}{rl}{y}_j& ={\stackrel{\mathbf{~}}{\mathit{\gamma }}}_j^{\prime }{\mathbf{y}}_r+{\stackrel{\mathbf{~}}{\mathbf{a}}}_j^{\prime }{\mathbf{z}}_r+{ϵ}_j\end{array}\end{array}$$

$$\begin{array}{}\text{(41)}& \begin{array}{rl}{\mathbf{y}}_r& =\mathbf{K}\mathbf{z}+{\tilde{\mathit{\phi }}}_j{ϵ}_j+{\mathit{ϵ}}_r\end{array}\end{array}$$

$$\begin{array}{}\text{(42)}& \begin{array}{rl}{f}_{({ϵ}_j,{\mathit{ϵ}}_r)}& \sim {\mathcal{N}}_L(\mathbf{0},\left(\begin{array}{cc}{\sigma }_{jj}& {\mathbf{0}}^{\prime }\\ \mathbf{0}& \mathbf{\Omega }\end{array}\right))\end{array}\end{array}$$ where $${\mathbf{y}}_r={\mathbf{S}}_j{\mathbf{y}}_o,{\stackrel{\mathbf{~}}{\mathbf{a}}}_j={\mathbf{S}}_r{\mathbf{a}}_j,{\mathbf{z}}_r={\mathbf{S}}_r\mathbf{z}$$ and ${\mathbf{S}}_r$ is the matrix selecting the non-null entries from ${\mathbf{a}}_j^{\prime }$ and the sub-set of predetermined variables playing an explicative role in the $j$-th equation.

Furthermore, $$\begin{array}{}\text{(43)}& \begin{array}{r}\mathbf{K}={\mathbf{S}}_j(\mathbf{I}-\mathbf{G}{)}^{-1}\mathbf{D},{\mathit{ϵ}}_r={\mathbf{S}}_j(\mathbf{I}-\mathbf{G}{)}^{-1}{\mathit{\zeta }}_o,\mathbf{\Omega }=E({\mathit{ϵ}}_r{\mathit{ϵ}}_r^{\prime })\end{array}\end{array}$$ Hence, the issue of the nature, unidirectional rather than bidirectional, of the equation at stake can be unfolded by testing a hypothesis in the given form $$\begin{array}{}\text{(44)}& \begin{array}{r}\{\begin{array}{l}{H}_0:{\stackrel{\mathbf{~}}{\mathit{\rho }}}_j=\mathbf{0}\\ H_1:{\stackrel{\mathbf{~}}{\mathit{\rho }}}_j\ne \mathbf{0}\end{array}\end{array}\end{array}$$

The Wald test takes the form $$\begin{array}{}\text{(45)}& W={{\widehat{\tilde{\mathit{\rho }}}}_j}^{\prime }(\hat{\mathbf{J}}{\hat{\mathbf{\Psi }}}^{-1}{\hat{\mathbf{J}}}^{\prime }{)}^{-1}{\widehat{\tilde{\mathit{\rho }}}}_j\end{array}$$ where ${\widehat{\tilde{\mathit{\rho }}}}_j$ is the maximum likelihood estimate of ${\tilde{\mathit{\rho }}}_j$ (see 4. in Appendix), and $\hat{\mathbf{J}}$, $\hat{\mathbf{\Psi }}$ are, respectively, the Jacobian matrix $$\begin{array}{}\text{(46)}& \hat{\mathbf{J}}=\frac{\partial {\tilde{\mathit{\rho }}}_j(\mathit{\theta })}{\partial \mathit{\theta }}{|}_{\mathit{\theta }=\hat{\mathit{\theta }}}\end{array}$$ and the information matrix $$\begin{array}{}\text{(47)}& \hat{\mathbf{\Psi }}=\frac{{\partial }^{2}l(\mathit{\theta })}{\partial \mathit{\theta }\partial {\mathit{\theta }}^{\prime }}{|}_{\mathit{\theta }=\hat{\mathit{\theta }}}\end{array}$$ evaluated in correspondence of the maximum likelihood estimate of the parameter vector $$\begin{array}{}\text{(48)}& {\mathit{\theta }}^{\prime }=[{\stackrel{\mathbf{~}}{\mathit{\gamma }}}_j^{\prime },{\stackrel{\mathbf{~}}{\mathbf{a}}}_j^{\prime },{\stackrel{\mathbf{~}}{\mathit{\phi }}}_j^{\prime },\text{vec}(\mathbf{K}),{\sigma }_{jj},\text{vech}(\mathbf{\Omega }){]}^{\prime }\end{array}$$ Under the null hypothesis $$\begin{array}{}\text{(49)}& W\underset{as.}{\approx }{\chi }_k^2\end{array}$$ where $k$ is the dimension of ${\tilde{\mathit{\rho }}}_j$.

If the Wald test provides evidence that the $j$-th equation is involved in a statistically significant feedback loop with other equations of the model, it is worth singling out the variables that are primarily responsible for the feedback at hand. They can be identified by checking the significance of each non-null element of ${\widehat{\tilde{\mathit{\rho }}}}_j$. Under the null that the $i$-th element of ${\widehat{\tilde{\mathit{\rho }}}}_j$ is non-zero, the Wald statistic, for testing the significance of the loop bridging the $i$-th and $j$-th endogenous, turns out to be $$\begin{array}{}\text{(50)}& W=({\mathbf{e}}_i^{\prime }{\widehat{\tilde{\mathit{\rho }}}}_j{)}^{\prime }[{\mathbf{e}}_i(\hat{\mathbf{J}}{\hat{\mathbf{\Psi }}}^{-1}{\hat{\mathbf{J}}}^{\prime }{)}^{-1}{\mathbf{e}}_i^{\prime }]({\widehat{\tilde{\mathit{\rho }}}}_j{\mathbf{e}}_i)\underset{as.}{\approx }{\chi }_1^2\end{array}$$

4 Detecting and testing causal and interdependent links in a model with SIRE

Investigating potential feedbacks with SIRE

The analysis developed in the previous sections allows the identification of the potential feedbacks operating in a model. By assuming the stochastic specification of the model as known, the investigation can be carried out by using binary matrices ${\mathbf{\Gamma }}^b$ and ${\mathbf{\Sigma }}^b$ without a preliminary estimation of the model. The causal structure, which emerges from this analysis, is implied by the theory underlying the model and mirrored by the topological properties of matrices $\mathbf{\Gamma }$ and $\mathbf{\Sigma }$. It is also important to point out that the feedback loops thus detected are only potential, because their effectiveness must find confirmation in empirical evidence.

We start by loading the SIRE package.

> install.packages("SIRE")
> library(SIRE)

The function causal_decompose() is devised for decomposing the matrix $\mathbf{\Gamma }$. If the structure of $\mathbf{\Sigma }$ is assumed as known by the user, the function takes the following arguments:

• data: not appropriate to simulated context, set to NULL.
• eq.system: the system of equations.
• resid.est: not appropriate to simulated context, set to NULL.
• instruments: not appropriate to simulated context, set to NULL.
• sigma.in: the binary matrix ${\mathbf{\Sigma }}^b$.

and provides the following output:

• eq.system: the system of equations given as input.
• gamma: the binary matrix $\mathbf{\Gamma }{}^b$.
• sigma: the binary matrix ${\mathbf{\Sigma }}^b$ given as input.
• C: the binary matrix of the coefficients associated to the endogenous variables involved in interdependent mechanisms operating at a systematic level.
• Psi1: the binary matrix of the coefficients associated to the endogenous variables involved in interdependent mechanisms induced by error correlation (if Sigma is not diagonal).
• Psi0: the binary matrix of the coefficients associated to the endogenous variables having a causal role.
• all.graph: the DAG object for the undecomposed path diagram (via the R package igraph; Amestoy (2017), (2017)).
• dec.graph: the DAG object for the decomposed path diagram.

Furthermore, if the error terms are assumed to be spherical, then the SIRE package simply splits $\mathbf{\Gamma }$ in two sub-matrices $\mathbf{C}{}^b$ and $\mathbf{\Psi }{}^b$, reflecting the interdependent and causal relations operating in the system at a deterministic level.

With regard to the system (13), the corresponding code is

> eq.system <- list(
+              eq1 = y1 ~ y5 + y7, eq2 = y2 ~ z,
+              eq3 = y3 ~ y11, eq4 = y4 ~ y3,
+              eq5 = y5 ~ y10, eq6 = y6 ~ y5 + y9,
+              eq7 = y7 ~ y6, eq8 = y8 ~ y12,
+              eq9 = y9 ~ y7, eq10 = y10 ~ y5,
+              eq11 = y11 ~ y12, eq12 = y12 ~ y4 + y11,
+              eq13 = y13 ~ y2 + y6)
> #fictitious Sigma matrix
> Sigma <- diag(length(eq.system))
> #function call
> decompose.A <- causal_decompose(eq.system , sigma.in = Sigma)

The output is comprised of matrices $\mathbf{C}{}^b$ and $\mathbf{\Psi }{}^b$ given in (14). The graphical representation of the system, given in Figure 2, is obtained with the tkplot() function of the R package igraph

> tkplot(decompose.A$dec.graph)

The following example refers to a matrix ${\mathbf{\Sigma }}^b$ specified as in (22)

> # indexes of non-null elements of Sigma
> sigma.idx <- cbind(
+  rbind(rep(1,5),c(4,5,8,10,12)),   #y1
+  rbind(rep(2,4),c(4,6,8,9)),       #y2
+  rbind(rep(3,4),c(6,7,11,13)),     #y3
+  rbind(rep(4,6),c(5,6,8,9,10,12)), #y4
+  rbind(rep(5,3),c(8,10,12)),       #y5
+  rbind(rep(6,5),c(7,8,9,11,13)),   #y6
+  rbind(rep(7,2),c(11,13)),         #y7
+  rbind(rep(8,3),c(9,10,12)),       #y8
+  rbind(rep(10,1),c(12)),           #y10
+  rbind(rep(11,1),c(13)))           #y11
> # fictitious Sigma matrix
> low.tri <- as.matrix(Matrix::sparseMatrix(i = sigma.idx[2,] , j = sigma.idx[1,], x = 1,
+                                           dims = rep(length(eq.system),2)))
> Sigma <- low.tri + t(low.tri) + diag(length(eq.system))
> # function call
> decompose.B <- causal_decompose(eq.system = eq.system,
+                             sigma.in = Sigma)

In this case, the package provides as output matrix $\mathbf{C}{}^b$ and splits matrix $\mathbf{\Psi }{}^b$ into sub-matrices ${\mathbf{\Psi }}_1^b$ and ${\mathbf{\Psi }}_0^b$, as in (23) and (24). The tkplot() function can still be used to obtain the pictures of the relations among the variables given in Figure 3.

The next section will show how to perform the decomposition with causal_decompose() if the structure of $\mathbf{\Sigma }$ is not known and the goal is to carry out estimation and feedback testing from observed data.

Finding significant feedbacks with SIRE: an application to Italian macroeconomic data

As pointed out in the previous section, empirical evidence aside, the results of a decomposition based on binary matrices ${\mathbf{\Gamma }}^b$ and ${\mathbf{\Sigma }}^b$ must be considered as preliminary since they show only the potential links acting in the system. The effectiveness of these links demands a confirmation based on a sound empirical-evidence argument. In fact, the lack of significance of one or more of the feedbacks thus detected can alter the nature of the connections among the endogenous variables found by the preliminary decomposition, which is based only on the topological properties of matrices $\mathbf{\Gamma }$ and $\mathbf{\Sigma }$. In order to show how effective feedbacks operating in a model can be detected and tested, we have applied the functionalities of SIRE to the Klein model (see Klein (1950), (1950), and Greene (2003), (2003)). This model, originally conceived for the US economy, has been recast for the Italian economy. The Italian macroeconomic variables, mirroring the US counterparts, are available at http://dati.istat.it/. The given model is composed of $n=60$ observations on a quarterly basis and six equations explaining the following endogenous variables: consumption expenses for Italian families [C], added value [CP], private wages from dependent employment [WP], gross investment [I], gross capital stock [K], gross domestic product [GDP]. The model is specified as follows $$\begin{array}{}\text{(51)}& \begin{array}{rl}\left[\begin{array}{c}{\text{C}}_t\\ {\text{I}}_t\\ {\text{WP}}_t\\ {\text{GDP}}_t\\ {\text{CP}}_t\\ {\text{K}}_t\end{array}\right]=& {\mathbf{a}}_0+\left[\begin{array}{cccccc}0& {\gamma }_{12}& 0& 0& {\gamma }_{15}& 0\\ 0& 0& 0& 0& 0& {\gamma }_{26}\\ 0& {\gamma }_{32}& 0& {\gamma }_{34}& 0& 0\\ {\gamma }_{41}& {\gamma }_{42}& 0& 0& 0& 0\\ 0& 0& {\gamma }_{53}& 0& 0& 0\\ 0& {\gamma }_{62}& 0& 0& 0& 0\end{array}\right]\left[\begin{array}{c}{\text{C}}_t\\ {\text{I}}_t\\ {\text{WP}}_t\\ {\text{GDP}}_t\\ {\text{CP}}_t\\ {\text{K}}_t\end{array}\right]+\\ & +\left[\begin{array}{cccc}{a}_{11}& 0& 0& 0\\ a_{21}& 0& 0& 0\\ 0& 0& a_{34}& 0\\ 0& 0& a_{44}& 0\\ 0& 0& 0& a_{55}\\ 0& a_{62}& 0& 0\end{array}\right]\left[\begin{array}{c}{\text{CP}}_{t-1}\\ {\text{K}}_{t-1}\\ {\text{GDP}}_{t-1}\\ {\text{T}}_t\end{array}\right]+\left[\begin{array}{c}{e}_{\text{C}}\\ e_{\text{I}}\\ e_{\text{WP}}\\ e_{\text{GDP}}\\ e_{\text{CP}}\\ e_{\text{K}}\end{array}\right]\end{array}\end{array}$$ where ${\mathbf{a}}_0$ is the intercept vector. As equation (51) shows, the set of predetermined variables includes one exogenous variable, taxes [${\text{T}}_t$], and three lagged endogenous variables, that is: the one-lagged added value [${\text{CP}}_{t-1}$], the one-lagged gross capital stock [${\text{K}}_{t-1}$], and the one-lagged gross domestic product [${\text{GDP}}_{t-1}$]. We first load the data into the R workspace.

> data(macroIT)

Following Greene (2003), the model equations have been estimated with 3SLS by using the R package systemfit (Henningsen and Hamann (2017), (2017)). The one-lagged capital stock [${\text{K}}_{t-1}$], [${\text{T}}_t$], [${\text{CP}}_{t-1}$], and [${\text{GDP}}_{t-1}$] have been employed as instrumental variables. Matrix $\mathbf{\Sigma }$, if the user does not specify its structure, is estimated by using the covariance matrix of the structural residuals. The function causal_decompose() can be also employed to estimate both the model via 3SLS and the $\mathbf{\Sigma }$ matrix, and yields three matrices: $\mathbf{C}$, ${\mathbf{\Psi }}_1$, and ${\mathbf{\Psi }}_0$. The first two include the coefficients associated to variables affected by feedback loops, operating either at a deterministic level or induced by error terms, the third contains the coefficients associated to variables playing a causal role in the system.

This version of causal_decompose() takes the following arguments:

• data: data frame containing all the variables in the equations.
• eq.system: list containing all the equations, as in systemfit.
• resid.est: denotes the method used to estimate $\mathbf{\Sigma }$, on the basis of 3SLS residuals; this method is specified in systemfit.
• instruments: set of instruments used to estimate the model, introduced either as a list or as a character vector, as in systemfit.
• sigma.in: not appropriate to empirical context, set to NULL.

The output of this function is a list containing the following objects:

• eq.system: the same list of equations provided as input.
• gamma, C, Psi0, Psi1, A, and Sigma: respectively matrices $\mathbf{C}$, ${\mathbf{\Psi }}_0$, ${\mathbf{\Psi }}_1$, $\mathbf{A}$, and $\mathbf{\Sigma }$.
• systemfit: the output of the systemfit() function used to estimate the model.
• all.graph: the DAG object for the undecomposed path diagram.
• dec.graph: the DAG object for the decomposed path diagram.
• path: the data-set containing all the paths/relations among the endogenous variables, along with their classification (i.e., causal, interdependent). The graph highlights which interdependent relations work at a systematic level and which are induced by the effect of correlations among residuals).

The code below performs the decomposition using the macroIT data

> #system of equations
> eq.system <- list(eq1 <- C ~  CP  + I +  CP_1  ,
+                   eq2 <- I ~ K + CP_1,
+                   eq3 <- WP ~  I + GDP +  GDP_1,
+                   eq4 <- GDP ~ C + I + GDP_1,
+                   eq5 <- CP ~   WP  + T,
+                   eq6 <- K ~ I + K_1)
> #instruments
> instruments <- ~ T +  CP_1 + GDP_1 + K_1
> #decomposition
> dec.macroIT <- causal_decompose(data = macroIT,
+                                 eq.system = eq.system,
+                                 resid.est = "noDfCor",
+                                 instruments = instruments)

Table 1 shows the results of the model estimation. Since some coefficients are not statistically significant (such as the coefficient associated to $[\text{I}]$ in the equation explaining $[\text{C}]$ and the coefficient associated to $[\text{GDP}]$ in the equation explaining $[\text{WP}]$), the model has been re-estimated and the coefficient matrix associated to the explicative endogenous variables decomposed again.

> #system of equations
> eq.system <- list(eq1 <- C ~ CP + CP_1  ,
+                   eq2 <- I ~ K,
+                   eq3 <- WP ~  I +  GDP_1,
+                   eq4 <- GDP ~ C + I + GDP_1,
+                   eq5 <- CP ~   WP  + T,
+                   eq6 <- K  ~ I + K_1)
> #instruments
> instruments <- ~ T +  CP_1 + GDP_1 + K_1
> #decomposition
> dec.macroIT.new <- causal_decompose(data = macroIT,
+                                     eq.system = eq.system,
+                                     resid.est = "noDfCor",
+                                     instruments = instruments)

The results of the last estimation process are shown in Table 2. Looking at the Theil inequality indexes (Theil (1961), (1961)) reported in the last column of the table, we can see that the estimated equations fit the data very well. In fact, all Theil indexes are close to zero.

Table 1: Macroeconomic model: preliminary estimates with 3SLS. The last column shows the Theil index for each model equation.
$\bullet$: significant at level $\alpha$=0.1
*: significant at level $\alpha$=0.05
**: significant at level $\alpha$=0.01
***: significant at level $\alpha$=0.001

	${\mathbf{a}}_0$	${\text{C}}_t$	${\text{I}}_t$	${\text{WP}}_t$	${\text{GDP}}_t$	${\text{CP}}_t$	${\text{K}}_t$	${\text{T}}_t$	${\text{GDP}}_{t-1}$	${\text{CP}}_{t-1}$	${\text{K}}_{t-1}$	Theil
${\text{C}}_t$	$\underset{(11.97)}{5.50}$	-	$\underset{(0.19)}{-0.1}$	-	-	${\underset{(0.09)}{10.1}}^{\ast \ast \ast }$	-	-	-	${\underset{(0.099)}{-0.35}}^{\ast \ast \ast }$	-	0.0073
${\text{I}}_t$	$\underset{(7.627)}{10.28}$	-	-	-	-	-	${\underset{(0.089)}{0.72}}^{\ast \ast \ast }$	-	-	$\underset{(0.04)}{0.005}$	-	0.0115
${\text{WP}}_t$	${\underset{(28.87)}{-320.17}}^{\ast \ast \ast }$	-	${\underset{(0.453)}{-1.93}}^{\ast \ast \ast }$	-	$\underset{(0.277)}{0.46}$	-	-	-	${\underset{(0.251)}{1.04}}^{\ast \ast \ast }$	-	-	0.0288
${\text{GDP}}_t$	$\underset{(10.732)}{4.55}$	${\underset{(0.098)}{1.02}}^{\ast \ast \ast }$	${\underset{(0.181)}{0.37}}^{\ast }$	-	-	-	-	-	${\underset{(0.075)}{0.31}}^{\ast \ast \ast }$	-	-	0.004
${\text{CP}}_t$	${\underset{(10.848)}{144.49}}^{\ast \ast \ast }$	-	-	${\underset{(0.045)}{0.49}}^{\ast \ast \ast }$	-	-	-	${\underset{(0.373)}{3.62}}^{\ast \ast \ast }$	-	-	-	0.0074
${\text{K}}_t$	${\underset{(1.599)}{-5.69}}^{\ast \ast \ast }$	-	${\underset{(0.073)}{0.68}}^{\ast \ast \ast }$	-	-	-	-	-	-	-	${\underset{(0.053)}{0.49}}^{\ast \ast \ast }$	0.0063

Table 2: Macroeconomic model: final estimates with 3SLS. The last column shows the Theil index for each model equation.
$\bullet$: significant at level $\alpha$=0.1
*: significant at level $\alpha$=0.05
**: significant at level $\alpha$=0.01
***: significant at level $\alpha$=0.001

	${\mathbf{a}}_0$	${\text{C}}_t$	${\text{I}}_t$	${\text{WP}}_t$	${\text{GDP}}_t$	${\text{CP}}_t$	${\text{K}}_t$	${\text{T}}_t$	${\text{GDP}}_{t-1}$	${\text{CP}}_{t-1}$	${\text{K}}_{t-1}$	Theil
${\text{C}}_t$	$\underset{(9.17)}{10.06}$	-	-	-	-	${\underset{(0.115)}{1.02}}^{\ast \ast \ast }$	-	-	-	${\underset{(0.111)}{-0.39}}^{\ast \ast \ast }$	-	0.0076
${\text{I}}_t$	${\underset{(2.202)}{11.22}}^{\ast \ast \ast }$	-	-	-	-	-	-	-	-	${\underset{(0.027)}{0.73}}^{\ast \ast }$	-	0.0114
${\text{WP}}_t$	${\underset{(26.387)}{-299.12}}^{\ast \ast \ast }$	-	${\underset{(0.444)}{-1.65}}^{\ast \ast \ast }$	-	-	-	-	-	${\underset{(0.134)}{1.39}}^{\ast \ast \ast }$	-	-	0.0304
${\text{GDP}}_t$	$\underset{(8.89)}{3.04}$	${\underset{(0.118)}{1.09}}^{\ast \ast \ast }$	${\underset{(0.112)}{0.39}}^{\ast \ast }$	-	-	-	-	-	${\underset{(0.076)}{0.28}}^{\ast \ast \ast }$	-	-	0.0042
${\text{CP}}_t$	${\underset{(9.554)}{142.70}}^{\ast \ast \ast }$	-	-	${\underset{(0.034)}{0.48}}^{\ast \ast \ast }$	-	-	-	${\underset{(0.299)}{3.68}}^{\ast \ast \ast }$	-	-	-	0.0073
${\text{K}}_t$	${\underset{(1.606)}{-5.53}}^{\ast \ast \ast }$	-	${\underset{(0.074)}{0.67}}^{\ast \ast \ast }$	-	-	-	-	-	-	-	${\underset{(0.053)}{0.49}}^{\ast \ast \ast }$	0.0062

The estimated covariance matrix of the structural error terms is given by $$\hat{\mathbf{\Sigma }}=\left[\begin{array}{cccccc}10.93& & & & & \\ -2.51& 2.61& & & & \\ 10.75& -5.04& 52.31& & & \\ -7.55& 1.55& 3.66& 7.15& & \\ -9.6& 4.27& -19.73& 6.07& 15.08& \\ 0.43& -0.68& 0.53& -0.09& -0.68& 0.81\end{array}\right]$$ while matrices $\mathbf{C}+{\mathbf{\Psi }}_1$ and ${\mathbf{\Psi }}_0$ turn out to be $$\begin{array}{}\text{(52)}& \begin{array}{rl}\mathbf{C}+{\mathbf{\Psi }}_1& =\left[\begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0.73\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0.67& 0& 0& 0& 0\end{array}\right]+\left[\begin{array}{cccccc}0& 0& 0& 0& 1.02& 0\\ 0& 0& 0& 0& 0& 0\\ 0& -1.65& 0& 0& 0& 0\\ 1.09& 0.39& 0& 0& 0& 0\\ 0& 0& 0.48& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right]\\ & =\left[\begin{array}{cccccc}0& 0& 0& 0& 1.02& 0\\ 0& 0& 0& 0& 0& 0.73\\ 0& -1.65& 0& 0& 0& 0\\ 1.09& 0.39& 0& 0& 0& 0\\ 0& 0& 0.48& 0& 0& 0\\ 0& 0.67& 0& 0& 0& 0\end{array}\right]\end{array}\end{array}$$ $$\begin{array}{}\text{(53)}& {\mathbf{\Psi }}_0=\mathbf{0}\end{array}$$ The matrix in Equation (52) embodies all the coefficients associated to variables involved in feedback loops, while matrix (53) includes those associated to variables playing a causal role. Looking at (52) we find a direct feedback between variables [I] and [K], while the variables of the pairs [I, WP], [I, GDP], [C, GDP], [CP, C], and [CP, WP] are directly linked (a black arrow connects the variables of each pair) as well as explained by equations with correlated errors. Accordingly, the variables of each pair may be internally connected by feedback loops.

The goal of our testing procedure will be to bring out which of these feedbacks, being significant, are truly effective.

Figure 5 depicts the links operating in this model, using the function tkplot() of the igraph package. In this figure, a unidirectional arrow denotes that a variable is explicative for another. If two variables are explicative one for the other, a direct feedback loop exists, depicted as two red arrows going in opposite directions. Instead, a red, dashed, curved, two-headed arrow between two variables indicates the existence of a feedback induced by error correlation.

> tkplot(dec.macroIT.new$dec.graph)

Figure 5: Path diagram of the macroeconomic model. Unidirectional arrows denote that one variable is explicative for another. The two red unidirectional arrows denote the presence of a direct feedback. The red, dashed, curved, double-headed arrows between pairs of variables denote feedback loops induced by error correlation.

Testing for feedback effects

The significance of these loops has been investigated by using the function feedback_ml() which performs the Wald test given in (50). The 3SLS parameter estimates have been used as preliminary estimates to obtain the maximum likelihood (ML) estimates of the parameters needed to build the test statistic. In particular, in order to reach the global maximum of the log-likelihood, the initial 3SLS parameter estimates have been randomly perturbed a certain number of times. The optimizer chosen for the scope is included in the Rsolnp package where the function gosolnp is specially designed for the randomization of starting values. The function feedback_ml() takes the following arguments:

• data: data frame containing all the variables in the equations.
• out.decompose: the output from the previous causal decomposition which is called by using the command causal_decompose().
• lb and ub: upper and lower bound of the parameter space (as in gosolnp).
• nrestarts, nsim and seed.in: parameters tuning the number of random initializations (as in gosolnp).

The output of this function is a list containing the following objects:

• rho.est: a data frame containing the estimated feedback loops for a given equation. The first column of this data frame, feedback eqn., provides the indexes of the equations involved in the feedback loop with the equation given in input, while the coefficients associated to the explicative endogenous for the equation in question are shown in the column rho.est.
• loglik: the estimated log-likelihood of the best model.
• theta.hessian: the estimated Hessian matrix $\hat{\mathbf{I}}$.
• rho.jacobian: the estimated Jacobian matrix $\hat{\mathbf{J}}$.
• wald: the value of the Wald test statistic $W$.

As an example, let us assume that the interest is in testing the significance of the feedbacks affecting the second equation, explaining the endogenous variable $[\text{I}]$. According to the previous analysis, this variable is connected to $[\text{K}]$ by a bidirectional link.

The Wald test for the significance of this feedback is performed by using the function feedback_ml() specified as follows

> test.E2=feedback_ml(data = macroIT,
+                     out.decompose = dec.macroIT.new,
+                     eq.id = 2,
+                     lb = min(dec.macroIT.new$Sigma) - 10,
+                     ub = max(dec.macroIT.new$Sigma) + 10,
+                     nrestarts = 10,
+                     nsim = 20000,
+                     seed.in = 1)

By visualizing the estimate of $\mathit{\rho }$ and the Wald statistic

> test.E2$rho.tbl
  Feedback eqn.   rho.est
1             6 0.1641469

> test.E2$wald
         [,1]
[1,] 4.115221

we can see that the existence of a feedback loop between [I] and [K] is confirmed.

Table 3 shows the results of the test for all the equations of the model. Looking at the $p$-values we conclude that all feedbacks are significant except the ones involving [CP] and [GDP]. For what concerns [CP], it is explained by [WP] without a feedback effect from the latter to the former. Regarding [GDP], which is affected by feedback effects, a deeper analysis is required in order to understand which of its two explicative variables [C] and [I] (if not both) are responsible for it. To this end, we have applied the Wald statistic given in (50) which leads us to conclude that only [C] is involved in a feedback loop with [GDP].

Table 3: Macroeconomic model: tests for feedback effects for the final model. Joint $W$ denotes the Wald statistic used to test the set of feedback loops affecting a given variable (see (45). Singular $W$ denotes the Wald statistic used to test the feedback effect between two specific variables (see (50)).

Equation	Feedback Variable	Joint $W$	$p$-value	Singular $W$	$p$-value
$\text{C}$	$\text{CP}$	386.6	$<.001$
$\text{I}$	$\text{K}$	4.115	0.042	-	-
$\text{WP}$	$\text{I}$	25.55	$<.001$	-	-
$\text{GDP}$	$\text{C}$	95.368	$<0.0001$	84.315	$<0.0001$
	$\text{I}$			0.352	0.553
$\text{CP}$	$\text{WP}$	0.046	0.831	-	-
$\text{K}$	$\text{I}$	19.595	<0.0001	-	-

In the end, the path diagram fully describing the recurrent and interdependent relationships in the model is displayed in Figure 6.

Figure 6: Path diagram of the modified macroeconomic model after testing for feedback effects. Black arrows denote causal link $({\mathbf{\Psi }}_0)$, red arrows denote interdependent links $(\mathbf{C})$, black arrows and red dashed arrows denote interdependent links induced by the correlation of the error terms $({\mathbf{\Psi }}_1)$.

5 Discussion

The set of functions worked out in the paper allows a system of simultaneous equations to be split into recursive and/or interdependent subsystems. The user can rely on causal_decompose() in two ways: to assess the presence of interdependent relations with a known structure of correlation among the error terms, or to estimate the whole model in presence of empirical data.

The significance of the feedback loops operating in the model is tested with a Wald test using the feedback_ml() function. The 3SLS parameter estimates are used as preliminary estimates to obtain the maximum likelihood ones, which are needed to build the test.

As for the rationale of our procedure, which rests on a properly devised test, it is worth taking into account the considerable concern raised recently in the statistical community about the use of significance testing (see Wasserstein and Lazar (2016), (2016)). In this connection, in order to avoid improper use of $p$-values and significance-related results, it may be worth addressing the issue of detecting feedback mechanisms in a simultaneous equations model with different approaches. Among them, the construction of confidence intervals and the employment of Bayesian methods look particularly promising for future investigation.

Moving now on more technical notes:

• The ML estimation is performed by concentrating the likelihood with respect to the 3SLS estimates of $\mathbf{A}$ in Equation (1), to reduce the computation required to otherwise re-estimate parameters that are unnecessary for the computation of the feedback effect.
• As far as the error covariance matrix $\mathbf{\Sigma }$ is concerned, in the current formulation of the test its estimate $\hat{\mathbf{\Sigma }}$ is not involved by itself in any testing sub-routine (in fact, all of its elements are retained after the 3SLS step), and computing the related matrix of standard errors is therefore of secondary importance. However, if a matrix normal distribution is hypothesized on $\mathbf{E}$, then the distribution of $\hat{\mathbf{\Sigma }}$ turns out to be a $L$-dimensional Wishart with $T$ degrees of freedom and scale matrix $\mathbf{\Sigma }$. Thus, the variance of its elements can be calculated straightforwardly (see Gupta and Nagar (1999), (1999)).

6 Proof of relevant formulas

In this Appendix we provide the proofs of some relevant formulas of the paper.

1. Let $\mathbf{\Sigma }$ and $\mathbf{R}$ be defined as in Section 2. Then, the proof that ${\mathbf{\Sigma }}^b(\mathbf{I}+\mathbf{R})$ is the binary matrix associated to $\mathbf{\Sigma }{(\mathbf{I}-\mathbf{\Gamma })}^{-1}$ is based on the following two theorems.

   If two conformable matrices, $\mathbf{A}$ and $\mathbf{B}$, are such that $$\mathbf{A}=\mathbf{A}\ast \mathbf{H}\mathbf{B}=\mathbf{B}\ast \mathbf{K}$$ then the binary matrix associated to $\mathbf{A}\mathbf{B}$ is $(\mathbf{H}\mathbf{K}{)}^b$.$◼$

   If a non-singular matrix $\mathbf{A}$ is such that $$\mathbf{A}=\mathbf{A}\ast \mathbf{H}$$ where $\mathbf{H}$ is a given binary matrix, then $$({\mathbf{A}}^{-1})\ast (\mathbf{I}+\sum _{n=1}^{N-1}{\mathbf{H}}^n{)}^b=({\mathbf{A}}^{-1})$$ where $N$ is the matrix dimension.$◼$

   Now, upon noting that $$\mathbf{(}\mathbf{I}\mathbf{-}\mathbf{\Gamma }\mathbf{)}=\mathbf{(}\mathbf{I}\mathbf{-}\mathbf{\Gamma }\mathbf{)}\ast (\mathbf{I}-{\mathbf{\Gamma }}^b),$$ reference to Theorem [th:t2] leads to conclude that $$\begin{array}{}\text{(54)}& \mathbf{(}\mathbf{I}\mathbf{-}\mathbf{\Gamma }{\mathbf{)}}^{\mathbf{-}\mathbf{1}}=\mathbf{(}\mathbf{I}\mathbf{-}\mathbf{\Gamma }{\mathbf{)}}^{\mathbf{-}\mathbf{1}}\ast \mathbf{(}\mathbf{I}\mathbf{+}\mathbf{R}\mathbf{)}\end{array}$$ Next, taking into account that ${\mathbf{\Sigma }}^b$ and $(\mathbf{I}+\mathbf{R})$ are the binary counterparts of the $\mathbf{\Sigma }$ matrix and $\mathbf{(}\mathbf{I}\mathbf{-}\mathbf{\Gamma }{\mathbf{)}}^{\mathbf{-}\mathbf{1}}$ reference to Theorem entails the following \[\boldsymbol{\Sigma}(\mathbf I-\boldsymbol\Gamma)=[\boldsymbol{\Sigma}(\mathbf I-\boldsymbol\Gamma)]\ast [\boldsymbol{\Sigma}^b(\mathbf I+\mathbf R)].\]

2. The proof that $\mathbf{C}$ and $\mathbf{F}$, defined as in Section 3, satisfy the following relationship $$\begin{array}{}\text{(55)}& \mathbf{C}\ast \mathbf{F}=\mathbf{C}\end{array}$$ hinges on a preliminary result given in the following theorem.

   The matrices $\mathbf{C}$ and $\mathbf{I}+\mathbf{R}$ satisfy the following relationship $${\mathbf{C}}^b\ast (\mathbf{I}+\mathbf{R})={\mathbf{C}}^b$$

   Proof

   Taking into account that the Hadamard product is both commutative $(\mathbf{A}\ast \mathbf{B}=\mathbf{B}\ast \mathbf{A})$ and idempotent for binary matrices $({\mathbf{A}}^b\ast {\mathbf{A}}^b={\mathbf{A}}^b)$, and being $\mathbf{\Gamma }$ hollow, the following holds $${\mathbf{\Gamma }}^b\ast \mathbf{I}=\mathbf{0},$$ simple computations yield $$\begin{array}{}\text{(56)}& \begin{array}{r}{\mathbf{C}}^b\ast (\mathbf{I}+\mathbf{R})={\mathbf{\Gamma }}^b\ast \mathbf{R}\ast (\mathbf{I}+\mathbf{R})={\mathbf{\Gamma }}^b\ast \mathbf{R}\ast \mathbf{I}+{\mathbf{\Gamma }}^b\ast \mathbf{R}\ast \mathbf{R}={\mathbf{C}}^b\end{array}\end{array}$$ $◼$
   Now, consider the following theorem (where the symbol $\mathbf{A}\ge \mathbf{0}$ denotes that all the elements of matrix $\mathbf{A}$ are non negative numbers):

   Let $\mathbf{B}\ge \mathbf{0}$ and ${\mathbf{A}}^b\ast {\mathbf{B}}^b={\mathbf{A}}^b$. If $\mathbf{C}\ge \mathbf{0}$, then $${\mathbf{A}}^b\ast {\mathbf{(}\mathbf{B}\mathbf{+}\mathbf{C}\mathbf{)}}^b={\mathbf{A}}^b$$ $◼$

   Given this premise, we can now prove (55). To this end, let us write ${\mathbf{\Sigma }}^b$ as follows $$\begin{array}{}\text{(57)}& (\mathbf{I}+\mathbf{\Delta })={\mathbf{\Sigma }}^b\end{array}$$ where $\mathbf{\Delta }\ast \mathbf{I}=\mathbf{0}$ is a hollow matrix, and note that, in light of (57) and (54), the binary matrix associated to $\mathbf{F}$ is, according to Theorem [th:t1], given by $${\mathbf{F}}^b=[(\mathbf{I}+\mathbf{\Delta })(\mathbf{I}+\mathbf{R}){]}^b$$ Next, use of Theorems [th:t3] and [th:t4], yields the following $${\mathbf{C}}^b\ast {\mathbf{F}}^b={\mathbf{C}}^b\ast [(\mathbf{I}+\mathbf{\Delta })(\mathbf{I}+\mathbf{R}){]}^b={\mathbf{C}}^b\ast [(\mathbf{I}+\mathbf{R})+\mathbf{\Delta }(\mathbf{I}+\mathbf{R}){]}^b={\mathbf{C}}^b$$ as $(\mathbf{\Delta }(\mathbf{I}+\mathbf{R}){)}^b\ge \mathbf{0}$. This, in turn, entails that $$\begin{array}{}\text{(58)}& {\mathbf{C}}^b+{\mathbf{\Psi }}_1^b\ast {\mathbf{F}}^b=({\mathbf{C}}^b+{\mathbf{\Psi }}_1^b)\ast {\mathbf{F}}^b={\mathbf{\Gamma }}^b\ast {\mathbf{F}}^b\end{array}$$ which means that $\mathbf{C}+{\mathbf{\Psi }}_1\ast \mathbf{F}$ and $\mathbf{\Gamma }\ast \mathbf{F}$ have the same topological structure. $◼$

3. Proof of (37) . Formula (37) can be proved as follows. First, note that matrix ${\mathbf{\Gamma }}^{\ast }$ weighting the current endogenous explicative variables in the model (30), (31) can be expressed as $$\begin{array}{}\text{(59)}& {\mathbf{\Gamma }}^{\ast }={\mathbf{P}}_j\mathbf{\Gamma }{\mathbf{P}}_j\end{array}$$ where ${\mathit{P}}_j$ is a permutation matrix obtained from an identity matrix by interchanging its first row with its $j$-th row. Then note that $${\mathbf{\Gamma }}^{\ast }=\left[\begin{array}{cc}0& {\mathit{\gamma }}_j^{\prime }\\ \mathit{\eta }& \mathbf{(}\mathbf{I}-{\mathbf{\Gamma }}_{\mathit{o}})\end{array}\right]$$ and that $$(\mathbf{I}-{\mathbf{\Gamma }}^{\ast }{)}^{-1}=\left[\begin{array}{cc}1+{\mathit{\gamma }}_j^{\prime }{\mathbf{L}}^{-1}\mathit{\eta }& {\mathit{\gamma }}_j^{\prime }{\mathbf{L}}^{-1}\\ {\mathbf{L}}^{-1}\mathit{\eta }& {\mathbf{L}}^{-1}\end{array}\right]$$ , where $$\mathbf{L}=\mathbf{I}-{\mathbf{\Gamma }}_o-\mathit{\eta }{\mathit{\gamma }}_j^{\prime }=(\mathbf{I}-\mathbf{G})$$ Accordingly $$\frac{1}{{\sigma }_{jj}}{\mathbf{M}}_1^{\prime }(\mathbf{I}-{\mathbf{\Gamma }}^{\ast }{)}^{-1}\mathbf{\Sigma }{\mathbf{e}}_1=(\mathbf{I}-\mathbf{G}{)}^{-1}\mathbf{d}={\mathit{\phi }}_j,$$ where ${\mathbf{e}}_j$ is the first elementary vector, $\mathbf{\Sigma }$, $\mathbf{G}$ and $\mathbf{d}$ are defined as in (32)and (36) respectively, and ${\mathbf{M}}_1$ is the selection matrix obtained from the identity matrix by deleting its first column. Now, taking into account that the following holds $$(\mathbf{I}-{\mathbf{\Gamma }}^{\ast })={\mathbf{P}}_j(\mathbf{I}-\mathbf{\Gamma }){\mathbf{P}}_j$$ in light of (59), and that the following proves true $$(\mathbf{I}-{\mathbf{\Gamma }}^{\ast }{)}^{-1}={\mathbf{P}}_j(\mathbf{I}-\mathbf{\Gamma }{)}^{-1}{\mathbf{P}}_j,$$ as ${\mathbf{P}}_j$ is both symmetric and orthogonal, some computations yield $$\begin{array}{}\text{(60)}& \begin{array}{rl}{\mathit{\phi }}_j& =\frac{1}{{\sigma }_{jj}}{\mathbf{M}}_1^{\prime }(\mathbf{I}-{\mathbf{\Gamma }}^{{\ast }^{\prime }}{)}^{-1}\mathbf{\Sigma }{\mathbf{e}}_1=\frac{1}{{\sigma }_{jj}}{\mathbf{M}}_1^{\prime }{\mathbf{P}}_j(\mathbf{I}-\mathbf{\Gamma }{)}^{-1}{\mathbf{P}}_j\mathbf{\Sigma }{\mathbf{P}}_j{\mathbf{P}}_j{\mathbf{e}}_1=\\ & =\frac{1}{{\sigma }_{jj}}{\mathbf{M}}_j^{\prime }(\mathbf{I}-\mathbf{\Gamma }{)}^{-1}\mathbf{\Sigma }{\mathbf{e}}_j=\frac{1}{{\sigma }_{jj}}{\mathbf{M}}_j^{\prime }{\mathbf{F}}^{\prime }{\mathbf{e}}_j\end{array}\end{array}$$ $◼$

4. Derivation of the log-likelihood for the model (40)-(42)
   The logarithm of the density in (42) is given by $$\begin{array}{}\text{(61)}& \ln f_{({ϵ}_j,{\mathit{ϵ}}_r)}=c-\frac{1}{2}\ln {\sigma }_{jj}-\frac{1}{2}\ln |\mathbf{\Omega }|-\frac{{ϵ}_j^2}{2{\sigma }_{jj}}-\frac{1}{2}{\mathit{ϵ}}_r^{\prime }{\mathbf{\Omega }}^{-1}{\mathit{ϵ}}_r\end{array}$$ where $c$ is a constant term. Now, upon noting that $$|\mathbf{J}|=\left|\frac{\partial ({ϵ}_j,{\mathit{ϵ}}_r{)}^{\prime }}{\partial (y_j,{\mathbf{y}}_r^{\prime })}\right|=1,$$ and assuming to operate with $N$ observations on the variables of interest, the log-likelihood function can be written as $$\begin{array}{}\text{(62)}& l=\sum _{t=1}^N{l}_{(y_j,{\mathbf{y}}_r^{\prime })}=k-\frac{N}{2}\ln {\sigma }_{jj}-\frac{N}{2}\ln |\mathbf{\Omega }|-\frac{{\mathit{\alpha }}^{\prime }\mathbf{H}\mathit{\alpha }}{2{\sigma }_{jj}}-\frac{1}{2}\text{tr}({\mathbf{\Xi }}^{\prime }{\mathbf{\Omega }}^{-1}\mathbf{\Xi }\mathbf{H})\end{array}$$ where $$\begin{array}{}\text{(63)}& \begin{array}{rl}{\mathit{\alpha }}^{\prime }& =[1,-{\stackrel{\mathbf{~}}{\mathit{\gamma }}}_j^{\prime },-{\stackrel{\mathbf{~}}{\mathbf{a}}}_j^{\prime }{\mathbf{S}}_j]\end{array}\end{array}$$

   $$\begin{array}{}\text{(64)}& \begin{array}{rl}\mathbf{\Xi }& =[-{\stackrel{\mathbf{~}}{\mathit{\phi }}}_j,\mathbf{I}+{\stackrel{\mathbf{~}}{\mathit{\phi }}}_j{\stackrel{\mathbf{~}}{\mathit{\gamma }}}_j^{\prime },{\stackrel{\mathbf{~}}{\mathit{\phi }}}_j{\stackrel{\mathbf{~}}{\mathbf{a}}}_j^{\prime }{\mathbf{S}}_j-\mathit{K}]\end{array}\end{array}$$

   $$\begin{array}{}\text{(65)}& \begin{array}{rl}{\mathit{\nu }}^{\prime }& =[y_j,{\mathbf{y}}_o,\mathbf{z}]\\ \mathbf{H}& =\left(\sum _{t=1}^N{\mathit{\nu }}_t{\mathit{\nu }}_t^{\prime }\right),\end{array}\end{array}$$ and $k$ is a constant term. Formula (62) can be obtained by noting that, in light of (40), the following holds $$\begin{array}{}\text{(66)}& {ϵ}_j=y_j-{\stackrel{\mathbf{~}}{\mathit{\gamma }}}_j^{\prime }{\mathbf{y}}_r-{\stackrel{\mathbf{~}}{\mathbf{a}}}_j^{\prime }{\mathbf{z}}_r=[1,-{\stackrel{\mathbf{~}}{\mathit{\gamma }}}_j^{\prime },-{\stackrel{\mathbf{~}}{\mathbf{a}}}_j^{\prime }{\mathbf{S}}_r]\left[\begin{array}{c}{y}_j\\ {\mathbf{y}}_r\\ \mathbf{z}\end{array}\right]={\mathit{\alpha }}^{\prime }\mathit{\nu }\end{array}$$ and that, according to (41), we have $$\begin{array}{}\text{(67)}& \begin{array}{rl}{\mathit{ϵ}}_r& ={\mathbf{y}}_r-\mathbf{K}\mathbf{z}-{\stackrel{\mathbf{~}}{\mathit{\phi }}}_j{ϵ}_j=\\ & ={\mathbf{y}}_r-\mathbf{K}\mathbf{z}-{\stackrel{\mathbf{~}}{\mathit{\phi }}}_j(y_j-{\stackrel{\mathbf{~}}{\mathit{\gamma }}}_j^{\prime }{\mathbf{y}}_r-{\stackrel{\mathbf{~}}{\mathbf{a}}}_j^{\prime }{\mathbf{S}}_r\mathbf{z})=\\ & =[-{\stackrel{\mathbf{~}}{\mathit{\phi }}}_j,\mathbf{I}+{\stackrel{\mathbf{~}}{\mathit{\phi }}}_j{\stackrel{\mathbf{~}}{\mathit{\gamma }}}_j^{\prime },{\stackrel{\mathbf{~}}{\mathit{\phi }}}_j{\stackrel{\mathbf{~}}{\mathbf{a}}}_j^{\prime }{\mathbf{S}}_r-\mathbf{K}]\left[\begin{array}{c}{y}_j\\ {\mathbf{y}}_r\\ \mathbf{z}\end{array}\right]=\mathbf{\Xi }\mathit{\nu }\end{array}\end{array}$$ This implies that $$\begin{array}{}\text{(68)}& {\mathit{ϵ}}_r^{\prime }{\mathbf{\Omega }}^{-1}{\mathit{ϵ}}_r=\text{tr}({\mathit{\nu }}^{\prime }{\mathbf{\Xi }}^{\prime }{\mathbf{\Omega }}^{-1}\mathbf{\Xi }\mathit{\nu })=\text{tr}({\mathbf{\Xi }}^{\prime }{\mathbf{\Omega }}^{-1}\mathbf{\Xi }\mathit{\nu }{\mathit{\nu }}^{\prime })\end{array}$$ $◼$