Nonparametric Instrumental Regression with Two-Way Fixed Effects

2.1 Nonparametric IV Regression with One-Way FE

Fève and Florens (2014) consider the a panel model given by

where Y i t ∈ R and Z i t ∈ R p are observable random variables, ξ i ∈ R is an unobserved individual effect, and U _it is an error term. Here, Z _i1, …, Z _iT ∀ i is supposed to be endogenous in that it might be correlated with both ξ _i and U _it .^[4] Estimating ξ _i  ∀ i would lead into the well-known incidental parameter problem (Lancaster 2000; Neyman and Scott 1948). Hence, to eliminate the FE they first apply first-difference, i.e. Y _it − Y _i,t−1, which however still does not fully solve the endogeneity issue as E(φ(Z _it ) − φ(Z _i,t−1)|U _it − U _i,t−1) ≠ 0. Supposing that there exists a valid set of instruments W ∈ R q , satisfying E(U _it − U _i,t−1|W) = 0, they then proceed by applying the Tikhonov regularization method (Tikhonov and Arsenin 1977) to solve the resulting problem below for φ, which is given by

where ( K φ ) ( w ) = E ( φ ( Z i t ) − φ ( Z i , t − 1 ) | W = w ) ∈ L W 2 and r ( w ) = E ( Y i t − Y i , t − 1 | W = w ) ∈ L W 2 .^[5] Here, K is a conditional mean operator projecting functions of Z onto the space of W, i.e. K : L Z 2 → L W 2 . Hence, (Kφ)(w) can be expressed by

which is a Fredholm equation of the first kind (Fredholm 1903).

Consider again equation (2), by iterating the projection onto the space of the endogenous variable Z, one obtains

where K* is an adjoint conditional expectation operator projecting functions of W onto the space of function of Z, i.e. K * : L W 2 → L Z 2 .

As described by the authors, simply inverting K*K for obtaining an estimate of φ results in an ill-posed inverse problem as the solution of (2) for φ requires the inversion of an integral which does not represent a continuous mapping. That is, while small changes in Kφ might induce small changes in r, this does not hold vice versa.^[6] Therefore, regularization methods are required, where, in this case, the Tikhonov regularization seeks to minimize

w.r.t. φ, where α > 0 is a regularization parameter to be chosen by the user, and with

Fève and Florens (2014) show that the solution to (5) is then given by

The intuition behind Tikhonov regularization is to add the constant α to the eigenvalues of K*K which in turn allows for inversion. Finally, the function φ is estimated by replacing the conditional mean objects K, K*, and r by appropriate nonparametric estimates.

As outlined, the here presented method relies on “global first-difference” transformation to eliminate the FE. However, the expansion of the approach to properly deal with a two-way FE model seems to be rather complicated. Hence, by the use of the below described nonparametric estimator for the conditional mean function by taking into account FE, I will develop a nonparametric instrumental estimator that is appealing by its simple implementation and higher flexibility to treat both one-way and two-way FE panel models.

2.2 Nonparametric Regression with One- and Two-Way FE

A method to take potential endogeneity from FE into account was recently presented by LMU. They consider a similar model with regard to (1), however, suppose that only ξ _i (but not U _it ) can be correlated with Z _i1, …, Z _iT . To consistently estimate the unknown function φ, they use what they call the local-within kernel estimator, which eliminates the FE locally. The motivation for the local-within kernel estimator lies in the fact, which they prove, that global-demeaning or global first-difference to remove the unobserved individual effect – i.e. by transformation the data as it is practiced in many studies as well as in the estimation procedure in Fève and Florens (2014) – introduces a non-degenerating bias, even for a large time horizon. Instead, removing the individual effect locally, i.e. by local-within or local-first-difference, avoids this bias.

3.1 Setup of the Estimator

Consider the panel model given by

with the observable random variables Y i t ∈ R and Z i t ∈ R (for simplicity with only one covariate). The function of interest φ is supposed to belong to the space of square integrable functions of Z, denoted by L z 2 . The endogenous explanatory variable Z _it is potentially correlated with the unobserved individual and time specific effects, ξ i ∈ R and δ t ∈ R , as well as with the error term U _it . If there exists an instrument W ∈ R satisfying E(U _it |W = w) = E(U _it |w) = 0, we can write

where, given this characterization, φ(Z _it ) is the solution of a Fredholm integral equation of the first kind.^[11] Given the joint density f(y, z, ξ, δ, w) and using the common abuse of notation f for denoting different densities, (20) can be expressed as

where the two middle integrals cannot be computed since ξ _i and δ _t are unobserved. Further, let K φ ̃ = E ( φ ( Z i t ) + ξ i + δ t | w ) , where K denotes again the conditional expectation operator, projecting functions of Z onto the space of functions of W, i.e. K : L z 2 → L w 2 (Centorrino, Feve, and Florens 2017) and φ ̃ ≡ φ ( Z i t ) + ξ i + δ t , and let r = E(Y _it |w). Then, equation (20) can be expressed by

Similar to what was shown in the context of the estimator presented by Fève and Florens (2014) discussed above, the challenge now is to solve (22) for φ ̃ – which is an ill-posed inverse problem – while recovering φ. To be more precise, the ill-posed inverse problem here again arises since the inversion of the operator K, to solve for φ, is discontinuous, which leads to inconsistent estimates when not regularized. To deal with the ill-posed inverse problem I make use of the Landweber–Fridman regularization (Fridman 1965; Landweber 1951), where I follow Racine (2019, p. 282). As will be shown, the motivation to use the Landweber–Fridman instead of the Tikhonov regularization is that it allows to combine nonparametric instrumental regression with the advantageous features of the local-within LMU estimator to control for the unobserved FE (i.e. to flexibly deal with one- or two-way FE).

Let K* be the adjoint conditional mean operator of K projecting functions of W onto the space of functions of Z _t , i.e. K * : L w 2 → L z 2 . In other words, the adjoint conditional mean operator K* is the reverse projection of K (for more details see for instance Florens, Johannes, and Van Bellegem (2011)).^[12] Consider equation (22), taking the scalar product with respect to K* and multiplying with a constant c yields

where, in contrast to the Tikhonov regularization used in Fève and Florens (2014), the Landweber–Fridman regularization does not aim to invert K*K but avoids its inversion using instead the iterative scheme obtained by equivalently writing

Using the Landweber–Fridman regularization, setting c < 1, φ ̃ can be obtained in an iterative manner, by

and so

In practice, the regularized estimator of φ is obtained by substituting the conditional mean objects r, K, and K* in equation (23) by consistent estimators. To consistently estimate these objects, the unobserved FE, ξ _i and δ _t , need to be controlled for – whose appearance is explicitly shown in the equivalent expression (24) – which, as shown below, will be effectively done by the use of the LMU two-way FE estimator defined in equation (17).^[13] The regularization procedure is then conducted for k = 1 , … , k ̄ iterations until convergence, with k ̄ the total number of iterations, determined by a stopping rule described in the following sub-section.

3.2 Implementation of the Estimation Algorithm

As mentioned above, to compute φ _k (z) various conditional mean objects need to be consistently estimated, i.e. controlling for FE, where the regularization procedure involves to be repeated (i.e. updated) over several iterations, where the last iteration, k ̄ , should yield φ ̂ k ̄ ( z ) ≈ φ ( z ) . To achieve this objective, the following estimation algorithm is employed:

Step 1. Compute the initial guess φ ̂ 0 ( z ) .

To start with the procedure, an initial guess is estimated by considering E(Y _it |z, ξ _i , δ _t ) = φ ₀(z) + ξ _i + δ _t , supposing E(ξ _i |z) ≠ 0 and E(δ _t |z) ≠ 0. That is, if we regressed Y _it on Z _it using conventional kernel regression methods we would obtain biased estimates for φ ₀(z). Instead, to account for FE the LMU estimator is applied, which consists in locally demeaning the dependent and explanatory variable by applying equation (15), obtaining Y i t * * ( z ) and Z i t * * ( z ) , and using local-linear kernel regression of Y i t * * ( z ) on Z i t * * ( z ) . In doing so, according to the LMU approach, we first obtain the (pointwise) gradient, i.e. φ ̂ 0 ′ ( z ) , whereupon the corresponding conditional mean function, φ ̂ 0 ( z ) , is recovered by equation (17).

Step 2. Compute E ̂ Y ̃ i t | w , ξ i , δ t with Y ̃ i t ≡ Y i t − φ ̂ 0 ( Z i t ) .

The previous step yielded and unbiased estimate of φ ̂ 0 ( Z i t ) where FE were accounted for. When subtracting φ ̂ 0 ( Z i t ) from Y _it it can be shown that the empirical model to be considered in this step consists in E Y ̃ i t | w , ξ i , δ t = φ Y ̃ | w 0 ( w ) + ξ i + δ t , where φ Y ̃ | w 0 ( w ) is the conditional mean function we aim to consistently estimate.^[14] That is, if the instrument is correlated with the FE, the LMU estimator has to be applied by first obtaining the locally demeaned variables and conduct local-linear kernel regression of Y ̃ i t * * ( w ) on W**(w), which yields an estimate of the gradient, φ ̃ Y ̃ | w ′ 0 ( w ) , and recover subsequently φ ̂ Y ̃ | w 0 ( w ) ≡ E ̂ ( Y ̃ i t | w , ξ i , δ t ) by (17).^[15]

Step 3. Compute E ̂ E ̂ Y ̃ i t | w , ξ i , δ t | z , ξ i , δ t .

Let Y ̃ ̂ i t ≡ E ̂ Y ̃ i t | w , ξ i , δ t , obtained from the previous step. The regression performed in this step builds on the empirical model given by E Y ̃ ̂ | z , ξ i , δ t = φ Y ̃ ̂ | z 0 ( z ) + ξ i + δ t , where φ Y ̃ ̂ | z 0 ( z ) is the conditional mean function we aim to consistently estimate. Since, again, E(ξ _i |z) ≠ 0 and E(δ _t |z) ≠ 0 is assumed, the LMU estimator is applied by computing the respective locally demeaned variables and use local-linear regression of Y ̃ ̂ i t * * ( z ) on Z i t * * ( z ) yielding first an estimate of the gradient φ Y ̃ ̂ | z ′ 0 ( z ) and recover subsequently φ ̂ Y ̃ ̂ | z 0 ( z ) ≡ E ̂ E ̂ Y ̃ i t | w , ξ i , δ t | z , ξ i , δ t by (17).

Step 4. Compute φ ̂ 1 ( z ) = φ ̂ 0 ( z ) + c E ̂ E ̂ Y ̃ i t | w , ξ i , δ t | z , ξ i , δ t .

For k = 1 we now have an estimate φ ̂ 1 ( z ) , while having effectively controlled for unobserved two-way FE.

Step 5. Repeat steps 2–4, replacing φ ̂ 0 ( z ) by φ ̂ 1 ( z ) to compute φ ̂ 2 ( z ) and so on.

Step 6. This is continued until the stopping criterion given by

stabilizes throughout the iterations.^[16]

To establish the stopping criterion E(Y _it |w, ξ _i , δ _t ) = φ _Y|w (Y|w) + ξ _i + δ _t is estimated by the LMU approach before the algorithm starts. That is, we aim to consistently estimate φ _Y|w (w) under the assumption that potentially E(ξ _i |w) ≠ 0 and E(δ _t |w) ≠ 0. Hence, the LMU estimator is employed by locally regressing in a first step the locally demeaned variables Y i t * * ( w ) on W**(w) (yielding an estimate of the gradient φ Y | w ′ ( w ) ) and, subsequently, φ ̂ Y | w ( w ) ≡ E ̂ ( Y i t | w , ξ i , δ t ) is recovered. The second object in the stopping criterion, E ̂ ( φ ̂ k ( Z i t ) | w ) , needs to be computed and updated in each iteration. Here, the LMU estimator is likewise employed, allowing to control for potential unobserved FE. More precisely, let φ ̂ k , i t ≡ φ ̂ k ( Z i t ) , and regress the locally demeaned variables φ ̂ k , i t * * ( w ) on W i t * * ( w ) and recover the corresponding conditional mean function subsequently.

4.1 Data Generating Process (DGP)

I extend the DGP presented in Darolles et al. (2011) to the case of endogeneity characterized by correlation between the explanatory variable and both the error term and the unobserved two-way FE. More precisely, I generate a panel dataset setting n = 100 and T = 20, yielding nT = 2000 observations. The unobserved two-way FE are generated by^[18]

Next, an auxiliary explanatory variable z ̃ i t is generated as a function of the FE, given by

From these draws, the instrumental variable w _it is generated by

where ρ _zw = 0.2 denotes the correlation coefficient between the explanatory and the instrumental variable, and v w , i t ∼ N ( 0,0.15 ) denotes and error term. Only now the generation of the endogenous explanatory variable is completed by

where v z , i t ∼ N ( 0,0.8 ) . Finally, the dependent variable y _it is obtained by

where φ ( z i t ) = z i t 2 refers to as the true DGP and u _it = ρ _uz v _z,it + ϵ _it , with ρ _uz = 0.8 and ϵ i t ∼ N ( 0,0.05 ) .

Table 1 presents the covariance matrix to illustrate the introduced endogeneity. As can be seen, the endogenous explanatory variable z _it is correlated with all other important components, i.e. with the instrument, w _it , with the unobserved individual and temporal effects, ξ _i and δ _t , as well as with the error term, u _it . Note that the instrument is also slightly correlated with the individual/temporal effects but uncorrelated with the error term, which is necessary to serve as a valid instrument.

4.2 Estimation Results

Figure 1 illustrates the regularized solution path to obtain the final estimate. The initial guess estimate, i.e. φ ̂ 0 ( z ) , is indicated by the bottom solid blue line. The iterative estimation procedure then approaches the true conditional mean function φ(z), indicated by the dashed red line. The estimation resulting from the last iteration, indicated by the green line, corresponding to φ ̂ k ̄ ( z ) , is then very close to the true curve, i.e. the DGP. For better visibility, Figure 2 only shows the initial guess (blue line) as well as the last estimate of the iteration (green line).

Figure 1:

The regularized solution path starting from the initial guess φ ̂ 0 ( z ) until reaching the final IV regression estimate φ ̂ k ̄ ( z ) .

Figure 2:

Comparison between the initial guess φ ̂ 0 ( z ) and the final IV regression estimate. φ ̂ k ̄ ( z ) .

As mentioned above, the number of iterations is determined by the stopping criterion given in (25), i.e. the iterations continue until there is no significant change in the value of this criterion between two successive iterations. Figure 3 shows the evolution of the value of the stopping criterion throughout the iterations. As illustrated, the criterion’s value decreases from the first iteration on until it flattens out. More precisely, after k ̄ = 44 iterations, the value of the stopping criterion only changes by 1 % or less, which is where the iteration procedure stops.

Figure 3:

Evolution of the function value of the stopping criterion.

Finally, to illustrate possible biases I compare different nonparametric estimators, based on the described DGP. An overview of the compared estimators is presented in Table 2.

In particular, I compare (i) simple local-linear estimation (LL), where no source of endogeneity is taken into account at all; (ii) the local-within FE estimator presented by Lee et al. (2019) (LMU), that is, only taking into account individual and temporal FE; (iii) the Landweber–Fridman procedure (L-F), i.e. only applying nonparametric IV regression without taking into account FE; and (iv) the here presented estimator, the Landweber–Fridman regularization in combination with the local-within estimator (L-F/LMU), i.e. nonparametric IV and FE regression;

Figure 4 shows the results applying the different estimators. It can be seen that all estimators, except the L-F/LMU estimator, lead to biased estimates, as none of the corresponding estimated conditional mean functions hits the true DGP, indicated by the dashed red line. As expected, the L-F/LMU estimator is closely located around the true curve, which suggests a good finite sample behaviour.

Figure 4:

Comparison between nonparametric estimators.

4.3 Bootstrap Inference and Finite Sample Behaviour

To allow for statistical inference of the presented estimator, I apply the wild residual block bootstrap (described in more detail in Appendix B) to estimate confidence intervals (CI). Here, the α/2 × 100 % and (1 − α/2) × 100 % pointwise percentiles are estimated from the empirical distribution of B = 400 bootstrap estimates φ ̂ 1 ( Z i t ) , … , φ ̂ B ( Z i t ) . For α = 0.05 this yields the pointwise 95 % confidence interval of φ ̂ ( Z i t ) . Figure 5 illustrates that the used bootstrap method provides very narrow confidence intervals suggesting statistically reliable estimates of φ ̂ ( Z i t ) , except at the boundaries where the data is scarce, which is typical for nonparametric kernel estimation. While the employed bootstrap method yields convincing results, providing proof of its validity in this context is not trivial and beyond the scope of this paper.

Figure 5:

Pointwise estimates of the regularized conditional mean function along with the bootstrapped 95 % confidence interval (CI), using 400 replications.

Further, to properly investigate the proposed estimator’s finite sample performance, I conduct a Monte Carlo simulation using different sample sizes. More specifically, building on the Monte Carlo simulation shown by Lee, Mukherjee, and Ullah (2019), I use four sets of samples of n ∈ {50, 100} and T ∈ {10, 20}, each generated from the above described DGP. Then, φ(Z _it ) from (19) is estimated based on 400 random samples (for each of the four sets) using the estimators L-F/LMU, L-F, LMU, and LL (see Table 2 for a description of the estimators). Further, for each of the estimators, the integrated mean squared errors (IMSE) and the integrated mean absolute errors (IMAE) are computed by averaging the pointwise MSE and MAE over z, where the RIMSE reports the root of IMSE. Typically, Monte Carlo simulations with nonparametric estimators are computationally time intensive, not only because of pointwise estimation, but also because, ideally, optimal bandwidths should be used for each repetition. Beyond that, the nonparametric instrumental estimators, based on the Landweber–Fridman regularization method, need various iterations themselves to yield the final estimate. The available computational power, however, does not allow me to conduct the Monte Carlo simulation using updated cross-validated optimal bandwidths for each of the repetitions. Instead, for each of the compared estimators and each of the four samples, I compute optimal bandwidths only for the first of the 400 repetitions and which will then be used for the remaining repetitions. While this might introduce some inaccuracy at the boundaries of the data, it does not change the quality of the results. This is because optimal bandwidths are only marginally affected by resampling from the same distributions or more generally from the same DGP. The employed Monte Carlo simulation does therefore allow to conjecture on the estimator’s finite sample behaviour.^[19]

Table 3 presents the results. For all datasets, the proposed L-F/LMU estimator, taking endogeneity in the explanatory variable and FE into account, reaches the lowest values both for the RIMSE and the IMAE, indicating the best estimation performance among the applied nonparametric estimators. Considering the RIMSE, it can be seen that increasing the sample size improves the here proposed L-F/LMU estimator’s performance. Especially when increasing the time dimension from T = 10 to T = 20 leads to a significant improvement. The IMAE shows a similar pattern, where the best performance of the L-F/LMU estimator is shown for the largest dataset with n = 100 and T = 20, suggesting convergence and a good finite sample behaviour of the estimator. Generally, the RIMSE is higher compared to the IMSE because the RIMSE squares the errors before averaging, which puts higher weights on large errors that occur particularly at the boundaries. This effect becomes somewhat amplified when not using up-dated (optimal) bandwidths for each of the 400 random samples. Hence, the IMAE might here be preferably considered to evaluate the estimators’ performance.

Among the other estimators, the LMU estimator shows the best performance, as it grasps a large part of the unobserved heterogeneity. Instead, the L-F estimator and the local-linear estimator do not show any significant improvement when increasing the sample size. It would be interesting to conduct the same simulation based on different DGP’s to confirm and further study the performance of the estimator.