4.2 Non orthogonality condition

As shown earlier, under the orthogonality condition the OLS estimators and the classical tests save their distributions asymptotically. In real world data, however, there exists sometimes some dependence between the regressors and the error term, even in large samples, due to several reasons. If so, the regressors are called endogenous, in contrary when they are exogenous under the orthogonality condition.

We will discuss the common cases where the linear model inevitably contains such dependence.

4.2.1 Omission of relevant variables

Suppose that the DGP model contains the following regressors in which all the assumptions A1.1 to A9 are satisfied:

\[\begin{equation*} y_t=\beta_1+\beta_2x_{2t}+\beta_3x_{3t}+\beta_4x_{4t}+\varepsilon_t \end{equation*}\]

And its best estimated model is the following one:

\[\begin{equation} y_t=\widehat\beta_1+\widehat\beta_2x_{2t}+\widehat\beta_3x_{3t}+\widehat\beta_4x_{4t}+e_t \tag{4.12} \end{equation}\]

However, instead of estimating that model, the investigator estimates the following one with the variable \(x_4\) excluded:

\[\begin{equation} y_t=\tilde\beta_1+\tilde\beta_2x_{2t}+\tilde\beta_3x_{3t}+\tilde e_t \tag{4.13} \end{equation}\]

Since, by definition, the residual term should contain all the remaining factors that possibly affect the response variable and not included as regressors, so the residual term \(\tilde e\) in (4.13) should contain the effect of both \(e\), in (4.12), and the excluded regressor \(x_4\). Formally written:

\[\begin{equation} \tilde e_t=\gamma x_{4t}+e_t \tag{4.14} \end{equation}\]

Indeed, \(x_4\) and \(\varepsilon\) are independent according to the assumptions stated above. But what about the correlation that may exist between \(x_4\) and the other regressors?

In the previous chapter in which we discussed this case in more detail, it was proved that in a model with a missed relevant variable, the coefficients of the other ones would be biased unless the missed one is orthogonal to the others, which is rare when it occurs.
As \(x_4\) is correlated with \(\tilde e\) from (4.14), and, in general, correlated with the other regressors, then \(\tilde e\) will be in turn correlated with those regressors in (4.13), and hence the orthogonality condition will be violated.

First of all, and before any estimation occurs, the investigator should collect all the relevant variables that explain the best its variable of interest. If one or more variables have been omitted, then some relevant explanation has been lost. In other words, all that any estimation method can do is to somehow reformulate the theoretical explanation into the most closer model to the true one so that if some relevant variables are missed, the estimated model loses a relevant distance from the true model whatever the method used.

Solving the orthogonality problem essentially depends on whether the omitted variable (\(x_4\) in our example) is completely unknown or known without observation. These two cases will be treated using the instrumental variables method (IV) that will be discussed shortly. Fortunately, There exists an easy alternative for the latter case if it is possible to find a variable, with observations, that is highly correlated with the omitted variable in question and supposed to have the same effect on the response variable. This variable is called proxy variable, and we can just put it into the original model in place of the omitted variable and fit the model with the OLS as usual. It is worthwhile to mention that all the methods that can be used here give only consistent estimators and never unbiased.

4.2.2 Measurement errors in the regressors

Measurement errors are common when collecting economic data. Some examples are answers to survey data for the respondent’s lack of certainty or lack of understanding of the questions. Some macroeconomic variables can not be determined precisely as they are estimated via other variables, etc. Let us consider our true model \(y=X\beta+\varepsilon\), with all the classical assumptions satisfied. Suppose that the regressors in \(X\) and the dependent variable \(y\) are measured inaccurately (with errors). Since these variables in \(X\) and \(y\) are slightly different from those of the true model, we will denote the true variables by \(X^*\) and \(y^*\), and can be related by the following:

\[\begin{equation} \begin{cases} X=X^*+\nu\\ y=y^*+\upsilon \tag{4.15} \end{cases} \end{equation}\]

Where \(\nu\) and \(\upsilon\) are random variables IID with mean zero and fixed variances \(\sigma^2_{\nu}\) and \(\sigma^2_{\upsilon}\) respectively, and assumed to be independent of \(X^*\) and \(y^*\). So they can be considered as white noises.

Substituting (4.15) into the true model, we get:

\[\begin{align*} &y^*+\upsilon=\big(X^*+\nu\big)\beta+\varepsilon\implies\\ &y^*=X^*\beta+\nu\beta+\varepsilon+\upsilon \end{align*}\]

Since each of the random variables, \(\nu\beta\), \(\varepsilon\), and \(\upsilon\) assumed to have zero mean and fixed variance, they can be thus put together into one single variable that will denote \(\eta\):

\[\begin{equation} \eta=\nu\beta+\varepsilon-\upsilon \tag{4.16} \end{equation}\]

Then the correct specification of the true model should be:

\[\begin{equation} y^*=X^*\beta+\eta \end{equation}\]

Let us now check the orthogonality assumption using the model \(y=X\beta+\varepsilon\) with the corrupted variables \(X\) and \(y\):

\[\begin{align*} E(X^t\varepsilon)&=E\bigg(X^t(\eta-\nu\beta+\upsilon)\bigg)\\ &=E(X^t\eta)-\overbrace{E(X^t\nu)}^{\not=0}\beta+E(X^t\upsilon)\\ &\not=0 \end{align*}\]

This equation can not be equal to zero since, from (4.15), the quantity \(E(X^t\nu)\not=0\), hence the orthogonality condition is no longer satisfied, and the OLS estimators are no longer consistent.

4.2.2.1 Instrumental variables method (IV)

For simplification, we will illustrate this method with the following simple linear model with demeaned variables (without intercept):

\[\begin{equation*} y_t=\beta x_t+\varepsilon_t \end{equation*}\]

Suppose now that the orthogonality assumption A5 is violated \(E(x_t\varepsilon_t)\not=0\). The intuitive interpretation of this restriction is that the variation of the explanatory variable \(x\), which is supposed to explain appropriately the variation of the dependent variable \(y_t\), does not come entirely from \(x\) itself. This is because a small part of that variation is caused by the error term.

In other words, the parameter \(\beta\) in the above model is supposed to represent only the single effect of \(x\) on \(y\) while here represents also another effect from \(\varepsilon_t\) that comes indirectly from the joint variation with the regressor. Therefore, the IV method principle is to isolate the variation that is not correlated with the error term and, hence get consistent estimators.

This method includes another explanatory variable to isolate the effect side of the original regressor from the error term (in other words, clean the regressor). Intuitively, the suggested variable is the one that has the highest correlation with the regressor and simultaneously is independent of the error term, which satisfies the orthogonality condition.

Let \(z\) is the instrumental variable that has the properties noted above. Using the probability limit, these properties can be interpreted as follows:

\[\begin{equation*} \begin{cases} plim\bigg(\frac{1}{n}\sum x_tz_t\bigg)\not=0\\ plim\bigg(\frac{1}{n}\sum z_t\varepsilon_t\bigg)=0 \end{cases} \end{equation*}\]

Multiplying both sides of our initial model by \(\frac{1}{n}z_t\) then summingn we get:

\[\begin{equation*} \sum \frac{1}{n}y_tz_t=\beta\sum \frac{1}{n}x_tz_t+\overbrace{\sum \frac{1}{n}z_t\varepsilon_t}^{\approx0} \end{equation*}\]

Using the properties noted above, the estimator of \(\beta\) will be:

\[\begin{equation} \beta_{IV}=\frac{\sum \frac{1}{n}y_tz_t}{\sum \frac{1}{n}x_tz_t}=\frac{\sum y_tz_t}{\sum x_tz_t} \tag{4.17} \end{equation}\]

Or with the original variables:

\[\begin{equation*} \beta_{IV}=\frac{\sum (y_t-\overline y)(z_t-\overline z)}{\sum (x_t-\overline x)(z_t-\overline z)} \end{equation*}\]

Where \(\beta_{IV}\) is called the instrumental variable estimator.

That formula shows that this instrument is used to compute the parameter estimates and does not replace the original regressor unless the correlation between them tends to be perfect, hence return to the non-orthogonality condition. In more detail, the instrument should have a reasonable correlation with the original regressor to save the proper variation of this regressor but not a very high to still independent of the error term. For instance, if the instrument has a perfect linear combination with the regressor such as \(z_t=\alpha x_t\), then the orthogonality condition will be violated and will have the same inconsistent estimators as follows:

\[\begin{align*} \beta_{IV}&=\frac{\sum (y_t-\overline y)(z_t-\overline z)}{\sum (x_t-\overline x)(z_t-\overline z)}\\ &=\frac{\sum (y_t-\overline y)(\alpha x_t-\alpha\overline x)}{\sum (x_t-\overline x)(\alpha x_t-\alpha\overline x)}\\ &=\frac{\alpha\sum (y_t-\overline y)(x_t-\overline x)}{\alpha\sum (x_t-\overline x)^2}\\ &=\widehat\beta \end{align*}\]

It can be said that the best correlation is that equivalent to the variation of the regressor that is independent of the error term.

Another interesting interpretation of the IV estimator is that it can be seen as the parameter estimate resulted from the regression of \(y_t\) on \(\widehat x_t\), where \(\widehat x_t\), in turn, is the fitted variable resulted from the regression of the original regressor \(x_t\) on the instrument \(z_t\).

In more detail, suppose that we have the following two regressions expressed in terms of the fitted variables, \(\widehat y_t=\widehat\gamma\widehat x_t\) and \(\widehat x_t=\widehat\alpha z_t\), where \(\widehat\gamma=\frac{\sum y_t\widehat x_t}{\sum \widehat x^2_t}\) and \(\widehat\alpha=\frac{\sum x_tz_t}{\sum z^2_t}\), The instrumental estimator will be:

\[\begin{align*} \beta_{IV}&=\frac{\sum y_tz_t}{\sum x_tz_t}\\ &=\frac{\widehat\alpha\sum y_tz_t}{\widehat\alpha\sum x_tz_t}\\ &=\frac{\sum y_t(\overbrace{\widehat\alpha z_t}^{=\widehat x_t})}{\sum x_t(\underbrace{\widehat\alpha z_t}_{=\widehat x_t})}\\ &=\frac{\sum y_t\widehat x_t}{\sum x_t\widehat x_t} \end{align*}\]

Now, if we denote by \(e\) the residuals resulted from the regression of \(x\) on \(z\), then by using the normal equation \(\sum z_te_t=0\), the denominator of the above expression will be:

\[\begin{align*} \sum x_t\widehat x_t&=\sum (\widehat x_t+e_t)\widehat x_t\\ &=\sum \widehat x_t^2+\sum e_t\widehat x_t\\ &=\sum \widehat x_t^2+\sum e_t(\widehat\alpha z_t)\\ &=\sum \widehat x_t^2+\widehat\alpha\underbrace{\sum e_tz_t}_{=0}\\ &=\sum \widehat x_t^2 \end{align*}\]

Hence:

\[\begin{equation} \beta_{IV}=\frac{\sum y_tx_t}{\sum x_tz_t}=\frac{\sum y_t\widehat x_t}{\sum \widehat x_t^2}=\widehat\gamma \tag{4.18} \end{equation}\]

This crucial result means that the best way to eliminate the non-orthogonality problem is to replace the original regressor with its fitted variable resulting from its regression on the instrument. However, this does not mean that the best-fitted values in \(\widehat x_t\) are those that are more closer to the original values in \(x_t\) as usual, so if it does then, it would be better to use the regressor itself without getting the headache of searching for its estimator. So here, the best estimation is the one based on the selected instrument that should be strongly correlated to the regressor and at the same time highly uncorrelated with the error term.

In more general, if we collect all the possible relevant variables that can significantly explain the variable \(x_t\) and each of which satisfies the orthogonality condition with the error of the original model \(\varepsilon\), then those variables can be used all together as instruments. That is why the number of instrumental variables can be greater than the number of regressors.

In conclusion, to satisfy the orthogonality condition, the instruments should be independent of the error term, so if they are, then the same holds for any of their combinations. The best combination thus is the one that provides the best explanation of the original regressor \(x_t\), which is that used by OLS estimation that gives the fitted values \(\widehat x_t\).

Now the previous analysis carries over the case of the general linear model. Suppose, for instance, that some regressors in the matrix \(X\) have been measured with errors, and let us consider again our original model:

\[\begin{equation*} y=X\beta+\varepsilon \end{equation*}\]

Let \(Z\) be an \((n\times k)\) matrix of the instrumental variables, which means that the number of instruments is equal to the number of regressors. Now, let us premultiply the original model by \(Z^t\) as follows:

\[\begin{equation*} Z^ty=Z^tX\beta+Z^t\varepsilon \end{equation*}\]

As the regressors are assumed to be stochastic, the instruments can also be assumed to be stochastic and, in addition, independent of the error term. Thus, the matrix \(Z\) should have the following properties:

• Orthogonality condition: \(plim\bigg(\frac{1}{n}Z^t\varepsilon\bigg)=0\).
• Stability condition with \(X\): \(plim\bigg(\frac{1}{n}Z^tX\bigg)=Q_{XZ}\).
• Stability condition with \(y\): \(plim\bigg(\frac{1}{n}Z^ty\bigg)=Q_{yZ}\).

The first property stands for the orthogonality condition. The second one means that the dot product of \(Z\) and \(X\) must converge in order to be invertible. The last one also ensures the convergence to a fixed matrix.

Now by including the plim operator in the above transformed model, we get:

\[\begin{equation*} plim\big(\frac{1}{n}Z^ty\big)=plim\big(\frac{1}{n}Z^tX\big)\beta+plim\big(\frac{1}{n}Z^t\varepsilon\big) \end{equation*}\]

Using the properties mentioned above, the last expression will be:

\[\begin{equation} Q_{Zy}=Q_{ZX}\beta_{IV}\implies\beta_{IV}=\big(Q_{ZX}\big)^{-1}Q_{Zy} \tag{4.19} \end{equation}\]

In practice, the matrices \(Q_{ZX}\) and \(Q_{Zy}\) are unknown, but they can be approximated respectively by \(\frac{1}{n}Z^tX\) and \(\frac{1}{n}Z^ty\). If so, the expression (4.19) will be rewritten:

\[\begin{align} \beta_{IV}&=\bigg(\frac{1}{n}Z^tX\bigg)^{-1}\bigg(\frac{1}{n}Z^ty\bigg)\notag\\ &=\big(Z^tX\big)^{-1}\big(Z^ty\big) \tag{4.20} \end{align}\]

Which is a consistent estimator, because:

\[\begin{align*} \beta_{IV}&=\big(Z^tX\big)^{-1}\big(Z^ty\big)\\ &=\big(Z^tX\big)^{-1}Z^t(X\beta+\varepsilon)\\ &=\beta+\big(Z^tX\big)^{-1}Z^t\varepsilon \end{align*}\]

Including the plim operator:

\[\begin{align*} plim(\beta_{IV})&=\beta+plim\bigg(\frac{1}{n}Z^tX\bigg)^{-1}\overbrace{plim\bigg(\frac{1}{n}Z^t\varepsilon\bigg)}^{=0}\\ &=\beta \end{align*}\]

As discussed before, this estimator vector can be obtained by regressing each regressor (whether is independent of the error term or not) on all the instruments and compute the fitted values in the first step, then regress the dependent variable on these fitted variables in the second step. This is the same method as IV but also called Two satge least squares 2SLS that will be discussed in more detail later on.

Distribution of \(\beta_{IV}\):

As the assumption of orthogonality is satisfied now and \(\beta_{IV}\) is determined asymptotically, then the central limit theorem can be used to get the asymptotic distribution. The expression (4.20) can be rewritten as follows:

\[\begin{align*} \beta_{IV}&=\big(Z^tX\big)^{-1}\big(Z^ty\big)\\ &=\big(Z^tX\big)^{-1}Z^t\big(X\beta+\varepsilon\big)\\ &=\beta+\big(Z^tX\big)^{-1}Z^t\varepsilon\\ \beta_{IV}-\beta&=\big(Z^tX\big)^{-1}Z^t\varepsilon\\ \sqrt n(\beta_{IV}-\beta)&=\bigg(\frac{1}{n}Z^tX\bigg)^{-1}\frac{1}{\sqrt n}Z^t\varepsilon \end{align*}\]

Using the plim and the previous results, the asymptotic distribution of \(\frac{1}{\sqrt n}Z^t\varepsilon\) will be:

\[\begin{equation*} \frac{1}{\sqrt n}Z^t\varepsilon\overset{d}{\rightarrow} \rm N(0, \sigma^2Q_{ZZ}) \end{equation*}\]

Hence the asymptotic distribution of \(\bigg(\frac{1}{n}Z^tX\bigg)^{-1}\frac{1}{\sqrt n}Z^t\varepsilon\) will be:

\[\begin{equation*} \bigg(\frac{1}{n}Z^tX\bigg)^{-1}\frac{1}{\sqrt n}Z^t\varepsilon\overset{d}{\rightarrow}\rm N\bigg(0,\sigma^2Q_{ZX}^{-1}Q_{ZZ}Q_{XZ}^{-1}\bigg) \end{equation*}\]

That is:

\[\begin{equation} \sqrt n(\beta_{IV}-\beta)\overset{d}{\rightarrow}\rm N\bigg(0,\sigma^2Q_{ZX}^{-1}Q_{ZZ}Q_{XZ}^{-1}\bigg) \tag{4.21} \end{equation}\]

Then we can derive the asymptotic distribution of \(\beta_{IV}\):

\[\begin{equation} \beta_{IV}\overset{d}{\rightarrow}N\bigg(\beta,\frac{\sigma^2}{n}Q_{ZX}^{-1}Q_{ZZ}Q_{XZ}^{-1}\bigg) \tag{4.22} \end{equation}\]

As before, the unknown matrices \(Q_{ZX}\), \(Q_{ZZ}\), and \(Q_{XZ}\) can be replaced by their approximations in the sample respectively \(\bigg(\frac{1}{n}Z^tX\bigg)\), \(\bigg(\frac{1}{n}Z^tZ\bigg)\), and \(\bigg(\frac{1}{n}X^tZ\bigg)\).The estimated variance matrix thus will be:

\[\begin{align} \widehat{Var}(\beta_{IV})&=\frac{1}{n}\widehat\sigma^2\bigg(\frac{1}{n}Z^tX\bigg)^{-1}\bigg(\frac{1}{n}Z^tZ\bigg)\bigg(\frac{1}{n}X^tZ\bigg)^{-1}\notag\\ &=\widehat\sigma^2\bigg(Z^tX\bigg)^{-1}\bigg(Z^tZ\bigg)\bigg(X^tZ\bigg)^{-1} \tag{4.23} \end{align}\]

Note that the estimated variance \(\widehat\sigma^2\) still computed by the usual formula \(\widehat\sigma^2=\frac{\sum e_{IV}^2}{n-k}\), or just \(\widehat\sigma^2=\frac{\sum e_{IV}^2}{n}\) since they are asymptotically equal. However, the residuals \(e_{IV}\) are computed by using \(\beta_{IV}\) (and not \(\widehat\beta\)), that is obtained from the model \(y=\beta_{IV}X+e_{IV}\).

Now it has been proved that the IV estimators are consistent but what about their efficiency compared to the OLS estimators (with the variance matrix \(Var(\widehat\beta)=\sigma^2(X^tX)^{-1}\))?. So let us compute the difference:

\[\begin{equation*} Var(\beta_{IV})-Var(\widehat\beta)=\sigma^2\bigg[\big(Z^tX\big)^{-1}\big(Z^tZ\big)\big(X^tZ\big)^{-1}-\big(X^tX\big)^{-1}\bigg] \end{equation*}\]

Then by using the following property:

\[\begin{equation*} A-B\geqslant 0\Leftrightarrow (B)^{-1}-(A^{-1})\geqslant0 \end{equation*}\]

The right hand side of the above formula will be (ignoring for now \(\sigma^2\)):

\[\begin{align*} X^tX-X^tZ\big(Z^tZ\big)^{-1}Z^tX&=X^t\bigg(I-Z\big(Z^tZ\big)^{-1}Z^t\bigg)X\\ &=X^tM_zX \end{align*}\]

As \(M_z\) is idempotent and symmetric matrix semidefinite shown in the previous chapter, the right hand side of the above formula will be:

\[\begin{equation*} X^tM_zX=X^tM^t_zM_zX=(M^t_zX)^t(M_zX)\geqslant0 \end{equation*}\]

Which is quadratic form, and hence:

\[\begin{equation*} Var(\beta_{IV})\geqslant Var(\widehat\beta) \end{equation*}\]

Unsurprisingly, this mathematical result was expected, because the IV estimation does not taking all the information contained in the regressors as the instruments drop part of the variations of the endogenous regressors related to the variation of the disturbances. Therefore, the cost of using the IV method to eliminate the bias caused by the non-orthogonality is having more uncertainty about the IV estimators (larger variances).

1- Let the matrix \(X\) has the following regressors \(X=(x_1,x_2,x_3,x_4,x_5)\).

Suppose that just the first three regressors \((x_1,x_2,x_3)\) are measured with errors (they do not satisfy the orthogonality condition), named appropriately endogenous, and the last ones are exogenous. So the proposed matrix of instruments should contain the exogenous variables as they satisfy the orthogonality condition, and other new variables as follows:\(Z=(z_1,z_2,z_3,x_4,x_5)\).

2- The exogenous variable values and their fitted values remain the same \((x_4=\widehat x_4,x_5=\widehat x_5)\), because they belong to the matrix of instruments \(Z\) on which they are regressed. That is they are perfectly fitted.

4.2.2.2 Tow-stage last squares (2SLS)

Now suppose that we have more instruments than regressors, that is \(Z\) is \((n\times m)\) matrix where \(m>k\). The matrix \((Z^tX)\) is \((m \times k)\) non invertible rectangular matrix, and hence the IV estimators can not be computed from the formula (4.20). However, using what we have discussed above about the two steps regression on which the IV method is based on, we regress first \(X\) on \(Z\) which gives the fitted values \(\widehat X=Z(Z^tZ)^{-1}Z^tX=H_ZX\), then in the next step, we regress \(y\) on \(\widehat X\) as follows (see 3.3 for the properties of \(H\)):

\[\begin{align} \beta_{IV}&=\big(\widehat X^t\widehat X\big)\widehat X^ty\notag\\ &=\big(X^tH_Z^tH_ZX\big)^{-1}X^tH^t_Zy\notag\\ &=\big(X^tH_ZX\big)^{-1}X^tH_Zy\notag\\ &=\Big(X^tZ\big(Z^tZ\big)^{-1}Z^tX\Big)^{-1}X^tZ\big(Z^tZ\big)^{-1}Z^ty \tag{4.24} \end{align}\]

\[\begin{align} Var(\beta_{IV})&=\big(\widehat X^t\widehat X\big)^{-1}\widehat X^t\varepsilon\varepsilon^t\widehat X\big(\widehat X^t\widehat X\big)^{-1}\notag\\ &=\sigma^2\big(\widehat X^t\widehat X\big)^{-1}\notag\\ &=\sigma^2\big(X^tH_ZX\big)^{-1} \tag{4.25} \end{align}\]

Obviously, this expression cannot be simplified since the matrix \(Z^tX\) is not square matrix (\(m>k\)). However, if \(m=k\) then the formula (4.24) reduces to \(\beta_{IV}=(Z^tX)^{-1}Z^ty\) obtained previously.

Alternatively, this result can also be obtained by the so called General least squares GLS that will be discussed later.

The results above hold if the assumptions of the fixed regressors and the orthogonality condition hold, assuming that all the other standard assumptions are satisfied, especially the assumptions of homoskedasticity and non autocorrelations. The IV estimation under heteroskedasticity and/or autocorrelation will be treated in more detail in the following chapter.

Using the hat matrix \(H\), the IV estimator variance in the formula (4.22) can be rewritten:

\[\begin{align*} Var(\beta_{IV})&=\frac{\sigma^2}{n}Q_{ZX}^{-1}Q_{ZZ}Q_{XZ}^{-1}\\ &=\frac{\sigma^2}{n}\Big(Q_{XZ}Q_{ZZ}^{-1}Q_{ZX}\Big)^{-1}\\ &=\frac{\sigma^2}{n}\Bigg(\frac{1}{n}X^tZ\Big(\frac{1}{n}Z^tZ\Big)^{-1}\frac{1}{n}Z^tX\Bigg)^{-1}\\ &=\sigma^2\Big(X^tH_zX\Big)^{-1} \end{align*}\]

4.2.2.3 Statistical tests for the the significance of the model

since the IV estimator is normally distributed as shown in (4.22), and the variance is estimated by \(\widehat\sigma^2=\frac{\sum e_{IV}^2}{n}\), then the individual significance of each parameter can be tested by the usual t-test, as discussed in the previous chapters, which follows asymptotically the normal distribution. The use of the F-test computed from the standard formula \(R^2=1-\frac{SSR}{SST}\), unlike t-student, might not remain appropriate since the \(R^2\) value can be negative, because the sum of the squared IV residuals \(SSR_{IV}\) can be larger than the total sum of squares \(SST\). Moreover, the use of the IV residuals \(SSR_{IV}\) to compare between the restricted model and the unrestricted one will no longer be useful since \(SSR^{UR}_{IV}\) can be larger than \(SSR^{R}_{IV}\) which might make F-test negative.

However, the F-test based on the comparison between restricted and unrestricted estimation discussed in the previous chapter can be applied for IV estimation but after a slight transformation that takes into account the above results.

Les us call the true model (3.44) \(y=X_I\beta_I+X_{II}\beta_{II}+\varepsilon\). The OLS estimated unrestricted model then will be \(y=X_I\widehat\beta_I+X_{II}\widehat\beta_{II}+e\), and the restrited one will be \(y=X_I\widehat\beta_R+e_R\), where \(\widehat\beta_{II}\) can be computed separately from the formula (3.47). Therefore, under the null hypothesis \(H0:\beta_{II}=0\), the estimator \(\widehat\beta_{II}\) is normally distributed with mean zero and variance matrix \(\sigma^2\big(X^t_{II}M_IX_{II}\big)^{-1}\), hence:

\[\begin{equation} \frac{\widehat\beta^t_{II}\big(X^t_{II}M_IX_{II}\big)\widehat\beta_{II}}{\sigma^2}\sim\chi(\rm g) \end{equation}\]

Where \(\rm g\) is the length of the vector \(\widehat\beta_{II}\). Dividing then that expression by \(\frac{e^te}{\sigma^2}\sim\chi(k-n)\), we get the following F-test:

\[\begin{equation*} F^c_{(\rm g, n-k)}=\frac{\widehat\beta^t_{II}\big(X^t_{II}M_IX_{II}\big)\widehat\beta_{II}\Big/\rm g}{e^te\Big/(n-k)} \end{equation*}\]

It was also shown that this test can be computed using the residuals that gives the same result:

\[\begin{equation*} F^c_{(\rm g, n-k)}=\frac{(e^t_Re_R-e^te)\Big/\rm g}{e^te\Big/(n-k)} \end{equation*}\]

Suppose now that \(H_Z=Z\big(Z^tZ\big)^{-1}Z^t\) can be rewritten in terms of another \((m\times n)\) matrix denoted \(\Gamma\) such that \(H_Z=\Gamma^t\Gamma\).

Let us premultiply our initial classical model by this matrix as follows:

\[\begin{equation} \Gamma y=\Gamma X \beta+\Gamma \varepsilon \tag{4.26} \end{equation}\]

Since the instruments satisfy the orthogonality condition, the same holds for these new transformed regressors \(\Gamma X\) with new errors vector \(\Gamma\varepsilon\), which means that the OLS method can be applied now on (4.26), provided that the transformed errors \(\Gamma \varepsilon\) satisfy the remaining classical assumptions. Since \(\varepsilon\rightarrow\rm N\big(0,\sigma^2I\big)\), those residuals then follow also the normal distribution \(\Gamma\varepsilon\rightarrow\rm N\big(0,\sigma^2\Gamma\Gamma^t\big)\), but in order to satisfy the required assumptions, the matrix \(\Gamma\) should have the property \(\Gamma\Gamma^t=I_m\). The question now is whether this matrix exists or not?

Using the properties of the hat matrix \(H\) (see 3.3), the mtrix \(H_Z\) will be:

\[\begin{equation*} \underset{(n\times n)}{H_Z}\underset{(n\times m)}{Z}=Z\big(Z^tZ\big)^{-1}Z^tZ=Z \end{equation*}\]

Since this matrix is positive semidefinite with \(rank(H_Z)=m\), it can be diagonalized as follows:

\[\begin{equation*} H_Z=VDV^t \end{equation*}\]

where D is a diagonal matrix that contains the eigenvalues of \(H_Z\), and \(V\) is orthogonal \(V^tV=I_m\). Substituting the last expression into the first one we get:

\[\begin{equation*} H_ZZ=VDV^tZ=Z \end{equation*}\]

Now, as \(Z\) has \(m\) columns, that equation requires that the matrix \(D\) should have \(m\) eigenvalues, each of which equals 1, and all the remaining ones \(m_2=n-m\) are zeroes such that:

\[\begin{equation*} H_Z=VDV^t=\begin{pmatrix}\underset{(n\times m)}{V_1} & \underset{(n\times m_2)}{V_2}\end{pmatrix}\begin{pmatrix}\underset{(m\times m)}{I} & \underset{(m\times m_2)}{0}\\\underset{(m_2\times m)}{0}&\underset{(m_2\times m_2)}{0}\end{pmatrix}\begin{pmatrix}\underset{(m\times n)}{V^t_1} \\ \underset{(m_2\times n)}{V^t_2}\end{pmatrix}=\underset{(n\times m)}{V_1}\underset{(m\times n)}{V^t_1} \end{equation*}\]

That is we take just the first \(m\) vectors from the matrix \(V\), and hence the required matrix \(\Gamma\) should be equal to \(\underset{(m\times n)}{V^t_1}\)

With that said, the estimation of the initial model by means of the IV method defined by:

\[\begin{equation} y=X\beta_{IV}+e^{IV} \tag{4.27} \end{equation}\]

is equivalent to the estimated transformed model by means of the OLS defined by:

\[\begin{equation} \Gamma y=\Gamma X\widehat\beta+e \tag{4.28} \end{equation}\]

That is if we premultiply (4.27) by the matrix \(\Gamma\) (\(\Gamma y=\Gamma X\beta_{IV}+\Gamma e^{IV}\)), we will have the same model with the same variables, and using the fact that the two estimators (IV) on (4.27) and OLS on (4.28) are equals, as discussed before, then we get \(e=\Gamma e^{IV}\). By the same, the IV and OLS residuals of the restricted model will be \(e_R=\Gamma e^{IV}_R\).

Combining the above results, F-test will be computed by the following:

\[\begin{align} F^c_{(\mathrm g, n-k)}&=\frac{\big(e^t_Re_R-e^te\big)\Big/\mathrm g}{e^te\Big/(n-k)}\notag\\ &=\frac{\Bigg[\Big(\Gamma e^{IV}_R\Big)^t\Big(\Gamma e^{IV}_R\Big)-\Big(\Gamma e^{IV}\Big)^t\Big(\Gamma e^{IV}\Big)\Bigg]\Big/\mathrm g}{\Big(\Gamma e^{IV}\Big)^t\Big(\Gamma e^{IV}\Big)\Big/(n-k)}\notag\\ &=\frac{\Bigg(e^{IV^t}_R\Gamma^t\Gamma e^{IV}_R-e^{IV^t}\Gamma^t\Gamma e^{IV}\Bigg)\Big/\mathrm g}{e^{IV^t}\Gamma^t\Gamma e^{IV}\Big/(n-k)}\notag\\ &=\frac{\Bigg(e^{IV^t}_RH_Z e^{IV}_R-e^{IV^t}H_Z e^{IV}\Bigg)\Big/\mathrm g}{e^{IV^t}H_Z e^{IV}\Big/(n-k)} \tag{4.29} \end{align}\]

4.2.2.4 Statistical tests of the IV estimation

As mentioned before, the IV estimators have larger standard errors than the OLS estimators, which means that their use thus must be strongly justified by the presence of the following requirements:

1. All or some of the regressors are endogenous.
2. The instruments should be sufficiently correlated to the endogenous regressors (instruments relevance).
3. The instruments should be orthogonal to the disturbance term (exogeneity of instruments).

these conditions can be checked by the following statistic tests.

Testing the endogeneity

As the OLS method implies orthogonality between regressors and regression residuals indicated by the normal equations \(X^te=0\), It will be thus impossible to use the residuals to test the endogeneity of the regressors. However, we can check the justification of the use of the IV method by simply comparing that estimation with the OLS estimation so that if they are significantly different, then the null hypothesis of no endogeneity should be rejected, and hence the use of IV is justified, and their estimators would be consistent (providing that these instruments used are relevant and exogenous). Otherwise, The OLS ones should be preferred due to their efficiency.

Since both estimators are asymptotically normally distributed, their difference is also asymptotically normally distributed. Therefore, under the null hypothesis, the limit distribution of the mean will be:

\[\begin{align*} plim\big(\beta_{IV}-\widehat\beta\big)&=plim(\beta_{IV})-plim(\widehat\beta)\\ &=\Bigg[\Big[\Big(\frac{1}{n}X^tZ\Big)\Big(\frac{1}{n}Z^tZ\Big)^{-1}\Big(\frac{1}{n}Z^tX\Big)\Big]^{-1}\Big[\Big(\frac{1}{n}X^tZ\Big)\Big(\frac{1}{n}Z^tZ\Big)^{-1}\overbrace{\Big(\frac{1}{n}Z^t\varepsilon\Big)}^{=0}\Bigg]\\ &-plim\Bigg[\Big(\frac{1}{n}X^tX\Big)^{-1}\overbrace{\Big(\frac{1}{n}X^t\varepsilon\Big)}^{=0}\Bigg]\\ &=0 \end{align*}\]

Not that under the alternative hypothesis (under endogeneity) \(plim(\widehat\beta)\not= \beta\).

Using the properties of \(H_Z\), the asymptotic variance of this difference is equal to:

\[\begin{align*} &E\Bigg[\Big(\beta_{IV}-\widehat\beta\Big)\Big(\beta_{IV}-\widehat\beta\Big)^t\Bigg]\\ &=\Bigg[\Big(X^tH_ZX\Big)^{-1}X^tH_Z-\Big(X^tX\Big)^{-1}X^t\Bigg]E\Big(\varepsilon\varepsilon^t\Big)\Bigg[H_ZX\Big(X^tH_ZX\Big)^{-1}-X\Big(X^tX\Big)^{-1}\Bigg]\\ &=\sigma^2\Bigg[\Big(X^tH_ZX\Big)^{-1}-\Big(X^tX\Big)^{-1}-\Big(X^tX\Big)^{-1}+\Big(X^tX\Big)^{-1}\Bigg]\\ &=\sigma^2\Bigg[\Big(X^tH_ZX\Big)^{-1}-\Big(X^tX\Big)^{-1}\Bigg]\\ &=Var\big(\beta_{IV}\big)-Var\big(\widehat\beta\big) \end{align*}\]

Replacing the expectation operator by the plim operator, the asymptotic variance will be:

\[\begin{equation*} var\Big(\beta_{IV}-\widehat\beta\Big)=\frac{\sigma^2}{n}plim\Big(\frac{1}{n}X^tH_ZX\Big)^{-1}-\frac{\sigma^2}{n}plim\Big(\frac{1}{n}X^tX\Big)^{-1} \end{equation*}\]

As a result, the distribution of the variable of the difference under the null hypothesis is \(\Big(\beta_{IV}-\widehat\beta\Big)\overset{Asy}\rightarrow\mathrm N\Bigg(0,Var\big(\beta_{IV}\big)-Var\big(\widehat\beta\big)\Bigg)\). By standardizing and squaring that variable, we get a chi-square variable:

\[\begin{equation} \Big(\beta_{IV}-\widehat\beta\Big)^t\Bigg[Var\big(\beta_{IV}\big)-Var\big(\widehat\beta\big)\Bigg]^{-1}\Big(\beta_{IV}-\widehat\beta\Big)\overset{Asy}\rightarrow\chi_{(k)} \tag{4.30} \end{equation}\]

A Large value for this statistic leads to the rejection of the null hypothesis meaning that the use of the IV estimation est is justified. However, the variances in (4.30) are unknown, but they can be replaced by their usual estimators \(\widehat{Var}\big(\beta_{IV}\big)\) and \(\widehat{Var}\big(\widehat\beta\big)\).

Alternatively, there exists a test called Wu Haussman test that is easier to apply as it uses only the regression residuals resulted from regressing the original regressors on the instruments. To understand its implementation, let us call the classical simple linear regression model:

\[\begin{equation} y_t=\beta_1+\beta_2x_t+\varepsilon_t \tag{4.31} \end{equation}\]

Then suppose that the regression of \(x_t\) on \(z_t\) is defined by:

\[\begin{equation} x_t=\widehat\alpha_1+\widehat\alpha_2z_t+\widehat\mu_t \tag{4.32} \end{equation}\]

Since \(z_t\) is assumed to be uncorrelated with \(\varepsilon_t\), the expression (4.32) is supposed to isolate the part of \(x_t\) (that is correlated with \(\varepsilon_t\)) in the error term \(\widehat\mu_t\). In other words, the only way for \(x_t\) to be correlated with \(\varepsilon_t\) is through \(\widehat\mu_t\). that means that the original error term must be explained by \(\widehat\mu_t\) such as:

\[\begin{equation} \varepsilon_t=\gamma\widehat\mu_t+v_t \tag{4.33} \end{equation}\]

Substituting this expression into (4.31), we get:

\[\begin{equation} y_t=\beta_1+\beta_2x_t+\gamma\widehat\mu_t+v_t \tag{4.34} \end{equation}\]

Now, this equation is ready to be estimated by OLS since the vector of errors \(v_t\) is orthogonal to both regressors \(x_t\) and \(\widehat\mu_t\). Then we can readily use the usual individual t-test to test whether \(\widehat\mu_t\) is significant or not. If it is, then this part of \(x_t\) (\(\widehat\mu_t\)) is indeed correlated with \(\varepsilon_t\) from (4.33), and hence the orthogonality condition is not satisfied.

The Haussman test can be straightforwardly extended to the case of multiple linear models. If a model has \(m\) potential endogenous regressors among the \(k\) regressors, then the test will be implemented as follows:

1. Searching for at least \(m\) instruments.
2. Regressing each potential endogenous variable on these instruments by OLS (estimating the reduced form) and computing the residuals.
3. Including all these residuals as new regressors in the initial model.
4. Test the joint distribution by \(W\) test:

\[\begin{equation*} W=\frac{\big(SSR_{\text{with residuals}}-SSR\big)\Big/\mathrm g}{SSR_{\text{with residuals}}\Big/(n-k-\mathrm g)}\overset{Asy}\rightarrow F(n-k-\mathrm g,\mathrm g) \end{equation*}\]

Where \(\mathrm g\) is the number of the suspected regressors.

Let us take an example of a model with 3 regressors \(y_t=\beta_1+\beta_2x_{2t}+\beta_3x_{3t}+\beta_4x_{4t}+\varepsilon_t\), and suppose that two regressors \(x_{3t}\) and \(x_{4t}\) are suspected to be endogenous. Thus, we should find at least two relevant instruments, say three ones, \(x_{2t}\) and two new ones \(z_{1t}\) and \(z_{2t}\). Then, performing (by OLS) the first stage regression and retrieve the computed residuals:

\[\begin{equation*} \begin{cases} x_{3t}=\widehat\gamma^1_1+\widehat\gamma^1_2x_{2t}+\widehat\gamma^1_3z_{1t}+\widehat\gamma^1_4z_{2t}+\widehat\mu_{1t}\\ x_{4t}=\widehat\gamma^2_1+\widehat\gamma^2_2x_{2t}+\widehat\gamma^2_3z_{1t}+\widehat\gamma^2_4z_{2t}+\widehat\mu_{1t} \end{cases} \end{equation*}\]

Then testing the joint significant of \((\widehat\mu_{1t},\widehat\mu_{2t}\), after having been included in the initial model, by the usual W-test (\(F(n-5, 2)\)):

\[\begin{equation*} y_t=\widehat\beta_1+\widehat\beta_2x_{2t}+\widehat\beta_3x_{3t}+\widehat\beta_4x_{4t}+\widehat\beta_5\widehat\mu_{1t}+\widehat\beta_6\widehat\mu_{2t}+e_t \end{equation*}\]

A large value of this test leads to rejecting the null hypothesis (\(H0:\widehat\mu_{1t}=\widehat\mu_{2t}=0\)), and hence the use of the IV estimation will be justified.

Alternatively, W-test can be replaced by LM-test as follows:

1. Estimate (by OLS) the original model \(y_t=\beta_1+\beta_2x_{2t}+\beta_3x_{3t}+\beta_4x_{4t}+\varepsilon_t\) and obtaining the the residuals \(e_t\).
2. Estimate the auxiliary equation \(e_t=\widehat\gamma_1 +\widehat\gamma_2x_{2t}+\widehat\gamma_3x_{3t}+\widehat\gamma_4x_{4t}+\widehat\gamma_5\widehat\mu_{1t}+\widehat\gamma_6\widehat\mu_{2t}+v_t\) then compute \(R^2\).
3. Compute \(LM=nR^2\). since this test follows asymptotically the \(\chi_{(\mathrm g)}\) under the null hypothesis, a large value leads to reject that hypothesis.

Testing the relevance of instruments

The instruments are called relevant if they are sufficiently correlated to the endogenous regressors. Thus, the easiest way to check whether they are relevant or not is by using the first stage regression in the IV estimation, in which we regress each potential endogenous regressor on the instruments. Then, using the usual F-test for each regression to test the joint significance of these instruments so that larger values lead to the rejection of the null hypothesis. That means that those instruments are not relevant.

However, it is not only required for F-test to be significant, but also should have larger values. With small values, the instruments are called weak instruments, which means that the IV estimation still give consistent estimators but with very large standard errors.That is:

\[\begin{align*} Var(\beta_{IV})&=\sigma^2\big(Z^tX\big)^{-1}(Z^tZ\big)(X^tZ\big)^{-1}\\ &=\sigma^2\big(X^tX\big)^{-1}\big(X^tX\big)\big(Z^tX\big)^{-1}(Z^tZ\big)(X^tZ\big)^{-1}\\ &=Var(\widehat\beta)\big(X^tX\big)\big(Z^tX\big)^{-1}(Z^tZ\big)(X^tZ\big)^{-1}\\ &=Var(\widehat\beta)\Big[X^tZ\big(Z^tZ\big)^{-1}Z^tX\big(X^tX\big)^{-1}\Big]^{-1}\\ &=Var(\widehat\beta)\Big[\big(X^tZ\big)\widehat\gamma\big(X^tX\big)^{-1}\Big]^{-1} \end{align*}\]

Since \(\widehat\gamma\) is the matrix of the estimates of the reduced form (regression of each vector of \(X\) on \(Z\)), so a small significance (weak instruments) means that the components of \(\widehat\gamma\) will be closer to zero (without losing the significance), and hence the variances of IV estimators will be a lot larger than usual (with relevant instruments). A small inconsistency with small uncertainty sometimes is preferred to consistency with large uncertainty.

As a rule of thumb, many econometricians suggest that the F-test value should be larger than 10 to treat the instruments as relevant (not weak) (Staiger and Stock 1997).

Testing the validity of the instruments:

An instrument is said to be valid if it satisfies the orthogonality condition. Can this property be checked by testing the correlation between these instruments and the error resulted from the IV regression \(e_{IV}\)? The answer is yes, unlike the OLS regression that implies the orthogonality condition between the residuals and the regressors.

The Sargan test based on this idea as follows:

1. Regress the errors \(e_{IV}\) on the instruments \(z\) (auxiliary regression) \(e_{IV}=Z\widehat\gamma+\zeta\).
2. Check the significance by the F-test or by LM-test defined by \(LM=nR^2\overset{Asy}\rightarrow\chi_{(m-k)}\), where \(R^2\) is calculated from the auxiliary regression. Larger values leads to the rejection of the null hypothesis of exogeneity of instruments.

Note that the degrees of freedom of \(\chi\) is \((m-k)\), which is equal to zero if the number of the instruments is equal to the number of the regressors. That is, this test cannot be applied in the just identified case. It is only applied if the model is over identified, that is why this test is also called over identified test.

Example 4.1 To illustrate the IV method, we will use the data SchoolingReturns provided in the R package ivreg. For simplifiacation, we use only the first eight variables as regressors to explain the variable wage.

In R:

library("ivreg")
data("SchoolingReturns", package = "ivreg")
data_iv <- SchoolingReturns[,1:8]
head(data)

Table 4.1: SchoolingReturns

wage	education	experience	ethnicity	smsa	south	age	nearcollege
548	7	16	afam	yes	no	29	no
481	12	9	other	yes	no	27	no
721	12	16	other	yes	no	34	no
250	11	10	other	yes	no	27	yes
729	12	16	other	yes	no	34	yes
500	12	8	other	yes	no	26	yes

This data has 3010 rows and 22 columns. The variables that we will use are the following:

• wage: Raw wages in 1976 (in cents per hour).
• education: Education in 1976 (in years).
• experience: Years of labor market experience.
• ethnicity: Factor indicating ethnicity. Is the individual African-American (“afam”) or not (“other”)?
• smsa : Factor. Does the individual reside in a SMSA (standard metropolitan statistical area) in 1976?
• south: Factor. Does the individual reside in the South in 1976?
• age: Age in 1976 (in years).
• nearcollege: Factor. Did the individual grow up near a 4-year college?

This data was introduced and studied by(Card.D 1993). As he suspected that education and experience are endogenous variables, he argued that proximity to a college nearcollege can be a good instrument for education, and age for experience.

First let us fit the model with OLS method:

model_no_iv_R <- lm(log(wage) ~ education+experience+ethnicity+smsa+south, data = data_iv)

# display the summary
result_no_iv_R <- tidy(model_no_iv_R, conf.int = TRUE)
result_no_iv_R

Table 4.2: The estimates with no IV in R

term	estimate	std.error	statistic	p.value	conf.low	conf.high
(Intercept)	4.9133308	0.0631212	77.839680	0	4.7895657	5.0370959
education	0.0738070	0.0035336	20.887072	0	0.0668784	0.0807356
experience	0.0393134	0.0021955	17.906109	0	0.0350085	0.0436183
ethnicityafam	-0.1882225	0.0177678	-10.593489	0	-0.2230607	-0.1533843
smsayes	0.1647411	0.0156919	10.498506	0	0.1339732	0.1955090
southyes	-0.1290528	0.0152285	-8.474422	0	-0.1589122	-0.0991935

We see that a year of education increases the log(wages) by nearly 7%, and the experience by 4%. However, those who are african-american or live in the south have lower wages.

Now let us use the IV regression where the exogenous variables ethnicity, smsa, and south will be added as instruments with nearcollege and squared age.

model_iv_R <- ivreg(log(wage) ~ education+experience+ethnicity+smsa+south | ethnicity+smsa+south+nearcollege+poly(age,2), data = data_iv)

# display the summary
summary(model_iv_R)

[out] 
[out] Call:
[out] ivreg(formula = log(wage) ~ education + experience + ethnicity + 
[out]     smsa + south | ethnicity + smsa + south + nearcollege + poly(age, 
[out]     2), data = data_iv)
[out] 
[out] Residuals:
[out]      Min       1Q   Median       3Q      Max 
[out] -1.92396 -0.26644  0.02182  0.27788  1.36688 
[out] 
[out] Coefficients:
[out]                Estimate Std. Error t value Pr(>|t|)    
[out] (Intercept)    3.826040   0.484136   7.903 3.79e-15 ***
[out] education      0.155740   0.036302   4.290 1.84e-05 ***
[out] experience     0.040596   0.002558  15.869  < 2e-16 ***
[out] ethnicityafam -0.069513   0.056038  -1.240 0.214906    
[out] smsayes        0.087650   0.038346   2.286 0.022338 *  
[out] southyes      -0.088287   0.024935  -3.541 0.000405 ***
[out] 
[out] Diagnostic tests:
[out]                                df1  df2 statistic  p-value    
[out] Weak instruments (education)     3 3003     8.008 2.58e-05 ***
[out] Weak instruments (experience)    3 3003  1612.707  < 2e-16 ***
[out] Wu-Hausman                       1 3003     6.672  0.00984 ** 
[out] Sargan                           1   NA     0.312  0.57624    
[out] ---
[out] Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[out] 
[out] Residual standard error: 0.4287 on 3004 degrees of freedom
[out] Multiple R-Squared: 0.06857,  Adjusted R-squared: 0.06702 
[out] Wald test: 157.1 on 5 and 3004 DF,  p-value: < 2.2e-16

As we see, the function ivreg displayed the results with the above-discussed tests. The estimate of education now is much larger (15.57%). Since the p_values in the Diagnostic test table are very small (2.58e-05, 2e-16), the instruments are relevant (not weak). The Wu-Hausman has a small p-value indicating the justification of using the IV regression. For testing the validity of the instruments, The Sargan test has a very large p-value (0.57624), indicating that they are valid instruments (exogenous).

In Python:

First let us move the data from R to python.

if not 'data_iv_py' in globals():
  data_iv_py = r.data_iv

For the IV regression in python we will import the python package linearmodels. To run the IV regression we use the function IV2SLS that has four arguments in order: dependent variable, the exogenous variables, the endogenous, and the last one is the instruments. We should also mention that this function by default uses robust standard errors, which are not yet covered. To use the classical standard errors we set the argument cov_type inside the function fit to unadjusted

# we use the vectorized fucntion log in numpy
import numpy as np
from linearmodels.iv import IV2SLS

# add consatant
data_iv_py['constant']=1

# creat log wage
data_iv_py['lwage'] = np.log(data_iv_py.wage)

# create the squared age
data_iv_py['age2']=np.exp2(data_iv_py.age)

#define the arguments
dep = 'lwage'
exog = ['constant', 'ethnicity' , 'smsa' , 'south']
endog = ['education', 'experience']
instr = ['nearcollege', 'age', 'age2']

# specify and fit the model in one go
model_iv_py=IV2SLS(data_iv_py[dep], data_iv_py[exog], data_iv_py[endog], data_iv_py[instr]).fit(cov_type='unadjusted')

# display the results
model_iv_py

IV-2SLS Estimation Summary

Dep. Variable:	lwage	R-squared:	-0.0656
Estimator:	IV-2SLS	Adj. R-squared:	-0.0674
No. Observations:	3010	F-statistic:	688.54
Date:	Thu, Sept 17 2026	P-value (F-stat)	0.0000
Time:	11:13:06	Distribution:	chi2(5)
Cov. Estimator:	unadjusted

Parameter Estimates

	Parameter	Std. Err.	T-stat	P-value	Lower CI	Upper CI
constant	3.4853	0.5340	6.5264	0.0000	2.4386	4.5320
ethnicity.other	0.0363	0.0683	0.5312	0.5953	-0.0976	0.1701
smsa.yes	0.0661	0.0462	1.4303	0.1526	-0.0245	0.1566
south.yes	-0.0769	0.0289	-2.6571	0.0079	-0.1336	-0.0202
education	0.1787	0.0449	3.9774	0.0001	0.0906	0.2667
experience	0.0410	0.0028	14.858	0.0000	0.0356	0.0464

Endogenous: education, experience
Instruments: nearcollege.yes, age, age2
Unadjusted Covariance (Homoskedastic)
Debiased: False
id: 0x20e5e924a40

Unlike the R function ivreg that outputs all the tests. Here we have to explicitly call those tests.

We can test if an endogenous variable is actually endogenous or not using the wu_hausman test as follows:

# for education
model_iv_py.wu_hausman(['education'])

[out] Wu-Hausman test of exogeneity
[out] H0: Variables education are exogenous
[out] Statistic: 8.8778
[out] P-value: 0.0029
[out] Distributed: F(1,3003)
[out] WaldTestStatistic, id: 0x20e5e8136b0

The test has rejected the null hypothesis so that the use of the IV regression is justified

We can also get the sargan test as follows:

model_iv_py.sargan

[out] Sargan's test of overidentification
[out] H0: The model is not overidentified.
[out] Statistic: 1.4879
[out] P-value: 0.2225
[out] Distributed: chi2(1)
[out] WaldTestStatistic, id: 0x20e5e927f20

4.2.3 Simultaneous equations

Another cause of non orthogonality condition is when the relation between the dependent variable and a particular regressor is reciprocal, which means they depend on each other. For instance, if the causal relationship between two variables \(y\) and \(x\) is bidirectional, then the classical simple linear model \(y_t=\beta_1+\beta_2x_t+\varepsilon_{1t}\) characterizes the causal effect of \(x_t\) on \(y_t\). To represent the other direction however, we should add the equation \(x_t=\alpha_1+\alpha_2y_t+\varepsilon_{2t}\). Grouping these two equations together forms a system called simultaneous equations model SEM as follows:

\[\begin{equation} \begin{cases} y_t=\beta_1+\beta_2x_t+\varepsilon_{1t}\\ x_t=\alpha_1+\alpha_2y_t+\varepsilon_{2t} \end{cases} \tag{4.35} \end{equation}\]

It is clear that \(x_t\) is endogenous in the first equation, because it depends on \(y_t\) from the second equation. Therefore, the orthogonality in the first model is violated, and hence the OLS estimators will no longer be consistent.

Let us investigate the general form of a SEM system with \(\mathrm g\) endogenous variables and \(k\) exogenous ones. So there will be \(\mathrm g\) equations as follows:

\[\begin{equation} \begin{cases} a_{11}y_{1t}+a_{12}y_{2t}+...+a_{1\mathrm g}y_{\mathrm gt}+b_{11}x_{1t}+b_{12}x_{2t}+...+b_{1k}x_{kt}=\varepsilon_{1t}\\ a_{21}y_{1t}+a_{22}y_{2t}+...+a_{2\mathrm g}y_{\mathrm gt}+b_{21}x_{1t}+b_{22}x_{2t}+...+b_{2k}x_{kt}=\varepsilon_{2t}\\ ......................................................\\ a_{\mathrm g1}y_{1t}+a_{\mathrm g2}y_{2t}+...+a_{\mathrm g\mathrm g}y_{\mathrm gt}+b_{\mathrm g1}x_{1t}+b_{\mathrm g2}x_{2t}+...+b_{\mathrm gk}x_{kt}=\varepsilon_{\mathrm g t} \end{cases} \tag{4.36} \end{equation}\]

This is system known as structural form of SEM.

To account for the constant term we set the first exogenous variable to one \(x_{1t}=1\). The above expression can be rewritten in matrix form as:

\[\begin{equation*} \begin{pmatrix} a_{11}&a_{12}&..&a_{1\mathrm g}\\ a_{21}&a_{22}&..&a_{2\mathrm g}\\ ..&..&..&..\\ a_{\mathrm g1}&a_{\mathrm g2}&..&a_{\mathrm g\mathrm g} \end{pmatrix} \begin{pmatrix} y_{1t}\\ y_{2t}\\ ..\\ y_{\mathrm g t} \end{pmatrix}+ \begin{pmatrix} b_{21}&b_{22}&..&b_{2k}\\ b_{11}&b_{12}&..&b_{1k}\\ ..&..&..&..\\ b_{\mathrm g1}&b_{\mathrm g2}&..&b_{\mathrm gk} \end{pmatrix} \begin{pmatrix} x_{1t}\\ x_{2t}\\ ..\\ x_{k t} \end{pmatrix}= \begin{pmatrix} \varepsilon_{1t}\\ \varepsilon_{2t}\\ ..\\ \varepsilon_{\mathrm g t} \end{pmatrix} \end{equation*}\]

More compactly:

\[\begin{equation} \underset{(\mathrm g,\mathrm g)}A\underset{(\mathrm g,1)}{Y_t}+\underset{(\mathrm g,k)}B\underset{(k,1)}{X_t}=\underset{(\mathrm g,1)}{\varepsilon_t} \tag{4.37} \end{equation}\]

We should point out that this matrix form does not embody all the equations over time as what is done with the general model in the previous chapter. that expression characterizes the system at the specific time \(t\) (as indicated by the subscript \(t\)).

The system (4.36) can be written in an another matrix form different than that in (4.37) as follows (used by well known researchers such as Green and Amemiya):

\[\begin{equation*} \underset{(1,\mathrm g)}{Y^t_t}\underset{(\mathrm g,\mathrm g)}{P}+\underset{(1,k)}{X^t_t}\underset{(k,\mathrm g)}{\Gamma}=\underset{(1,\mathrm g)}{\varepsilon^t_t} \end{equation*}\]

It seems that the two forms are different, but in fact they are the same because:

\[\begin{align*} AY_t+BX_t=\varepsilon_t&\implies\big(AY_t+BX_t\big)^t=\varepsilon^t_t\\ &\implies Y^t_tA^t+X^t_tB^t=\varepsilon^t_t=Y^t_tP+X^t_t\Gamma \end{align*}\]

That is \(B^t=\Gamma\) and \(A^t=P\)

To allow each endogenous variable be expressed only in terms of the exogenous ones, the matrix \(A\) must be non singular such that:

\[\begin{equation} Y_t=-A^{-1}BX_t+A^{-1}\varepsilon_t=\Gamma X_t+v_t \tag{4.38} \end{equation}\]

This system is called reduced form.

With this form, each endogenous variable is explained in terms of the exogenous variables, which are now orthogonal to the errors \(\varepsilon_t\), and also to the new errors \(v_t\). Therefore, the OLS method can be applied to (4.38) and give consistent estimators provided that all the classical assumptions are satisfied (\(A2\), \(A3\), \(A4\)) for each equation. That is:

\[\begin{equation*} \begin{cases} E(v_{it})=E(\varepsilon_{it})=0\quad t\not=t^{'}\\ E(v_{it}v_{it^{'}})=E(\varepsilon_{it}\varepsilon_{it^{'}})=0\\ E(\varepsilon_{it}^2)=\sigma^2_{\varepsilon,i}=\sigma^2_i\\ E(v_{it}^2)=\sigma^2_{v,i} \end{cases} \end{equation*}\]

However, nothing has been said about any dependence between the errors of different equations \(E(\varepsilon_{it}\varepsilon_{jt})\quad , \quad i;j=1,2,..,\mathrm g\) because they do not affect the consistency of the OLS estimators (but they may affect the efficiency). Thus, the unknown contemporaneous covariance matrix for \(\varepsilon_t\) denoted by:

\[\begin{equation*} \underset{\mathrm g\times \mathrm g}\Sigma= \begin{pmatrix} \sigma_{11}&\sigma_{12}&..&\sigma_{1\mathrm g}\\ \sigma_{21}&\sigma_{22}&..&\sigma_{2\mathrm g}\\ ..&..&..&..\\ \sigma_{\mathrm g 1}&\sigma_{\mathrm g 2}&..&\sigma_{\mathrm g\mathrm g} \end{pmatrix} \end{equation*}\]

In general this matrix is not diagonal matrix, and the same holds for the covariance matrix of \(v_t\) derived from \(\Sigma\) by \(\Omega_v=A^{-1}\Sigma (A^{-1})^t\). With that done, the estimated model of (4.38) will be:

\[\begin{equation} Y_t=\widehat\Gamma X_t+e_t \tag{4.39} \end{equation}\]

4.2.3.1 Identification

We know that \(-A^{-1}B\) is \((\mathrm g\times k)\) matrix, so \(\widehat\Gamma\) contains \((\mathrm g\times k)\) elements that will be used to derive the estimates of each parameter in A and B. But as the total number of the parameters of interest is \((\mathrm g\times\mathrm g)+(\mathrm g\times k)\), which is greater than that of the estimated elements of \(\widehat\Gamma\), then, the derivation of the original parameters from this matrix will be impossible unless some restrictions on parameters are provided. If so that model will be called overidentified, or at least justidentified, and hence it can be estimated.

The estimation problem in this type of model appears only when we want to compute the parameters of the original matrices \(A\) and \(B\) using the model estimates in \(\widehat\Gamma\). That matrix has \(\mathrm g\times k\) elements characterizing the number of relations (equations) in the system. Whereas, \(A\) and \(B\) together have \(\mathrm g^2+\mathrm g k\) unknown elements to compute. Therefore, to get a solution for a system of \(\mathrm g\times k\) equations and \(\mathrm g^2+\mathrm g k\) variables, we should set at least \(\mathrm g^2\) constraints on the variables.

Example 4.2 Let us take the popular macroeconomic system with three endogenous variables (Consumption \(C\), Income \(Y\), and Investment \(I\)), and one exogenous variable (Interest rate \(r\)) as follows:

\[\begin{equation} \begin{cases} C_t=a_1+a_2Y_t+\varepsilon_{1t}\\ I_t=b_1+b_2r_t+\varepsilon_{2t}\\ Y_t=C_t+I_t \end{cases} \tag{4.40} \end{equation}\]

It can be rewritten in the form of (4.36) as follows:

\[\begin{equation} \begin{cases} C_t+0I_t-a_2Y_t+0r_t-a_1u_t=\varepsilon_{1t}\\ 0C_t+I_t+0Y_t-b_2r_t-b_1u_t=\varepsilon_{2t}\\ -C_t-I_t+Y_t+0r_t+0u_t=0 \end{cases} \tag{4.41} \end{equation}\]

Or in matrix form:

\[\begin{equation} \begin{pmatrix} 1&0&-a_2\\ 0&1&0\\ -1&-1&1 \end{pmatrix} \begin{pmatrix} C_t\\ I_t\\ Y_t \end{pmatrix}+ \begin{pmatrix} 0&-a_1\\ -b_2&-b_1\\ 0&0 \end{pmatrix} \begin{pmatrix} r_t\\ u_t \end{pmatrix}= \begin{pmatrix} \varepsilon_{1t}\\ \varepsilon_{2t}\\ 0 \end{pmatrix} \tag{4.42} \end{equation}\]

Where \(u_t=1\) stands for the constant term.

The matrices of the structural form (4.36) are:

\[\begin{equation*} A=\begin{pmatrix} 1&0&-a_2\\ 0&1&0\\ -1&-1&1 \end{pmatrix}\quad , \quad B=\begin{pmatrix} 0&-a_1\\ -b_2&-b_1\\ 0&0 \end{pmatrix} \end{equation*}\]

To get the reduced form, the matrix \(A\) must be invertible, which is the case if \(a_2\not=1\). Consequently, the reduced form will be:

\[\begin{equation} \begin{pmatrix} C_t\\ I_t\\ Y_t \end{pmatrix}= \begin{pmatrix} \frac{1}{1-a_2}&\frac{-a_2}{1-a_2}&\frac{-a_2}{1-a_2}\\ 0&1&0\\ \frac{-1}{1-a_2}&\frac{1}{1-a_2}&\frac{1}{1-a_2} \end{pmatrix} \begin{pmatrix} 0&-a_1\\ -b_2&-b_1\\ 0&0 \end{pmatrix} \begin{pmatrix} r_t\\ u_t \end{pmatrix}+ \begin{pmatrix} \frac{1}{1-a_2}&\frac{-a_2}{1-a_2}&\frac{-a_2}{1-a_2}\\ 0&1&0\\ \frac{-1}{1-a_2}&\frac{1}{1-a_2}&\frac{1}{1-a_2} \end{pmatrix} \begin{pmatrix} \varepsilon_{1t}\\ \varepsilon_{2t}\\ 0 \end{pmatrix} \tag{4.43} \end{equation}\]

Simplifying by solving the first member of the right hand side, we get:

\[\begin{equation} \begin{pmatrix} C_t\\ I_t\\ Y_t \end{pmatrix}= \begin{pmatrix} \frac{a_2b_2}{1-a_2}&\frac{a_2b_1-a_1}{1-a_2}\\ -b_2&-b_1\\ \frac{-b_2}{1-a_2}&\frac{a_1-b_2}{1-a_2} \end{pmatrix} \begin{pmatrix} r_t\\ u_t \end{pmatrix}+ \begin{pmatrix} v_{1t}\\ v_{2t}\\ v_{3t} \end{pmatrix} \tag{4.44} \end{equation}\]

Where

\[\begin{equation*} \begin{pmatrix} v_{1t}\\ v_{2t}\\ v_{3t} \end{pmatrix}= \begin{pmatrix} \frac{1}{1-a_2}&\frac{-a_2}{1-a_2}&\frac{-a_2}{1-a_2}\\ 0&1&0\\ \frac{-1}{1-a_2}&\frac{1}{1-a_2}&\frac{1}{1-a_2} \end{pmatrix} \begin{pmatrix} \varepsilon_{1t}\\ \varepsilon_{2t}\\ 0 \end{pmatrix} \end{equation*}\]

Since \(\varepsilon_t\approx IID(0,\Omega)\), then \(v_t\approx IID(0,A^{-1}\Omega(A^{-1})^t)\).

We see that each equation has two parameters to estimate (one for the exogenous variable and the constant term), and since the orthogonality condition is now satisfied, the OLS method can be applied to each equation separately.

Let us suppose that the estimated equations of (4.44) are:

\[\begin{equation} \begin{pmatrix} C_t\\ I_t\\ Y_t \end{pmatrix}= \begin{pmatrix} \widehat\alpha_1&\widehat\alpha_2\\ \widehat\beta_1&\widehat\beta_2\\ \widehat\gamma_1&\widehat\gamma_2 \end{pmatrix} \begin{pmatrix} r_t\\ u_t \end{pmatrix}+ \begin{pmatrix} e_{1t}\\ e_{2t}\\ e_{3t} \end{pmatrix} \tag{4.45} \end{equation}\]

Since the reduced form estimates are determined and expressed in terms of the original structural form parameters from (4.44), we expect to deduce a unique value for each one (just-identified). However, as the number of the reduced form estimates might differ from that of the original parameters, then the computation can give several solutions for the same parameter (over-identified), or even no solution for some parameters (under-identified ).

Therefore, each equation of the SEM faces three cases of identification:

1. Under-identified: the number of the restrictions is not enough to deduce all its original parameter estimates.
2. Just-identified: the number of the restrictions is the minimum required to deduce all the original parameter estimates.
3. Over-identified: the number of the restriction is greater than the minimum required.

That for each equation. The identification for the whole system depends on the following.

1. The model is under-identified if at least one of the equations is under-identified.
2. The model is just-identified if all the equations are just identified.
3. The model is over-identified if all the equations are either just-identified or over-identified.

As the identification based on the required number of restrictions, the investigator should look for at least the minimum number required using the theory behind the problem under study or extra information. There exist many types of restrictions. The most important ones are called zero restriction (or exclusion restriction), which are met when a particular variable (either endogenous or exogenous) is excluded from an equation (its coefficient equals zero). Other restrictions, called linearity combination restriction, can exist when a linear combination of some equations has the same variables as a particular equation of the SEM.

Now let us examine how to add restriction to the SEM. We start by transforming our structural model (4.37) as follows:

\[\begin{equation} \underset{\mathrm g\times(\mathrm g+k)}\Pi W_t=\big(A,B\big) \begin{pmatrix} Y_t\\ X_t \end{pmatrix}=\varepsilon_t \tag{4.46} \end{equation}\]

Thus, the \(i^{th}\) equation can be extracted from that expression as follows:

\[\begin{equation} \pi_iW_t=\varepsilon_{it} \tag{4.47} \end{equation}\]

Where \(\pi_i\) is the \(i^{th}\) row of \(\Pi\).

Now, as we did with W-test discussed in the previous chapter, the restrictions of the \(i^{th}\) equation should satisfy the following:

\[\begin{equation} \pi_i\underset{(\mathrm g+k)\times q}R=0 \tag{4.48} \end{equation}\]

Where \(q\) is the number of restrictions in the \(i^{th}\) equation.

Combine this last expression with that deduced from (4.38) that relates the reduced form parameters to those of the structural form \[\begin{equation*} A\Gamma + B=\big(A,B\big) \begin{pmatrix} \Gamma\\ I_k \end{pmatrix}=\Pi\Phi=0 \end{equation*}\]

we get:

\[\begin{equation} \pi_i\big(\Phi,R\big)=0 \tag{4.49} \end{equation}\]

Where \(\big(\Phi,R\big)\) is \((\mathrm g+k)\times(k+q)\) matrix. Eventually, this expression represents a system that has \((k+q)\) equations and \((\mathrm g+k)\) variables. To get thus a solution for the \(i^{th}\) equation (and hence satisfy the identification condition), the rank of that matrix must be at least equal to \(\mathrm g+k-1\):

\[\begin{equation} Rank\big(\Phi,R\big)\geqslant\mathrm g+k-1 \tag{4.50} \end{equation}\]

This is a necessary and sufficient condition called the rank condition.

Consequently, this condition can be applied to check the identification of all the model equations. In other words, an equation of the model is identified, if the matrix that contains all the variable coefficients (endogenous and exogenous) excluded from that equation has at least one determinant of order \((\mathrm g-1)\) equal to zero.

However, this condition is too critical to be satisfied, especially when the number of variables is larger. Fortunately, there exists another condition named the order condition, but with the inconvenience that it is only necessary and not sufficient. It is expressed by:

\[\begin{align} &k+q\geqslant\mathrm g+k-1\notag\\ &q\geqslant\mathrm g-1 \tag{4.51} \end{align}\]

That means, for a given equation to be satisfied, the number of restrictions must at least be equal to the number of equations minus 1.

Since the type of restriction that we face most in real-world problems is the exclusion restriction, practitioners set a rule of thumb to get the identification using the order condition by the following:

Let \(\mathrm g^{'}\) and \(k^{'}\) be the number of the excluded endogenous variables and the excluded exogenous ones respectively so that the number of restrictions is \(q=\mathrm g^{'}+k^{'}\), then the order condition will be:

\[\begin{equation} \mathrm g^{'}+k^{'} \geqslant\mathrm g-1 \tag{4.52} \end{equation}\]

The equation under study thus is identified if the number of all the excluded variables (endogenous and exogenous) is at least equal to the total number of the system equations minus one. So, the equation is over-identified if \(\mathrm g^{'}+k^{'} \geqslant\mathrm g-1\) and just-identified if \(\mathrm g^{'}+k^{'} =\mathrm g-1\). However, as this condition is necessary but not sufficient, we should always check the rank condition, when it is possible, that should solve the identification with certainty.

Let us call our previous example. As the reduced form parameters have been estimated, we expect to deduce the structural form estimates by solving the following system:

\[\begin{equation} \begin{cases} \widehat\alpha_1=\frac{\widehat a_2\widehat b_2}{1-\widehat a_2}\quad , \quad \widehat\alpha_2=\frac{\widehat a_2\widehat b_1-\widehat a_1}{1-\widehat a_2}\\ \widehat\beta_1=-\widehat b_2\quad , \quad \widehat\beta_2=-\widehat b_1\\ \widehat\gamma_1=\frac{-\widehat b_2}{1-\widehat a_2}\quad , \quad \widehat\gamma_2=\frac{\widehat a_1-\widehat b_2}{1-\widehat a_2} \end{cases} \tag{4.53} \end{equation}\]

Solving for the original parameters, we get:

• Unique value for \(\widehat b_1\): \(\widehat b_1=-\widehat \beta_2\).
• Unique value for \(\widehat b_2\): \(\widehat b_2=-\widehat \beta_1\).
• Tow possible values for \(\widehat a_2\): \(\widehat a_2=\frac{\widehat\gamma_1-\widehat\beta_1}{\widehat\gamma_1}\) and \(\widehat a_2=\frac{1}{\widehat\alpha_1-\widehat\beta_1}\).
• Four possible values for \(\widehat a_1\): \(\widehat a_1=\frac{\widehat\gamma_2\widehat\beta_1}{\widehat\gamma_1}\), \(\widehat a_1=\frac{-\widehat\beta_2\widehat\gamma_1+\widehat\beta_2\widehat\beta_1-\widehat\alpha_2\widehat\beta_1}{\widehat\gamma_1}\), \(\widehat a_1=\frac{-\widehat\beta_1^2-\widehat\alpha_1\widehat\beta_1+\widehat\gamma_2}{\widehat\alpha_1-\widehat\beta_1}\), and \(\widehat a_1=\frac{-\widehat\beta_2-\widehat\alpha_2\widehat\alpha_1+\widehat\alpha_2\widehat\beta_1+\widehat\alpha_2}{\widehat\alpha_1-\widehat\beta_1}\)

As each of the first two parameters, \(b_1\) and \(b_2\) of the second equation \(I_t=b_1+b_2r_t+\varepsilon_{2t}\) in the system (4.40) has a unique value, it can be said that that equation has been successfully estimated (because it is just identified). Whereas, each parameter of the first equation \(C_t=a_1+a_2Y_t+\varepsilon_{1t}\) has several possible estimates (because it is over-identified). The third equation \(Y_t=C_t+I_t\) however, has no parameter to be estimated and hence not concerned by the identification process.

Now let us check the identification by applying the rank condition and the order condition discussed above on each equation. So we have: \(\mathrm g=3\), \(k=2\),

\[\begin{equation*} \underset{(3,5)}\Pi=(A,B)= \begin{pmatrix} 1&0&-a_2&0&-a_1\\ 0&1&0&-b_2&-b_1\\ -1&-1&1&0&0 \end{pmatrix} \end{equation*}\],

\[\begin{equation*} \underset{(5,2)}\Phi=\begin{pmatrix}\Gamma\\I_2\end{pmatrix}= \begin{pmatrix} \frac{a_2b_2}{1-a_2}&\frac{a_2b_1-a_1}{1-a_2}\\ -b_2&-b_1\\ \frac{-b_2}{1-a_2}&\frac{a_1-b_2}{1-a_2}\\ 1&0\\ 0&1 \end{pmatrix} \end{equation*}\]

And recall our SEM (4.41) for ease of illustration:

\[\begin{equation*} \begin{cases} C_t+0I_t-a_2Y_t+0r_t-a_1u_t=\varepsilon_{1t}\\ 0C_t+I_t+0Y_t-b_2r_t-b_1u_t=\varepsilon_{2t}\\ -C_t-I_t+Y_t+0r_t+0u_t=0 \end{cases} \end{equation*}\]

as shown above, each equation has its own restriction matrix \(R\). So for the first equation, we have two exclusion restrictions (the coefficients of \(I_t\) and \(r_t\) are both equal to zero), so \(R_t\) has two columns (\(q=2\)), and \(\pi_1\) is the first line of the matrix \(\Pi\) as follows:

\[\begin{equation*} \pi_1= \begin{pmatrix} 1&0&-a_2&0&-a_1 \end{pmatrix}\quad , \quad R_t= \begin{pmatrix} 0&0\\1&0\\0&0\\0&1\\0&0 \end{pmatrix} \end{equation*}\]

Therefore, the matrix \((\Phi,R)\) will be:

\[\begin{equation*} (\Phi,R_2)= \begin{pmatrix} \frac{a_2b_2}{1-a_2}&\frac{a_2b_1-a_1}{1-a_2}&0&0\\ -b_2&-b_1&1&0\\ \frac{-b_2}{1-a_2}&\frac{a_1-b_2}{1-a_2}&0&0\\ 1&0&0&1\\ 0&1&0&0 \end{pmatrix} \end{equation*}\]

As we see, the rank of this matrix is \(Rank(\Phi,R)=4\) which satisfy the rule (4.50) since \(\mathrm g+k-1=3+2-1=4\). By the same way, the matrix \((\Phi,R)\) for the second equation will be:

\[\begin{equation*} (\Phi,R_2)= \begin{pmatrix} \frac{a_2b_2}{1-a_2}&\frac{a_2b_1-a_1}{1-a_2}&1&0\\ -b_2&-b_1&0&0\\ \frac{-b_2}{1-a_2}&\frac{a_1-b_2}{1-a_2}&0&1\\ 1&0&0&0\\ 0&1&0&0 \end{pmatrix} \end{equation*}\]

Here the rank condition is also satisfied. Now as both equations are identified, then the model is said to be identified and it can be estimated by one of the methods that will be discussed further.

There exist an easier alternative to check the rank condition as follows(koutsoyiannis.A 1977): The equation of interest is identified if and only if we can find at least one determinant of order \(\mathrm g -1\) composed of the coefficients of the variables excluded from this equation is different from zero. To show how to get those determinants, we follow the below steps using our example again.

1. we take out the whole matrix of the parameters of the system from (4.41):

\[\begin{equation*} \begin{bmatrix} 1&0&-a_2&0&-a_1\\ 0&1&0&-b_2&-b_1\\ -1&-1&1&0&0& \end{bmatrix} \end{equation*}\]

2. To check the rank condition of a particular equation, say the first one, we remove all the columns that have a non zero values in the first row (associated to the first equation):

\[\begin{equation*} \begin{bmatrix} 0&0&\\ 1&-b_2&\\ -1&0& \end{bmatrix} \end{equation*}\]

3. We remove the row of that equation:

\[\begin{equation*} \begin{bmatrix} 1&-b_2&\\ -1&0& \end{bmatrix} \end{equation*}\]

This determinant is:

\[\begin{equation*} \begin{vmatrix} 1&-b_2&\\ -1&0& \end{vmatrix}=-b_2\not=0 \end{equation*}\]

Which means that the first equation is identified (provided that the true \(b_2\not=0\), in other words, the associated variable should be a relevant one). Doing the same with the other equations.

For the second equation then, we will get:

\[\begin{equation*} \begin{vmatrix} 1&-a_2\\ -1&1 \end{vmatrix}=1-a_2\not=0 \end{equation*}\]

The second equation is also identified (provided that the unknown parameter \(a_2\not=1\))

Lastly, the third equation:

\[\begin{equation*} \begin{vmatrix} 0&-a_1\\ -b_2&-b_1 \end{vmatrix}=-b_2a_1\not=0 \end{equation*}\]

The third one is also identified (provided that neither of the parameters \(b_2\) or \(a_1\) is equal to zero), hence the whole system is now identified.

For the order condition, the first equation has two excluded variables, one endogenous (\(I_t\)) and the other one exogenous (\(r_t\)), that is (\((\mathrm g^{'}=2)\geqslant (\mathrm g-1=2)\)). The same holds for the second equation that has two excluded variables \((C_t,Y_t)\). As noted above, the order condition is a weak alternative for the rank condition. It is used only when the latter is not available.

4.2.3.2 Estimation methods

If the simultaneous equations model has some equations under-identified, then those equations can never be estimated as discussed before. However, if all the equations are just-identified, then the estimates of the reduced form by means of OLS will give a unique estimate for each structural parameter. This estimation method is called Indirect list square (ILS), and its estimators are consistent provided that all the classical assumptions are satisfied for the reduced form. Unfortunately, it is rare when all the SEM equations are just-identified. The most models are over-identified, which requires another estimation methods.

If we rely on the unbiasedness criterion instead of the consistency one (in the case of a finite sample with fixed regressors), the structural parameters evaluated by ILS method may not be unbiased even the reduced estimators are unbiased. This is because these two types of parameters may be related by non linear function, as in most cases, and hence the unbiasedness of the original parameters may not be able to be satisfied. For instance, if \(h(x)\) is non linear then \(E\big(h(x)\big)\not=h\big(E(x)\big)\). Using the above example \(E(\widehat\alpha_1)\not=\frac{E(\widehat a_2)E(\widehat b_2)}{1-E(\widehat a_2)}\). In contrast, the consistency will be maintained for the original parameters since it is based on the properties of the probability limit, for instance, \(plim(\widehat\alpha_1)=\frac{plim(\widehat a_2)plim(\widehat b_2)}{1-plim(\widehat a_2)}\).

Two-stage least squares (2SLS):

Since the ILS method, in most cases, gives either no estimates or more different estimates for the same parameter, we can use the 2SLS discussed above in the IV method as an alternative because this method can eliminate the bias resulted from the non-orthogonality condition.

The most critical inconvenience of the IV method resides in the difficulty of finding the instruments. Fortunately, the SEM provides the required instruments, which are the same exogenous variables in the system. They satisfy the orthogonality condition and are correlated with the other endogenous variables through the interactions among the system equations.

The 2SLS can be applied to each equation of the system separately, so as discussed before, the procedure of this method requires the two following steps:

• Step 1: transform the SEM to the reduced form then applying the OLS method on each equation separately. In other words, regressing each variable suspected being endogenous on all the exogenous variables of the system.
• Step 2: regress the endogenous variable of the SEM equation of interest on its exogenous variables along with the fitted values (estimated from step 1) of the endogenous included in this equation.

For instance, suppose that we want to estimate the second equation of (4.36)

\[\begin{equation*} a_{21}y_{1t}+a_{22}y_{2t}+...+a_{2\mathrm g}y_{\mathrm gt}+b_{21}x_{1t}+b_{22}x_{2t}+...+b_{2k}x_{kt}=\varepsilon_{2t} \end{equation*}\]

Or with the appropriate shape:

\[\begin{equation} y_{2t}=\alpha_{21}y_{1t}+\alpha_{23}y_{3t}+...+\alpha_{2\mathrm g}y_{\mathrm gt}+\beta_{21}x_{1t}+\beta_{22}x_{2t}+...+\beta_{2k}x_{kt}+v_{2t} \tag{4.54} \end{equation}\]

The first step of estimation is to estimate all the reduced form parameters by applying the OLS method on each equation separately and retrieve the fitted values as follows:

\[\begin{equation*} \begin{cases} \widehat y_{1t}=\widehat\gamma_{11}x_{1t}+\widehat\gamma_{12}x_{2t}+...+\widehat\gamma_{1k}x_{kt}\\ \widehat y_{2t}=\widehat\gamma_{21}x_{1t}+\widehat\gamma_{22}x_{2t}+...+\widehat\gamma_{2k}x_{kt}\\ ......................................................\\ \widehat y_{\mathrm g t}=\widehat\gamma_{\mathrm g1}x_{1t}+\widehat\gamma_{\mathrm g2}x_{2t}+...+\widehat\gamma_{\mathrm gk}x_{kt} \end{cases} \end{equation*}\]

In matrix form, this system can be rewritten:

\[\begin{equation*} \widehat y_t=X_t\overbrace{(X^t_tX_t)^{-1}X^t_tY_t}^{\widehat\gamma} \end{equation*}\]

Then in the second step, we replace the original endogenous variables of (4.54) by its corresponding fitted values from the last expression to finally get the estimates of the original parameters as follows:

\[\begin{equation} y_{2t}=\widehat\alpha_{21}\widehat y_{1t}+\widehat\alpha_{23}\widehat y_{3t}+...+\widehat\alpha_{2\mathrm g}\widehat y_{\mathrm gt}+\widehat\beta_{21}x_{1t}+\widehat\beta_{22}x_{2t}+...+\widehat\beta_{2k}x_{kt}+v_{2t} \end{equation}\]

Doing the same with the other equations, we will obtain consistent estimators.

Limited information maximum likelihood:

It is similar to the previous method in the sense of the single-equation estimation, but in addition, it assumes the normality distribution of the errors. The asymptotic distribution of its estimators was first derived by Anderson and Rubin (1950)(Anderson 2005).

Consider again our SEM:

\[\begin{equation*} AY_t+BX_t=\varepsilon_t \end{equation*}\]

Where \(E(\varepsilon_t\varepsilon^t_t)=\Sigma\).

The reduced form thus will be:

\[\begin{equation*} Y_t=A^{-1}BX_t+A^{-1}\varepsilon_t=\Pi X_t+v_t \end{equation*}\]

Notice that \(E(v_tv^t_t)=A^{-1}\Sigma(A^{-1})^t=\Omega\).

Now let us suppose that we want to estimate the first equation:

\[\begin{equation*} a_{11}y_{1t}+a_{12}y_{2t}+...+a_{1\mathrm g}y_{\mathrm g t}+b_{11}x_{1t}+b_{12}x_{2t}+...+b_{1k}x_{kt}=\varepsilon_{1t} \end{equation*}\]

Dividing by \(a_{11}\) for simplification purposes as follows:

\[\begin{equation*} y_{1t}+\alpha_{12}y_{2t}+...+\alpha_{1\mathrm g}y_{\mathrm g t}+\beta_{11}x_{1t}+\beta_{12}x_{2t}+...+\beta_{1k}x_{kt}=u_{1t} \end{equation*}\]

Then more compactly in matrix form:

\[\begin{equation} \alpha^tY_t+\beta^tX_t=u_t \tag{4.55} \end{equation}\]

Note that the first element of the vector \(\alpha\) now is equal to 1, and if we apply the same transformations to all the remaining equations, the diagonal elements of the matrix \(A\) are all equal to 1.

Now we should highlight the included variables in our equation (4.55) from the excluded ones as follows:

\[\begin{equation*} \begin{pmatrix} \alpha_{I}^t&0 \end{pmatrix} \begin{pmatrix} Y_{It}\\Y_{Et} \end{pmatrix}+ \begin{pmatrix} \beta_{i}^t&0 \end{pmatrix} \begin{pmatrix} X_{it}\\X_{et} \end{pmatrix}= \begin{pmatrix} u_{It}\\u_{Et} \end{pmatrix} \end{equation*}\]

Which can be reduced to:

\[\begin{equation} \alpha^tY_{It}+\beta^tX_{it}=u_{It} \tag{4.56} \end{equation}\]

Where \(I\) and \(i\) stand for the included endogenous and the included exogenous respectively. Similarly \(E\) and \(e\) corresponding to those excluded from the associated equation. Following these notations, \(\mathrm g_I\) is the number of the included endogenous and \(k_i\) those for the exogenous ones.

One might ask why we use these notations that differ between the excluded and the included variables since the equation of interest does not contain at all those indexed by \(E\) and \(e\). The reason, as we will show shortly, is that this estimation requires the results of the reduced form which contain theses excluded variables.

Now for the reduced form, we can partition the matrix of the regression coefficients using the above notations we get: \[\begin{equation*} \Pi= \begin{pmatrix} \Pi_{Ii}&\Pi_{Ie}\\ \Pi_{Ei}&\Pi_{Ee} \end{pmatrix} \end{equation*}\]

and hence the reduced form can be written as:

\[\begin{equation*} \begin{pmatrix} Y_{It}\\ Y_{Et} \end{pmatrix}= \begin{pmatrix} \Pi_{Ii}&\Pi_{Ie}\\ \Pi_{Ei}&\Pi_{Ee} \end{pmatrix} \begin{pmatrix} X_{it}\\ X_{et} \end{pmatrix}+ \begin{pmatrix} v_{It}\\ v_{Et} \end{pmatrix} \tag{4.57} \end{equation*}\]

and we can define the covariance matrix of \(v_{It}\) as \(E(v_{It}v_{it}^t)=\Omega_{II}\).

According to what was discussed in the identification problem, the matrix \(\Pi\) can be readily estimated by OLS method, but the problem appears when we want to derive the parameters of the structural equation from the relation \(\alpha^t\Pi=\beta^t\).

For our structural equation of interest, this relation can be expressed in terms of the excluded and the included variables as follows:

\[\begin{equation*} \begin{pmatrix} \alpha_I^t&0 \end{pmatrix} \begin{pmatrix} \Pi_{Ii}&\Pi_{Ie}\\ \Pi_{Ei}&\Pi_{Ee} \end{pmatrix}= \begin{pmatrix} \alpha_I^t\Pi_{Ii}&\alpha_I^t\Pi_{Ie} \end{pmatrix}=\begin{pmatrix} \beta_i^t&0 \end{pmatrix} \end{equation*}\]

That is:

\[\begin{equation} \alpha_I^t\Pi_{Ii}=\beta_i^t \tag{4.58} \end{equation}\] \[\begin{equation} \alpha_I^t\Pi_{Ie}=0 \tag{4.59} \end{equation}\]

Where the first expression (4.58) groups \(k_i\) equations and the second one (4.59) \(k_e=k-k_i\) equations.

Since OLS gives consistent estimators for \(\Pi_{Ii}\) and \(\Pi_{Ie}\), in the first step, \(\alpha_I^t\) will be computed from (4.59), then it will be substituted in the second step in (4.58) to compute \(\beta_i^t\). This is exactly the principle of the 2SLS.

However, since we have \(\mathrm g_I\) unknown parameters (to estimate) in the system (4.59) and \(\Pi_{Ie}\) is \(\mathrm g_I \times k_e\) matrix, then to get a possible solution we must have \(k_e\geqslant \mathrm g_I-1\) in which our structural equation is identified.

This means that the solution of (4.59) requires the rank of the matrix \(\Pi_{Ie}\) to be equal to \(\mathrm g_I\) or \(\mathrm g_1-1\). This system looks like a homogeneous one, but in fact, it is not since the first element of \(\alpha_I\) is one. To solve thus this system, it should be rewritten in terms of the estimates \(\widehat\Pi_{Ie}\) as follows:

\[\begin{equation} \alpha_I^t\widehat\Pi_{Ie}=(1,\alpha_I^{*t}) \begin{pmatrix} \widehat\pi\\\widehat\Pi^{*t}_{Ie} \end{pmatrix}\implies\alpha_I^{*t}\widehat\Pi^{*}_{Ie}=-\widehat\pi^t \tag{4.60} \end{equation}\]

Where \(\alpha_I^{*}\) stands for the coefficients of the remaining endogenous variables in our structural equation, \(\widehat\Pi^{*}_{Ie}\) is \((\mathrm g_1-1)\times k_1\) matrix, and \(\widehat\pi\) is the first row of \(\widehat\Pi_{Ie}\).

Now to get a solution, the following condition must be satisfied:

\[\begin{equation} Rank(\widehat\Pi_{Ie})=\mathrm g_I-1 \tag{4.61} \end{equation}\]

Therefore, the maximization of the likelihood function will be subjected to that condition.

However, as that maximization requires too complex computation, we use another alternative called Least variance ratio.

The idea behind this technique is very simple. Assuming that \(\alpha_I\) is known, then we can perform two OLS regressions. The first one is the regression of the transformed response \(Y^*_{It}=\alpha^t_IY_{It}\) on the included regressors \(X_{it}\), and the second one is the regression of the same response \(Y^*_{It}\) on the whole set of the regressors of the system \(X\) (included and excluded). Then compute the sum squared errors of each regression, denoted \(SSR\) for the first, and \(SSR^*\) for the second. Consequently, the sum squared errors never decreases when we add new regressors, and decreases significantly if those regressors are relevant. That is \(SSR^* \leqslant SSR\). But since our structural equation is supposed to have all the relevant variables, the two sums should not be significantly different \(SSR^*\approx SSR\).

Let us rearrange (4.56) as follows:

\[\begin{equation} \alpha^tY_{It}=-\beta^tX_{it}+u_{It} \tag{4.62} \end{equation}\]

The first sum squared errors of this equation thus will be:

\[\begin{align*} SSR&=\alpha^t_IY_{It}(\alpha^t_IY_{It})^t-\alpha^t_IY_{It}X_{it}^t(X_{it}X_{it}^t)^{-1}X_{it}(\alpha^t_IY_{It})^t\\ &=\alpha^t_IY_{It}Y_{It}^t\alpha_I-\alpha^t_IY_{It}X_{it}^t(X_{it}X_{it}^t)^{-1}X_{it}Y_{It}^t\alpha_I\\ &=\alpha^t_I\Big(Y_{It}Y_{It}^t-Y_{It}X_{it}^t(X_{it}X_{it}^t)^{-1}X_{it}Y_{It}^t\Big)\alpha_I\\ &=\alpha^t_IW_I\alpha_I \end{align*}\]

Similarly, the second one will be:

\[\begin{align*} SSR^*&=\alpha^t_I\Big(Y_{It}Y_{It}^t-Y_{It}X^t(XX^t)^{-1}XY_{It}^t\Big)\alpha_I\\ &=\alpha^t_IW\alpha_I \end{align*}\]

now the best estimates of \(\alpha_I\) should minimize the following ratio:

\[\begin{equation*} L=\frac{SSR^*}{SSR}=\frac{\alpha^t_IW_I\alpha_I}{\alpha^t_IW\alpha_I} \end{equation*}\]

Setting the derivatives to zero, we get:

\[\begin{align*} \frac{\partial L}{\partial \alpha_I}&=\frac{\alpha^t_IW\alpha_I\big(2W_I\alpha_I\big)-\alpha^t_IW_I\alpha_I\big(2W\alpha_I\big)}{\big(\alpha^t_IW\alpha_I\big)^2}\\ &=\frac{\frac{\alpha^t_IW\alpha_I}{\alpha^t_IW\alpha_I}\big(2W_I\alpha_I\big)-\frac{\alpha^t_IW_I\alpha_I}{\alpha^t_IW\alpha_I}\big(2W\alpha_I\big)}{\alpha^t_IW\alpha_I}\\ &=\frac{2\big(W_I-LW\big)\alpha_I}{\alpha^t_IW\alpha_I}\\ &=0 \end{align*}\]

Thus, the estimates of \(\alpha_I\) is the solution of:

\[\begin{equation} \big(W_I-LW\big)\alpha_I=0 \tag{4.63} \end{equation}\]

As we see, this system admits a non trivial solution if the determinant is equal to zero:

\[\begin{equation} |W_I-LW|=0 \tag{4.64} \end{equation}\]

This last expression is a polynomial in \(L\), we should thus take its smallest root \(\lambda_{min}\), which is the minimum value of \(L\) (since the criterion is the minimization of \(L\)) and substitute it into (4.63), we get:

\[\begin{equation} \big(W_I-\lambda_{min}W\big)\alpha_I=0 \tag{4.65} \end{equation}\]

Note that the first element of \(\alpha_I\) is equal to one, which ensures getting a unique solution. Finally, substituting these estimates into (4.62), we get the estimates of \(\beta_i\).

Three stage least squares (3SLS):

This method was proposed for the first time by(Zellner and Theil 1962). Unlike the 2SLS that estimates each equation separately, This method estimates the whole system as one block to take into account a possible presence of contemporaneous correlations between errors across the system equations. That is why, the former method is classified within the class of the limited information methods since it does not take all the information contained in the other system equations. Whereas the latter one belongs to the full information methods.

Since this method is based on Generalized least squares GLS and Feasible generalized least squares FGLS (that will be discussed in more detail in the following chapters), we will briefly describe it without going into detail.

Let us rewrite the structural form of the system in (4.36) as follows:

\[\begin{equation*} \begin{cases} y_{1t}=-\frac{a_{12}}{a_{11}}y_{2t}-...-\frac{a_{1\mathrm g}}{a_{11}}y_{\mathrm g t}-\frac{b_{11}}{a_{11}}x_{1t}-\frac{b_{12}}{a_{11}}x_{2t}-...-\frac{b_{1k}}{a_{11}}x_{kt}+\frac{\varepsilon_{1t}}{a_{11}}\\ y_{2t}=-\frac{a_{21}}{a_{22}}y_{1t}-...-\frac{a_{2\mathrm g}}{a_{22}}y_{\mathrm g t}-\frac{b_{21}}{a_{22}}x_{1t}-\frac{b_{22}}{a_{22}}x_{2t}-...-\frac{b_{2k}}{a_{22}}x_{kt}+\frac{\varepsilon_{2t}}{a_{22}}\\ ...........................................................\\ y_{\mathrm g t}=-\frac{a_{\mathrm g1}}{a_{\mathrm g\mathrm g}}y_{1t}-...-\frac{a_{\mathrm g(\mathrm g-1)}}{a_{\mathrm g\mathrm g}}y_{\mathrm g t}-\frac{b_{\mathrm g1}}{a_{\mathrm g\mathrm g}}x_{1t}-\frac{b_{\mathrm g2}}{a_{\mathrm g\mathrm g}}x_{2t}-...-\frac{b_{\mathrm gk}}{a_{\mathrm g\mathrm g}}x_{kt}+\frac{\varepsilon_{\mathrm gt}}{a_{\mathrm g\mathrm g}} \end{cases} \end{equation*}\]

As we see each endogenous variable is expressed in terms of the remaining endogenous and the full set of the original regressors. For instance, the \(i^{th}\) equation will be written as follows:

\[\begin{equation} y_{it}=\alpha_{i1}y_{1t}+...+\alpha_{i\mathrm g}y_{\mathrm gt}+\beta_{i2}x_{2t}+...+\beta_{ik}x_{kt}+u_{it} \tag{4.66} \end{equation}\]

Note that this equation does not have the \(y_{it}\) in the right hand side, which means that the number of the regressors now is \(\mathrm g+k-1\).

If this equation is stacked using the subscript \(t\), the above equation can be rewritten:

\[\begin{equation} \underset{(n\times 1)}{y_{i}}=\underset{(n\times (\mathrm g+k-1))}{Z_i}\underset{((\mathrm g+k-1)\times 1)}{\gamma_i}+\underset{(n\times 1)}{u_i} \tag{4.67} \end{equation}\]

Where \(Z_i=(Y_i,X_i)=(y_{1t},...,y_{\mathrm g t},x_{1t},...,x_{kt})\) (without \(y_{it}\)), and \(\gamma_i=\begin{pmatrix}\alpha_i\\\beta_i\end{pmatrix}\)

Similarly, we can stack all the equations in more general in matrix form as follows:

\[\begin{equation} \underset{(n\mathrm g\times 1)} {\begin{pmatrix} y_1\\..\\y_{\mathrm g} \end{pmatrix}}= \underset{(n\mathrm g\times \mathrm g(\mathrm g+k-1))} {\begin{pmatrix} Z_1&0&0\\0&..&0\\0&0&Z_{\mathrm g} \end{pmatrix}} \underset{(\mathrm g(\mathrm g+k-1)\times 1)} {\begin{pmatrix} \gamma_1\\..\\\gamma_{\mathrm g} \end{pmatrix}}+ \underset{(n\mathrm g\times 1)} {\begin{pmatrix} u_1\\..\\u_{\mathrm g} \end{pmatrix}} \tag{4.68} \end{equation}\]

Or more compactly:

\[\begin{equation} y=Z\gamma+u \tag{4.69} \end{equation}\]

If we estimate this equation by the OLS method we will obtain inconsistent estimators, because the endogenous variables in \(y_i\)’s in the regressors set are correlated with the disturbances \(u_i\)’s in which the orthogonality is not satisfied.

Furthermore, even if the OLS estimators are consistent and the errors of each equation satisfy the classical assumptions (no autocorrelation and homoskedasticity) \(E(u_iu_i^t)=\sigma_i^2I_n\), the possible significant correlations across equations \(E(u_{it}u_{jt})=\sigma_{ij}\) may undermine the efficiency of this type of estimators. Because the variance matrix of the whole system expressed by (Using the Kronecker products):

\[\begin{equation} E(uu^t)=\begin{pmatrix} \sigma_1^2I_n&\sigma_{12}I_n&..&\sigma_{1\mathrm g}I_n\\ \sigma_{21}^2I_n&\sigma_{2}I_n&..&\sigma_{2\mathrm g}I_n\\ ..&..&..&..\\ \sigma_{\mathrm g1}^2I_n&\sigma_{\mathrm g2}I_n&..&\sigma_{\mathrm g}I_n \end{pmatrix}=\Sigma_u \otimes I_n \tag{4.70} \end{equation}\]

Could be not diagonal.

However, the efficiency estimators that should be applied under these conditions is the \(GLS\) ones that will be discussed in the following chapter. If we consider the classical model \(y=X\beta+\varepsilon\), the GLS estimators then will be defined by:

\[\begin{equation*} \begin{cases} b_{GLS}=\big(X^t\Omega_{\varepsilon}^{-1}X\big)^{-1}X^t\Omega_{\varepsilon}^{-1}y\\ Var(b_{GLS})=\big(X^t\Omega_{\varepsilon}^{-1}X\big)^{-1} \end{cases} \end{equation*}\]

Where the non diagonal \(\Omega_{\varepsilon}\) is the variance matrix of the errors when they do not satisfy the assumptions and homoskedasticity.

As required for the OLS, the orthogonality condition is also required for the GLS estimators, which is unfortunately not satisfied since our model contains endogenous regressors in the matrix \(Z\). However, these endogenous regressors \(y_i\) can be replaced by their fitted values \(\widehat y_i\) using the 2SLS estimators such that the GLS estimators will be:

\[\begin{align} \widehat\gamma_{GLS}&=\bigg(\widehat Z^t\big(\Sigma_u \otimes I_n\big)^{-1}\widehat Z\bigg)^{-1}\widehat Z^t\big(\Sigma_u \otimes I_n\big)^{-1}y\\ &=\bigg(\widehat Z^t\big(\Sigma_u^{-1} \otimes I_n\big)\widehat Z\bigg)^{-1}\widehat Z^t\big(\Sigma_u^{-1} \otimes I_n\big)y \tag{4.71} \end{align}\]

Since \(\widehat Z=H_XZ=X(X^tX)^{-1}X^tZ\), we can than obtain the same result by pre-multiplying our model by the \((gk\times gn)\) matrix \((I_g\otimes X^t)\) as follows:

\[\begin{equation} (I_\mathrm g\otimes X^t)y=(I_\mathrm g\otimes X^t)Z\gamma+(I_\mathrm g\otimes X^t)u \tag{4.72} \end{equation}\]

Then the GLS estimators for this transformed model will be:

\[\begin{align*} \gamma_{GLS}&=\bigg[Z^t\big(I_g\otimes X\big)\big(\Sigma_u\otimes X^tX\big)^{-1}\big(I_g\otimes X^t\big)Z\bigg]^{-1}Z^t\big(I_g\otimes X\big)\big(\Sigma_u\otimes X^tX\big)^{-1}\big(I_g\otimes X^t\big)y\\ &=\bigg[Z^t\bigg(\Sigma^{-1}\otimes X\big(X^tX\big)^{-1}X^t\bigg)Z\bigg]^{-1}Z^t\bigg(\Sigma^{-1}\otimes X\big(X^tX\big)^{-1}X^t\bigg)y\\ &=\bigg[Z^t\bigg(\Sigma^{-1}\otimes H_X\bigg)Z\bigg]^{-1}Z^t\bigg(\Sigma^{-1}\otimes H_X\bigg)y \end{align*}\]

Where \(\Sigma_u\otimes X^tX=(I_\mathrm g\otimes X^t)(\Sigma_u\otimes I_n)(I_\mathrm g\otimes X)\) is the variance matrix of the transformed errors \((I_g\otimes X)u\).

To detail the last formula, first, for simplification, let us denote the elements of \(\Sigma_u^{-1}\) by \(\sigma^{ij}\), then by using the properties \(H_X=H^t_X=H^2_X\), we get:

\[\begin{align*} Z^t\big(\Sigma^{-1}_u\otimes H_X\big)Z&= \begin{pmatrix} Z^t_1&0&..&0\\0&Z^t_2&..&0\\..&..&..&0\\0&0&..&Z^t_{\mathrm g} \end{pmatrix} \begin{pmatrix} \sigma^{11}H_X&\sigma^{12}H_X&..&\sigma^{1\mathrm g}H_X\\\sigma^{21}H_X&\sigma^{22}H_X&..&\sigma^{2\mathrm g}H_X\\..&..&..&..\\\sigma^{\mathrm g1}H_X&\sigma^{\mathrm g2}H_X&..&\sigma^{\mathrm g\mathrm g}H_X \end{pmatrix} \begin{pmatrix} Z_1&0&..&0\\0&Z_2&..&0\\..&..&..&0\\0&0&..&Z_{\mathrm g} \end{pmatrix}\\ &=\begin{pmatrix} \sigma^{11}Z^t_1H_X&\sigma^{12}Z^t_1H_X&..&\sigma^{1\mathrm g}Z^t_1H_X\\\sigma^{21}Z^t_2H_X&\sigma^{22}Z^t_2H_X&..&\sigma^{2\mathrm g}Z^t_2H_X\\..&..&..&..\\\sigma^{\mathrm g1}Z^t_{\mathrm g}H_X&\sigma^{\mathrm g2}Z^t_{\mathrm g}H_X&..&\sigma^{\mathrm g\mathrm g}Z^t_{\mathrm g}H_X \end{pmatrix} \begin{pmatrix} Z_1&0&..&0\\0&Z_2&..&0\\..&..&..&0\\0&0&..&Z_{\mathrm g} \end{pmatrix}\\ &=\begin{pmatrix} \sigma^{11}Z^t_1H_XZ_1&\sigma^{12}Z^t_1H_XZ_2&..&\sigma^{1\mathrm g}Z^t_1H_XZ_{\mathrm g}\\\sigma^{21}Z^t_2H_XZ_1&\sigma^{22}Z^t_2H_XZ_2&..&\sigma^{2\mathrm g}Z^t_2H_XZ_{\mathrm g}\\..&..&..&..\\\sigma^{\mathrm g1}Z^t_{\mathrm g}H_XZ_1&\sigma^{\mathrm g2}Z^t_{\mathrm g}H_XZ_2&..&\sigma^{\mathrm g\mathrm g}Z^t_{\mathrm g}H_XZ_{\mathrm g} \end{pmatrix} \end{align*}\]

Using the equality \(Z^t_iH_XZ_j=\widehat Z^t_i\widehat Z^t_j\), the above matrix can be simplified as follows:

\[\begin{align*} Z^t\big(\Sigma_u^{-1}\otimes H_X\big)Z&= \begin{pmatrix} \sigma^{11}\widehat Z^t_1\widehat Z_1&\sigma^{12}\widehat Z^t_1\widehat Z_2&..&\sigma^{1\mathrm g}\widehat Z^t_1\widehat Z_\mathrm g\\ \sigma^{21}\widehat Z^t_2\widehat Z_1&\sigma^{22}\widehat Z^t_2\widehat Z_2&..&\sigma^{2\mathrm g}\widehat Z^t_2\widehat Z_\mathrm g\\ \sigma^{\mathrm g1}\widehat Z^t_\mathrm g\widehat Z_1&\sigma^{\mathrm g2}\widehat Z^t_\mathrm g\widehat Z_2&..&\sigma^{\mathrm g\mathrm g}\widehat Z^t_\mathrm g\widehat Z_\mathrm g \end{pmatrix}\\ &=\bigg(\widehat Z^t\big(\Sigma_u^{-1}\otimes I_n\big)\widehat Z\bigg)^{-1} \end{align*}\]

By following the same previous illustration, we will also have \(Z^t\big(\Sigma_u^{-1}\otimes H_X\big)y=\widehat Z^t\big(\Sigma_u^{-1}\otimes I_n\big)y\). This means that this estimator is equal to that of (4.71).

Finally, the asymptotic variance of this estimator will be:

\[\begin{equation} Var(\widehat\gamma_{GLS})=\bigg(\widehat Z^t\big(\Sigma_u^{-1}\otimes I_n\big)\widehat Z\bigg)^{-1} \tag{4.73} \end{equation}\]

The inconvenience of this formula is that the matrix \(\Sigma_u\) is rarely when it is known, that is why we call the feasible least squares estimators FGLS instead, which merely replaces that unknown matrix by its consistent estimator computed from the 2SLS.

In summary, to compute the FGLS we follow these steps:

1. Compute the 2SLS estimates for each equation \(i\) by:

\[\begin{equation*} \gamma^{GLS}_i=\begin{pmatrix}\alpha^{2SLS}_i\\\beta^{2SLS}_i\end{pmatrix}=\begin{pmatrix}\widehat y^t\widehat y_i&\widehat y^tX_i\\X^t_i\widehat y_i&X^t_iX_i\end{pmatrix}\begin{pmatrix}\widehat y^t_iy_i\\X^t_iy_i\end{pmatrix} \end{equation*}\]

Then retrieve the residuals \(e_i=y_i-Z_i\gamma^{2SLS}_i\).

2. compute the estimated covariance by:

\[\begin{align*} \widehat\sigma_{ij}&=\frac{\big(y_i-Z_i\gamma_i^{2SLS}\big)\big(y_j-Z_j\gamma_j^{2SLS}\big)}{\sqrt{(n-\mathrm g_i-k_i+1)\times (n-\mathrm g_j-k_j+1)}}\\ &\approx\underbrace{\frac{\big(y_i-Z_i\gamma_i^{2SLS}\big)\big(y_j-Z_j\gamma_j^{2SLS}\big)}{n}}_{\text{in arge samples}} \end{align*}\]

That will be the elements of \(\widehat\Sigma_u\). consequently, the FGLS estimator will be defined by:

\[\begin{equation} \widehat\gamma_{FGLS}=\bigg(\widehat Z^t\big(\widehat\Sigma_u^{-1}\otimes I_n\big)\widehat Z\bigg)^{-1}\widehat Z^t\big(\widehat\Sigma_u^{-1}\otimes I_n\big)y \tag{4.74} \end{equation}\]

And its asymptotic variance by:

\[\begin{equation} Var\big(\widehat\gamma_{FGLS}\big)=\bigg(\widehat Z^t\big(\widehat\Sigma_u^{-1}\otimes I_n\big)\widehat Z\bigg)^{-1} \tag{4.75} \end{equation}\]

The 2SLS is preferred to 3SLS if:

• No cross correlation is detected in the system which allows the matrix to be reduced to a constant scalar.
• The sample size \(n\) is not large enough compared to: \(\mathrm g-1\) of the fitted variables, \(k\) of the exogenous regressors, and \(\frac{1}{2}\mathrm g(\mathrm g+1)\) of the parameters of the variance matrix \(\Sigma_u\).

full information maximum likelihood FIML

If the errors are assumed to be at least asymptotically normally distributed \(u\overset{Asy}{\rightarrow}\mathrm N(0,\Sigma_u\otimes I_n)\), we can then call the maximum likelihood method, and when applied on the whole system we would have the following function:

\[\begin{align*} L(\gamma,\Sigma)&=\frac{1}{(2\pi)^{\frac{n\mathrm g}{2}}|\Sigma\otimes I_n|^{\frac{n}{2}}}exp\bigg(-\frac{1}{2}u^t(\Sigma\otimes I_n)^{-1}u\bigg)\\ &=\frac{1}{(2\pi)^{\frac{n\mathrm g}{2}}|\Sigma\otimes I_n|^{\frac{n}{2}}}exp\bigg(-\frac{1}{2}(y-Z\gamma)^t(\Sigma\otimes I_n)^{-1}(y-Z\gamma)\bigg)\bigg\lvert \frac{du}{dy}\bigg\rvert \end{align*}\]

Where the last member \(\bigg\lvert \frac{du}{dy}\bigg\rvert\) stands for the Jacobien.

Remember that the error term is that of the transformed (normalized) system as showed in the equation (4.66). However, to find the jacobien we need to express the transformed system at the time point \(t\) in matrix form as in (4.37). Using the relationship between the two errors (again from (4.66)) \(u_{it}=\frac{\varepsilon_{it}}{a_{it}}\), which can be expressed in matrix form \(u_t=diag\begin{pmatrix}\frac{1}{a_{11}}&\frac{1}{a_{22}}&..&\frac{1}{a_{\mathrm g\mathrm g}}\end{pmatrix}\varepsilon_t\), hence \(\varepsilon_t=diag\begin{pmatrix}a_{11}&a_{22}&..&a_{\mathrm g\mathrm g}\end{pmatrix}u_t=Du_t\). Substituting this result into (4.37)and pre-multiplying by \(D^{-1}\), we get:

\[\begin{equation*} D^{-1}AY_t+D^{-1}BX_t=u_t \end{equation*}\]

staking all the equations according the subscripts \(t\), the Jacobian then will be:

\[\begin{equation*} \bigg\lvert\frac{du}{dy}\bigg\rvert=\bigg\lvert D^{-1}A\bigg\rvert^n \end{equation*}\]

Where, the elements of the matrix \(D^{-1}A\) are merely the endogenous regressors coefficients of the original structural form:

\[\begin{equation*} D^{-1}A=\begin{pmatrix} a_{11}&a_{12}&..&a_{1\mathrm g}\\ a_{21}&a_{22}&..&a_{2\mathrm g}\\ ..&..&..&..\\ a_{\mathrm g1}&a_{\mathrm g2}&..&a_{\mathrm g\mathrm g} \end{pmatrix} \end{equation*}\]

Now substituting the Jacobian and moving on to the log-likelihood, we get:

\[\begin{equation*} l\big(\gamma, \Sigma\big)=-\frac{n\mathrm g}{2}ln(2\pi)-\frac{n}{2}ln\bigg\lvert\Sigma\bigg\rvert+nln\bigg\lVert D^{-1}A\bigg\rVert-\frac{1}{2}\big(y-Z\gamma\big)^t\big(\Sigma\otimes I_n\big)^{-1}\big(y-Z\gamma\big) \end{equation*}\]

Where, \(\bigg\lVert D^{-1}A\bigg\rVert\) is the absolute value of the \(D^{-1}A\).

Note that the parameters contained in \(\gamma\) are transformed from those of \(D^{-1}A\).

As we can see this function is non linear in the parameters, which requires some non linear method that will be discussed later on.

Example 4.3 For the sake of simplicity, we will simulate an SEM with relevant and non-relevant exogenous ones. For this simulation, we will use R, then copy the data to python. The population SEM will be the following:

\[\begin{equation*} \begin{cases} y_1=4+1.2y_2+3.5x_1+0.5x_2 \\ y_2=1.5+0.7y_1-1.7y_3+3.5x_1+0.5x_2-0.5x_3\\ y_3=2.5+1.7y_1+2.5x_1+0.9x_3 \end{cases} \end{equation*}\]

This system has three endogenous variables, and three exogenous ones in total.

In R:

To simulate this model we will use the R package lavaan.

# install the library
#install.packages('lavaan', dependencies = TRUE)

# call the library 
library('lavaan')

[out] This is lavaan 0.7-2
[out] lavaan is FREE software! Please report any bugs.

[out] 
[out] Attaching package: 'lavaan'

[out] The following object is masked from 'package:psych':
[out] 
[out]     cor2cov

# specify the true model
true_model <- ' y1 ~ 4 + 1.2*y2 + 3.5*x1
                y2 ~ 1.5 + 3.5*x1 + 0.5*x2 + 2.5*x3
                y3 ~ 2.5 + 1.7*y1 + 2.5*x1 + 0.9*x3
              '

# generate 200samples from that model
set.seed(111)
data_sem_R <- simulateData(true_model, sample.nobs=100L)
head(data_sem_R)

[out]           y1         y2        y3         x1         x2         x3
[out] 1 -1.8251544 -0.3464216 -4.247240 -0.2416494  0.4114033  0.4909933
[out] 2  3.2776508 -0.2213862  5.789617  0.4629265  0.3642892 -0.8104085
[out] 3  2.4688658  1.8364487  5.227607  0.2622169 -0.6626252  0.7814895
[out] 4 18.9686847  9.7234294 39.296428  2.0474443  2.4942754  1.0559264
[out] 5  1.7220537 -0.3892297  3.045059  0.2283535  0.8169536 -0.4507162
[out] 6 -0.9845232 -1.6411584 -2.217717  0.2613143 -0.7530267 -0.8084898

Now let us fit a SEM model with the correct specification..

# specify the model
sem_model <- ' 
              y1 ~ 1 + y2 + x1
              y2 ~ 1 + x1 + x2 + x3
              y3 ~ 1 + y1 + x1 + x3
             '
# firt the model
model_fitted <- sem(sem_model, data=data_sem_R)

We can plot the diagramm of that model using another package called semPlot.

suppressPackageStartupMessages(library("semPlot"))
semPaths(model_fitted, title = FALSE, curvePivot = TRUE)

Figure 4.1: the model path

# display the results
summary(model_fitted)

[out] lavaan 0.7-2 ended normally after 1 iteration
[out] 
[out]   Estimator                                         ML
[out]   Optimization method                           NLMINB
[out]   Number of model parameters                        14
[out] 
[out]   Number of observations                           100
[out] 
[out] Model Test User Model:
[out]                                                       
[out]   Test statistic                                 2.081
[out]   Degrees of freedom                                 4
[out]   P-value (Chi-square)                           0.721
[out] 
[out] Parameter Estimates:
[out] 
[out]   Standard errors                             Standard
[out]   Information                                 Expected
[out]   Information saturated (h1) model          Structured
[out] 
[out] Regressions:
[out]                    Estimate  Std.Err  z-value  P(>|z|)
[out]   y1 ~                                                
[out]     y2                1.194    0.033   35.988    0.000
[out]     x1                3.354    0.151   22.257    0.000
[out]   y2 ~                                                
[out]     x1                3.606    0.085   42.313    0.000
[out]     x2                0.352    0.099    3.563    0.000
[out]     x3                2.474    0.084   29.392    0.000
[out]   y3 ~                                                
[out]     y1                1.705    0.067   25.511    0.000
[out]     x1                2.574    0.520    4.955    0.000
[out]     x3                0.928    0.219    4.240    0.000
[out] 
[out] Intercepts:
[out]                    Estimate  Std.Err  z-value  P(>|z|)
[out]    .y1                0.141    0.094    1.500    0.134
[out]    .y2               -0.117    0.091   -1.278    0.201
[out]    .y3               -0.093    0.099   -0.943    0.346
[out] 
[out] Variances:
[out]                    Estimate  Std.Err  z-value  P(>|z|)
[out]    .y1                0.882    0.125    7.071    0.000
[out]    .y2                0.811    0.115    7.071    0.000
[out]    .y3                0.976    0.138    7.071    0.000

The sem function uses by default the maximum likelihood estimation ML. However, it has a bunch of other methods like GLS in case of autoccorrelation and/or heteroskedasticity, WLS weighted least squares, etc.

The test statistic used is the Chi-square, we see that it has value 22.884 (small p-value). It is mainly used to compare different SEM’s.

This model is over-identified with 4 degrees for freedom allowing the addition of more variables in the SEM.

We see that all the estimates are significant, except for the intercepts.

We have also the covariance matrix of the errors across the SEM equations.

Suppose now that we do not know the correct specification of our model, so that we add the irrelevant variable \(x_2\) to the first equation, and remove the relevant variable \(x_3\) from the second equation. This time will display some tests in the output

# specify the model
sem_model2 <- ' 
              y1 ~ 1 + y2 + x1 + x2 
              y2 ~ 1 + x1 + x2
              y3 ~ 1 + y1 + x1 + x3
             '
# fit the model
model_fitted2 <- sem(sem_model2, data=data_sem_R)
summary(model_fitted2, fit.measures=TRUE)

[out] lavaan 0.7-2 ended normally after 2 iterations
[out] 
[out]   Estimator                                         ML
[out]   Optimization method                           NLMINB
[out]   Number of model parameters                        14
[out] 
[out]   Number of observations                           100
[out] 
[out] Model Test User Model:
[out]                                                       
[out]   Test statistic                               228.663
[out]   Degrees of freedom                                 4
[out]   P-value (Chi-square)                           0.000
[out] 
[out] Model Test Baseline Model:
[out] 
[out]   Test statistic                              1368.336
[out]   Degrees of freedom                                12
[out]   P-value                                        0.000
[out] 
[out] User Model versus Baseline Model:
[out] 
[out]   Comparative Fit Index (CFI)                    0.834
[out]   Tucker-Lewis Index (TLI)                       0.503
[out] 
[out] Loglikelihood and Information Criteria:
[out] 
[out]   Loglikelihood user model (H0)               -521.016
[out]   Loglikelihood unrestricted model (H1)       -406.685
[out]                                                       
[out]   Akaike (AIC)                                1070.032
[out]   Bayesian (BIC)                              1106.504
[out]   Sample-size adjusted Bayesian (SABIC)       1062.288
[out] 
[out] Root Mean Square Error of Approximation:
[out] 
[out]   RMSEA                                          0.749
[out]   90 Percent confidence interval - lower         0.669
[out]   90 Percent confidence interval - upper         0.834
[out]   P-value H_0: RMSEA <= 0.050                    0.000
[out]   P-value H_0: RMSEA >= 0.080                    1.000
[out] 
[out] Standardized Root Mean Square Residual:
[out] 
[out]   SRMR                                           0.139
[out] 
[out] Goodness of Fit Index:
[out] 
[out]   Goodness of Fit Index (GFI)                    1.000
[out]   90 Percent confidence interval - lower         1.000
[out]   90 Percent confidence interval - upper         1.000
[out] 
[out] Parameter Estimates:
[out] 
[out]   Standard errors                             Standard
[out]   Information                                 Expected
[out]   Information saturated (h1) model          Structured
[out] 
[out] Regressions:
[out]                    Estimate  Std.Err  z-value  P(>|z|)
[out]   y1 ~                                                
[out]     y2                1.194    0.034   35.535    0.000
[out]     x1                3.354    0.152   22.014    0.000
[out]     x2                0.001    0.104    0.007    0.995
[out]   y2 ~                                                
[out]     x1                3.684    0.264   13.929    0.000
[out]     x2                0.490    0.306    1.600    0.110
[out]   y3 ~                                                
[out]     y1                1.705    0.028   60.551    0.000
[out]     x1                2.574    0.237   10.863    0.000
[out]     x3                0.928    0.092   10.063    0.000
[out] 
[out] Intercepts:
[out]                    Estimate  Std.Err  z-value  P(>|z|)
[out]    .y1                0.141    0.095    1.486    0.137
[out]    .y2                0.058    0.282    0.205    0.837
[out]    .y3               -0.093    0.099   -0.943    0.346
[out] 
[out] Variances:
[out]                    Estimate  Std.Err  z-value  P(>|z|)
[out]    .y1                0.882    0.125    7.071    0.000
[out]    .y2                7.816    1.105    7.071    0.000
[out]    .y3                0.976    0.138    7.071    0.000

Before interpreting this output we should introduce the main three fit indexes (in addition to the popular one, AIC, BIC, and LogLik) RMSA, CFI, and TLI, that are used in practice instead of Chi-square that is less informative in large samples. The first one is an absolute measure of fit unlike the last ones are used to compare between models. For instance, they assess whether adding or removing some variables could improve the model.

Confirmatory Factor Index CFI:

In the model_fitted2, we see a test named \(CFI\) with the value \(0.844\). By default, the sem function compared the current model named User Model with baseline model that regresses every dependent variable on only the intercept. This index calculates the ratio of the deviation of the current model from the baseline model against the deviation of the saturated model (the just identified model) from the latter one. It is calculated as follows:

\[\begin{equation*} CFI=\frac{\delta(baseline)-\delta(current)}{\delta(baseline)} \end{equation*}\]

Where \(\delta=\chi^2-df\). The \(CFI\) will be then:

\[\begin{equation*} CFI=\frac{(1355.354-12)-(213.770-4)}{(1355.354-12)}=0.844 \end{equation*}\]

Values closer to 1 indicate a good fit.

Tucker Lewis Index TLI:

It works like the previous one with a slight difference. It is defined by:

\[\begin{align*} CFI&=\frac{\chi(baseline)\big/df(baseline)-\chi(current)\big/df(current)}{\chi(baseline)\big/df(baseline)-1} \\ &= \frac{1355.354\big/12-213.770\big/4}{1355.354\big/12-1} \\ &= 0.523 \end{align*}\]

values closer to 1 indicate a good fit.

The root mean square error of approximation RMSA:

This index uses the \(\delta=\chi^2-df\) as a measure of misspecification. The greater the \(\delta\) the more accurately specified the model. It is defined by:

\[\begin{align*} RMSA&=\sqrt{\frac{\delta}{df(N-1)}} \\ &=\sqrt{\frac{213.770-4}{4(100-1)}} \\ &=0.728 \end{align*}\]

The cutoff criteria are:

• \(RMSA\leqslant 0.05\) good fit.
• \(0.05\leqslant RMSA\leqslant 0.08\) reasonable approximation, neither good fit nor poor fit.
• \(RMSA\geqslant 0.1\) poor fit.

By inspecting those indexes, RMSA indicates poor fits. CFI and TLi are not high enough. This is inline with the fact that we have removed a relevant variable from the model. Let us check those values from the first correct model.

summary(model_fitted, fit.measures=TRUE)

[out] lavaan 0.7-2 ended normally after 1 iteration
[out] 
[out]   Estimator                                         ML
[out]   Optimization method                           NLMINB
[out]   Number of model parameters                        14
[out] 
[out]   Number of observations                           100
[out] 
[out] Model Test User Model:
[out]                                                       
[out]   Test statistic                                 2.081
[out]   Degrees of freedom                                 4
[out]   P-value (Chi-square)                           0.721
[out] 
[out] Model Test Baseline Model:
[out] 
[out]   Test statistic                              1368.336
[out]   Degrees of freedom                                12
[out]   P-value                                        0.000
[out] 
[out] User Model versus Baseline Model:
[out] 
[out]   Comparative Fit Index (CFI)                    1.000
[out]   Tucker-Lewis Index (TLI)                       1.004
[out] 
[out] Loglikelihood and Information Criteria:
[out] 
[out]   Loglikelihood user model (H0)               -407.725
[out]   Loglikelihood unrestricted model (H1)       -406.685
[out]                                                       
[out]   Akaike (AIC)                                 843.450
[out]   Bayesian (BIC)                               879.923
[out]   Sample-size adjusted Bayesian (SABIC)        835.707
[out] 
[out] Root Mean Square Error of Approximation:
[out] 
[out]   RMSEA                                          0.000
[out]   90 Percent confidence interval - lower         0.000
[out]   90 Percent confidence interval - upper         0.110
[out]   P-value H_0: RMSEA <= 0.050                    0.802
[out]   P-value H_0: RMSEA >= 0.080                    0.114
[out] 
[out] Standardized Root Mean Square Residual:
[out] 
[out]   SRMR                                           0.001
[out] 
[out] Goodness of Fit Index:
[out] 
[out]   Goodness of Fit Index (GFI)                    1.000
[out]   90 Percent confidence interval - lower         1.000
[out]   90 Percent confidence interval - upper         1.000
[out] 
[out] Parameter Estimates:
[out] 
[out]   Standard errors                             Standard
[out]   Information                                 Expected
[out]   Information saturated (h1) model          Structured
[out] 
[out] Regressions:
[out]                    Estimate  Std.Err  z-value  P(>|z|)
[out]   y1 ~                                                
[out]     y2                1.194    0.033   35.988    0.000
[out]     x1                3.354    0.151   22.257    0.000
[out]   y2 ~                                                
[out]     x1                3.606    0.085   42.313    0.000
[out]     x2                0.352    0.099    3.563    0.000
[out]     x3                2.474    0.084   29.392    0.000
[out]   y3 ~                                                
[out]     y1                1.705    0.067   25.511    0.000
[out]     x1                2.574    0.520    4.955    0.000
[out]     x3                0.928    0.219    4.240    0.000
[out] 
[out] Intercepts:
[out]                    Estimate  Std.Err  z-value  P(>|z|)
[out]    .y1                0.141    0.094    1.500    0.134
[out]    .y2               -0.117    0.091   -1.278    0.201
[out]    .y3               -0.093    0.099   -0.943    0.346
[out] 
[out] Variances:
[out]                    Estimate  Std.Err  z-value  P(>|z|)
[out]    .y1                0.882    0.125    7.071    0.000
[out]    .y2                0.811    0.115    7.071    0.000
[out]    .y3                0.976    0.138    7.071    0.000

We see the values now are perfect.

Lastly, the interesting function that helps us to decide what variable should include is the following.

# display the first 5 rows
modindices(model_fitted2, sort=TRUE, maximum.number = 5)

[out]    lhs op rhs     mi    epc sepc.lv sepc.all sepc.nox
[out] 31  y2  ~  x3 89.625  2.474   2.474    0.552    0.515
[out] 30  y2  ~  y3 86.179  2.499   2.499    9.258    9.258
[out] 45  x3  ~  y2  1.570  0.006   0.006    0.028    0.030
[out] 33  y3  ~  x2  0.686 -0.091  -0.091   -0.005   -0.005
[out] 26  y2 ~~  y3  0.625  0.709   0.709    0.257    0.257

The most important value is the modification index (mi) that gives the Chi-square change if we add the associated path. So the function correctly finds the right path that we should include back the relevant variable \(x3\) into the second equation.

In Python:

Fortunately, Python has a package called semopy with approximately the same syntax as lavaan.

if not 'data_sem_py' in globals():
  data_sem_py = r.data_sem_R

First, We call the library, then fit the misspecified model discussed above. Note that the constant term will not be created automatically, like in lavaan, so it should be included explicitly.

import semopy

# specify the model
sem_model2_py = ''' 
              y1 ~   y2 + x1 + x2 
              y2 ~   x1 + x2
              y3 ~   y1 + x1 + x3
             '''
# To fit the models with intercepts, we use ModelMeans method
model_fitted2_py = semopy.ModelMeans(sem_model2_py)
model_res = model_fitted2_py.fit(data=data_sem_py)

# display the results
model_res

[out] SolverResult(fun=481.18803457264676, success=True, n_it=96, x=array([ 3.35951313,  0.02038922,  0.14339896,  3.63536879,  0.90764208,
[out]        -0.15175066,  2.29279811,  0.78308009, -0.10807555,  1.19270083,
[out]         1.74166029,  0.88267892,  7.23447175,  0.95870977]), message='Optimization terminated successfully', name_method='SLSQP', name_obj='FIML')

In addition to the estimation methods in R package lavaan, this package has the Full Information Maximum Likelihood FIML (which is not implemented yet in lavaan). As we see the default method is MLM, but we can change it to, say GLS, through the argument obj in the fit function. Note that the optimization method is from scipy package (the default SLSQP).

To show the estimates, we should type the following.

model_ins=model_fitted2_py.inspect()
model_ins

[out]    lval  op rval  Estimate  Std. Err    z-value       p-value
[out] 0    y1   ~   y2  1.192701  0.034930  34.145506  0.000000e+00
[out] 1    y1   ~   x1  3.359513  0.154938  21.682995  0.000000e+00
[out] 2    y1   ~   x2  0.020389  0.105833   0.192655  8.472291e-01
[out] 3    y2   ~   x1  3.635369  0.254151  14.303997  0.000000e+00
[out] 4    y2   ~   x2  0.907642  0.289071   3.139856  1.690311e-03
[out] 5    y3   ~   x1  2.292798  0.235655   9.729451  0.000000e+00
[out] 6    y3   ~   x3  0.783080  0.092621   8.454650  0.000000e+00
[out] 7    y3   ~   y1  1.741660  0.028014  62.170020  0.000000e+00
[out] 8    y1   ~    1  0.143399  0.095076   1.508255  1.314893e-01
[out] 9    y2   ~    1 -0.151751  0.271767  -0.558384  5.765819e-01
[out] 10   y3   ~    1 -0.108076  0.097965  -1.103203  2.699390e-01
[out] 11   y1  ~~   y1  0.882679  0.124830   7.071068  1.537437e-12
[out] 12   y2  ~~   y2  7.234472  1.023109   7.071068  1.537437e-12
[out] 13   y3  ~~   y3  0.958710  0.135582   7.071068  1.537437e-12

Since all the intercepts are insignificant, we should remove them.

# To fit the models without intercepts, we use Model method
model_fitted3_py = semopy.Model(sem_model2_py)
model_res2 = model_fitted3_py.fit(data=data_sem_py)

To get the indexes discussed above, we type the following.

model_index=semopy.calc_stats(model_fitted3_py)
model_index.T

[out]                      Value
[out] DoF              10.000000
[out] DoF Baseline     18.000000
[out] chi2            220.012161
[out] chi2 p-value      0.000000
[out] chi2 Baseline  1370.634352
[out] CFI               0.844738
[out] GFI               0.839482
[out] AGFI              0.711067
[out] NFI               0.839482
[out] TLI               0.720529
[out] RMSEA             0.460580
[out] AIC              17.599757
[out] BIC              46.256629
[out] LogLik            2.200122

As we see, this package has extras indexes over lavaan package.