--- title: "Interactive Fixed Effects for Balanced and Unbalanced Panels" author: "Binzhi Chen" date: "`r Sys.Date()`" output: rmarkdown::html_vignette: toc: true toc_depth: 3 number_sections: true vignette: > %\VignetteIndexEntry{Interactive Fixed Effects for Balanced and Unbalanced Panels} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include=FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 4 ) ``` # Introduction ## Why interactive fixed effects? Standard two-way fixed effects (TWFE) controls for unobserved heterogeneity through an additive structure: $$y_{it} = \alpha_i + \xi_t + X_{it}'\beta + u_{it}$$ The unit effect $\alpha_i$ absorbs any time-invariant unit characteristic; the time effect $\xi_t$ absorbs any unit-invariant time shock. This works well when heterogeneous responses to common shocks are negligible. Many empirical settings violate this assumption. Consider: - **Globalisation shocks**: countries are exposed differently depending on their trade openness (a unit characteristic that interacts with a common world-demand shock). - **Monetary policy**: firms with different financial leverage respond asymmetrically to a common interest-rate change. - **COVID-19**: regional outcomes depend on a national shock weighted by pre-existing health-system capacity. In each case, the unobserved confounder has the product form $\lambda_i' F_t$: a unit-specific loading $\lambda_i$ multiplied by a common factor $F_t$. If $X_{it}$ is correlated with $\lambda_i' F_t$, TWFE is inconsistent even as $N, T \to \infty$. ## The IFE model Bai (2009) proposes replacing the additive error structure with: $$y_{it} = \alpha_i + \xi_t + X_{it}'\beta + \lambda_i' F_t + u_{it} \qquad (1)$$ where: - $\beta \in \mathbb{R}^p$: structural parameters of interest, - $F_t \in \mathbb{R}^r$: common factors ($T \times r$ matrix $\hat{F}$, normalised $\hat{F}'\hat{F}/T = I_r$), - $\lambda_i \in \mathbb{R}^r$: unit-specific factor loadings ($N \times r$ matrix $\hat{\Lambda}$), - $u_{it}$: idiosyncratic error with $\mathbb{E}[u_{it} \mid X, F, \Lambda] = 0$. The $r$ interactive terms $\lambda_{ik} F_{kt}$ generalise both additive unit effects ($r=1$, $F_t = 1$) and additive time effects ($r=1$, $\lambda_i = 1$). With $r$ unrestricted, (1) subsumes TWFE as the special case $r = 0$. ## What `xtife` provides `xtife` is a pure base-R package (no compiled code, no external dependencies) offering: | Feature | Function | Reference | |---|---|---| | Balanced IFE: analytical SEs | `ife()` | Bai (2009) | | Balanced IFE: static bias correction | `ife(..., bias_corr=TRUE)` | Bai (2009, §7) | | Balanced IFE: dynamic regressors | `ife(..., method="dynamic")` | Moon & Weidner (2017) | | Factor number selection (balanced) | `ife_select_r()` | Bai & Ng (2002); Bai (2009) | | Unbalanced IFE: exact inference | `ife_unbalanced()` | Su, Wang & Wang (2025) | | Unbalanced IFE: factor selection | `ife_select_r_unb()` | Su, Wang & Wang (2025, §3.3) | --- # Balanced Panel IFE ## Data The package includes the `cigar` dataset: 46 US states $\times$ 30 years (1963--1992) of cigarette sales and prices from Baltagi (1995). ```{r data} library(xtife) data(cigar) dim(cigar) head(cigar[, c("state", "year", "sales", "price")], 3) ``` The panel is balanced: every state appears in every year. ## Quick start ```{r quickstart} fit <- ife(sales ~ price, data = cigar, index = c("state", "year"), r = 2, force = "two-way", se = "cluster") print(fit) ``` The `force` argument controls additive fixed effects: `"none"` | `"unit"` | `"time"` | `"two-way"` (default). Key output fields: ```{r components} fit$coef # named coefficient vector fit$se # standard errors fit$ci # 95% confidence intervals (matrix: coef x 2) fit$converged # TRUE if outer loop converged fit$n_iter # number of outer iterations dim(fit$F_hat) # T x r estimated factors dim(fit$Lambda_hat) # N x r estimated loadings ``` ## Estimation algorithm The Bai (2009) estimator alternates between two steps until convergence: **Step 1 — Factor extraction.** Given $\hat\beta$, form $W_{it} = \ddot y_{it} - \ddot X_{it}'\hat\beta$ (residuals after additive FE demeaning). The top-$r$ eigenvectors of $WW'/NT$ (scaled to satisfy $\hat F' \hat F / T = I_r$) give $\hat F$. The loadings follow by regression: $\hat\Lambda = W'\hat F / T$. **Step 2 — Coefficient update.** Project the demeaned regressors onto the orthogonal complement of $\hat F$: $$\tilde X_{it} = \ddot X_{it} - \hat F_t (\hat F' \hat F)^{-1} \hat F' \ddot X_i$$ Then $\hat\beta = (\tilde X'\tilde X)^{-1}\tilde X'\ddot y$ (Frisch-Waugh-Lovell). The outer loop converges when $\max|\hat\beta^{(k+1)} - \hat\beta^{(k)}| < \text{tol}$ (default $10^{-9}$). ## Standard error types ```{r se_compare} fit_std <- ife(sales ~ price, data = cigar, index = c("state", "year"), r = 2, se = "standard") fit_rob <- ife(sales ~ price, data = cigar, index = c("state", "year"), r = 2, se = "robust") fit_cl <- ife(sales ~ price, data = cigar, index = c("state", "year"), r = 2, se = "cluster") data.frame( se_type = c("standard", "robust (HC1)", "cluster"), coef = round(c(fit_std$coef, fit_rob$coef, fit_cl$coef), 5), se = round(c(fit_std$se, fit_rob$se, fit_cl$se), 5), t_stat = round(c(fit_std$tstat, fit_rob$tstat, fit_cl$tstat), 3) ) ``` All three use the sandwich formula $\widehat{\mathrm{Var}}(\hat\beta) = A^{-1} B A^{-1}$ where $A = \tilde X'\tilde X$ and: | `se=` | $B$ matrix | Degrees of freedom | |---|---|---| | `"standard"` | $\hat\sigma^2 A$ | $NT - p - r(N+T-r) - \text{FE dof}$ | | `"robust"` | $\sum_{it} \hat u_{it}^2 \tilde x_{it}\tilde x_{it}'$ (HC1) | same | | `"cluster"` | $\sum_i \psi_i \psi_i'$, $\psi_i = \sum_t \hat u_{it}\tilde x_{it}$ | same | **Rule of thumb**: use `"cluster"` for most macro panels; `"robust"` when clustering by unit is not appropriate; `"standard"` only as a benchmark. ## Comparison with TWFE Setting $r = 0$ reduces `ife()` to the standard TWFE estimator: ```{r twfe_compare} fit0 <- ife(sales ~ price, data = cigar, index = c("state", "year"), r = 0, force = "two-way") # Manual within-transformation cigar$y_dm <- cigar$sales - ave(cigar$sales, cigar$state) - ave(cigar$sales, cigar$year) + mean(cigar$sales) cigar$x_dm <- cigar$price - ave(cigar$price, cigar$state) - ave(cigar$price, cigar$year) + mean(cigar$price) lm0 <- lm(y_dm ~ x_dm - 1, data = cigar) cat(sprintf("ife (r=0): %.7f\n", fit0$coef["price"])) cat(sprintf("lm TWFE: %.7f\n", coef(lm0)["x_dm"])) cat(sprintf("diff: %.2e\n", abs(fit0$coef["price"] - coef(lm0)["x_dm"]))) ``` The difference is at machine precision, confirming algebraic equivalence. TWFE gives $\hat\beta \approx -0.38$; IFE with $r=2$ gives $\approx -0.52$: the gap reflects the OVB from ignoring interactive confounders. ## Factor number selection `ife_select_r()` fits models for $r = 0, \ldots, r_{\max}$ and reports five information criteria. ```{r select_r, eval=FALSE} # Not run during build (takes ~20 s); set verbose=FALSE in scripts sel <- ife_select_r(sales ~ price, data = cigar, index = c("state", "year"), r_max = 6, force = "two-way", verbose = FALSE) attr(sel, "suggested") # named vector: r chosen by each IC ``` The five criteria: | Criterion | Penalty | Notes | |---|---|---| | IC1 | $\hat\sigma^2 (N+T-r)\log(NT)/(NT)$ | Bai & Ng (2002) | | IC2 | $\hat\sigma^2 (N+T-r)\log(\min(N,T)^2)/(NT)$ | Bai & Ng (2002) | | IC3 | $\hat\sigma^2 (N+T-r)\log(\min(N,T))/(NT)$ | Bai & Ng (2002) | | IC(BIC) | $V(r)\log(NT)/(NT)$ | Bai (2009, §4.3) | | PC | $V(r) + r \cdot \hat\sigma^2 (N+T)\log(NT)/(NT)$ | Bai (2009) | **Recommendation**: IC(BIC) for small-to-moderate panels ($\min(N,T) < 60$); IC1 or IC2 for large panels. ## Asymptotic bias correction (balanced panel) ### Static regressors — Bai (2009) The IFE estimator has asymptotic bias $O(1/N + 1/T)$ from the incidental parameters in $\hat\Lambda$ and $\hat F$. The analytical correction is: $$\hat\beta^\dagger = \hat\beta - \hat B / N - \hat C / T$$ where $\hat B$ and $\hat C$ are closed-form expressions involving $\hat F$, $\hat\Lambda$, and the regression residuals (Bai 2009, Theorem 2). ```{r bias_static} fit_bc <- ife(sales ~ price, data = cigar, index = c("state", "year"), r = 2, se = "standard", bias_corr = TRUE) cat(sprintf("Raw IFE: %.5f\n", fit_bc$coef_raw["price"])) cat(sprintf("Corrected:%.5f\n", fit_bc$coef["price"])) cat(sprintf("B_hat: %.5f (B_hat/N = %.5f)\n", fit_bc$B_hat["price"], fit_bc$B_hat["price"] / fit_bc$N)) cat(sprintf("C_hat: %.5f (C_hat/T = %.5f)\n", fit_bc$C_hat["price"], fit_bc$C_hat["price"] / fit_bc$T)) ``` The correction is most important when $T/N$ is non-negligible. For the cigar panel ($N=46$, $T=30$, $T/N \approx 0.65$) it shifts the coefficient by about 0.007. ### Dynamic regressors — Moon and Weidner (2017) When regressors include lagged outcomes or are otherwise only *weakly* exogenous (correlated with past errors), the static bias correction is insufficient. The `method = "dynamic"` option applies the Moon and Weidner (2017) double-projection algorithm and three-term bias correction: $$\hat\beta^* = \hat\beta + \hat W_X^{-1}\!\left(\frac{\hat B_1}{T} + \frac{\hat B_2}{N} + \frac{\hat B_3}{T}\right)$$ $\hat B_1$ (Nickell-type) is near zero for approximately exogenous regressors. ```{r bias_dynamic} fit_dyn <- ife(sales ~ price, data = cigar, index = c("state", "year"), r = 2, se = "standard", method = "dynamic", bias_corr = TRUE, M1 = 1L) cat(sprintf("Dynamic BC coef: %.5f\n", fit_dyn$coef["price"])) cat(sprintf("B1/T (Nickell): %.5f\n", fit_dyn$B1_hat["price"] / fit_dyn$T)) ``` For `price` (approximately exogenous), $\hat B_1 / T \approx 0$, so the static and dynamic corrected estimates nearly coincide. --- # Unbalanced Panel IFE ## Motivation Real-world panels are rarely balanced. Survey non-response, merger events, data-collection gaps, and rotating samples all cause missing $(i,t)$ cells. Simply deleting units with incomplete records (list-wise deletion) induces selection bias; imputation by column means distorts the factor structure. Su, Wang and Wang (2025, hereafter SWW2025) extend the Bai (2009) estimator to genuinely unbalanced panels. Their paper provides: 1. A two-step estimation procedure: nuclear-norm regularisation (NNR) for consistent initialisation of the factor structure, followed by an Alternating Maximisation (AM) outer loop that iterates between updating $\beta$ and $(\hat F, \hat\Lambda)$; within each AM iteration the factors are updated via the EM algorithm of Bai (2009, Appendix B); 2. Analytical standard errors yielding $\sqrt{NT}$-asymptotically normal inference; 3. Analytical bias correction for both strictly and weakly exogenous regressors. `xtife` implements the complete SWW2025 framework. ## The unbalanced IFE model Let $d_{it} \in \{0,1\}$ indicate whether $(i,t)$ is observed. The model is: $$y_{it} = X_{it}'\beta + \lambda_i' F_t + u_{it} \qquad \text{for all } (i,t) \text{ with } d_{it}=1 \qquad (2)$$ There are no additive fixed effects in (2). Individual and time heterogeneity is absorbed by the factor structure: a unit effect corresponds to a factor with $F_t = 1$ for all $t$; a time effect corresponds to a factor with $\lambda_i = 1$ for all $i$. In practice, increase $r$ by 1 to absorb a unit additive effect, or by 2 to absorb both unit and time additive effects (two-way heterogeneity). ### Identification The normalisation $\hat F' \hat F / T = I_r$ and $\hat\Lambda' \hat\Lambda$ diagonal identify $F$ and $\Lambda$ up to an $r \times r$ rotation. Only the column space of $F$ (the factor space) is identified, which suffices for inference on $\beta$. ## Quick start with unbalanced data ```{r unbalanced_setup} # Create a 10% randomly missing version of the cigar panel set.seed(42) cigar_unb <- cigar[sample(nrow(cigar), size = floor(0.9 * nrow(cigar))), ] cat(sprintf("Observations: %d / %d (%.0f%% fill)\n", nrow(cigar_unb), nrow(cigar), 100 * nrow(cigar_unb) / nrow(cigar))) ``` ```{r unbalanced_fit} fit_unb <- ife_unbalanced(sales ~ price, data = cigar_unb, index = c("state", "year"), r = 2L, se = "cluster") print(fit_unb) ``` Key output fields mirror those of `ife()`: ```{r unbalanced_fields} fit_unb$coef # estimated beta fit_unb$se # standard errors fit_unb$ci # 95% confidence intervals fit_unb$converged # outer-loop convergence fit_unb$n_obs # number of observed (i,t) pairs fit_unb$N # number of units fit_unb$TT # number of distinct time periods dim(fit_unb$F_hat) # TT x r dim(fit_unb$Lambda_hat) # N x r ``` ## Estimation algorithm ### Outer alternating-maximisation (AM) loop `ife_unbalanced()` solves the unbalanced least-squares problem: $$(\hat\beta, \hat F, \hat\Lambda) = \arg\min_{\beta, F, \Lambda} \sum_{i,t} d_{it}\left(y_{it} - X_{it}'\beta - \lambda_i' F_t\right)^2$$ subject to $F'F/T = I_r$. The outer loop alternates: 1. **Factor step**: given $\hat\beta$, form $W_{it} = y_{it} - X_{it}'\hat\beta$ and apply the EM inner loop (below) to recover $(\hat F, \hat\Lambda)$. 2. **Coefficient step**: project out the factor space from each unit's observed data and run pooled OLS on the projected regressors $\tilde X_{it} = X_{it} - \hat F_t(\hat F_i'\hat F_i)^{-1}\hat F_i' X_i$, where the inner product is restricted to unit $i$'s observed periods. ### EM inner loop (Bai 2009, Appendix B) — used within each AM iteration Within each AM outer iteration, $(\hat F, \hat\Lambda)$ are updated by the EM algorithm of Bai (2009, Appendix B), which treats the missing $W_{it}$ as latent variables imputed from the current factor estimate: 1. **E-step**: impute missing cells with the current fitted values: $W_{it}^{\text{fill}} = \hat\Lambda_i'\hat F_t$ for all unobserved $(i,t)$. 2. **M-step**: apply standard SVD factor extraction to the completed $T \times N$ matrix $W^{\text{fill}}$. The EM loop converges when $\max_{(i,t): d_{it}=1} |\hat\Lambda_i^{(k+1)}\!\!'\hat F_t^{(k+1)} - \hat\Lambda_i^{(k)}\!\!'\hat F_t^{(k)}| < \text{tol}_{\text{EM}}$ (convergence on fitted values at observed cells; the rotational indeterminacy of $F$ makes direct comparison of $\hat F$ meaningless). ### Initialisation Two options for the initial $(\hat\beta^{(0)}, \hat F^{(0)}, \hat\Lambda^{(0)})$: **`init = "ols"` (default):** plain OLS on observed cells gives $\hat\beta^{(0)}$; the EM inner loop is then run from $(\hat F, \hat\Lambda) = 0$. **`init = "nnr"` (nuclear-norm regularisation):** solves the nuclear-norm penalised problem $$(\hat\beta^{(0)}, \hat\Theta^{(0)}) = \arg\min_{\beta, \Theta} \tfrac{1}{2}\!\sum_{it} d_{it}(y_{it} - \Theta_{it} - X_{it}'\beta)^2 + \nu_{NT}\|\Theta\|_*$$ via soft-impute iteration. The penalty $\nu_{NT}$ is cross-validated over $c \cdot \sqrt{\max(N,T)}$ with $c \in \{0.01, 0.1, 1, 10\}$. NNR initialisation is recommended when the panel is severely unbalanced (fill rate $< 60\%$) or has a rich factor structure ($r \geq 3$). For mildly unbalanced panels, both initialisations converge to the same fixed point. ```{r init_compare} # Both converge to the same coefficient for mild unbalancedness fit_ols <- ife_unbalanced(sales ~ price, data = cigar_unb, index = c("state", "year"), r = 2, init = "ols") fit_nnr <- ife_unbalanced(sales ~ price, data = cigar_unb, index = c("state", "year"), r = 2, init = "nnr") cat(sprintf("OLS init: coef = %.5f, %d outer iterations\n", fit_ols$coef["price"], fit_ols$n_iter)) cat(sprintf("NNR init: coef = %.5f, %d outer iterations\n", fit_nnr$coef["price"], fit_nnr$n_iter)) ``` ## Exact standard errors (SWW2025 Theorem 4.1) The exact sandwich variance of $\hat\beta$ requires projecting each regressor out of both the factor space and the loading space, giving *doubly-projected* regressors $\hat x_{itk}$: $$\hat x_{itk} = x_{itk} - \hat\delta_{ki}'\hat F_t - \hat\omega_{kt}'\hat\lambda_i$$ where $\hat\delta_{ki}$ and $\hat\omega_{kt}$ solve the alternating least-squares system (SWW2025, eq. 4.3): $$\hat\delta_{ki} = \Bigl(\sum_t d_{it}\hat F_t\hat F_t'\Bigr)^{-1} \!\sum_t d_{it}\hat F_t(x_{itk} - \hat\lambda_i'\hat\omega_{kt})$$ $$\hat\omega_{kt} = \Bigl(\sum_i d_{it}\hat\lambda_i\hat\lambda_i'\Bigr)^{-1} \!\sum_i d_{it}\hat\lambda_i(x_{itk} - \hat F_t'\hat\delta_{ki})$$ The sandwich variance is then $$\widehat{\mathrm{Var}}(\hat\beta) = \hat W_X^{-1} \hat\Omega_X \hat W_X^{-1}$$ with $\hat W_X = \frac{1}{NT}\sum_{it} d_{it}\hat x_{it}\hat x_{it}'$ and $\hat\Omega_X$ depending on the `se` argument: | `se=` | $\hat\Omega_X$ | Assumption | |---|---|---| | `"standard"` | $\frac{1}{NT}\sum_{it} d_{it}\hat u_{it}^2\hat x_{it}\hat x_{it}'$ | Homoskedastic errors | | `"robust"` | same formula, HC1-scaled | Heteroskedastic errors | | `"cluster"` | $\frac{1}{NT}\sum_i \psi_i\psi_i'$, $\psi_i = \sum_t d_{it}\hat u_{it}\hat x_{it}$ | Serial correlation within units | | `"hac"` | Bartlett kernel over $(t,s)$ pairs: $\frac{1}{NT}\sum_i\sum_{t,s}\Gamma\!\left(\frac{|s-t|}{L_T}\right)d_{it}d_{is}\hat u_{it}\hat u_{is}\hat x_{it}\hat x_{is}'$ | Serial correlation over time within units | The Bartlett kernel is $\Gamma(u) = (1-u)\mathbf{1}\{u \leq 1\}$; the bandwidth defaults to $L_T = \lfloor 2T^{1/5}\rfloor$. ```{r se_unb} fit_std <- ife_unbalanced(sales ~ price, data = cigar_unb, index = c("state", "year"), r = 2, se = "standard") fit_rob <- ife_unbalanced(sales ~ price, data = cigar_unb, index = c("state", "year"), r = 2, se = "robust") fit_cl <- ife_unbalanced(sales ~ price, data = cigar_unb, index = c("state", "year"), r = 2, se = "cluster") fit_hac <- ife_unbalanced(sales ~ price, data = cigar_unb, index = c("state", "year"), r = 2, se = "hac") data.frame( se_type = c("standard", "robust", "cluster", "hac"), coef = round(c(fit_std$coef, fit_rob$coef, fit_cl$coef, fit_hac$coef), 5), se = round(c(fit_std$se, fit_rob$se, fit_cl$se, fit_hac$se), 5) ) ``` **When to use `"hac"`**: when the errors $u_{it}$ are serially correlated (e.g., in dynamic models or when residuals show autocorrelation). For strictly exogenous regressors with i.i.d. errors, `"cluster"` is sufficient. ## Factor number selection (SWW2025 §3.3) `ife_select_r_unb()` uses the singular value thresholding (SVT) rule applied to the nuclear-norm regularised matrix $\hat\Theta^{(0)}$: $$\hat r = \sum_{s=1}^{N \wedge T} \mathbf{1}\!\left\{ \sigma_s\!\left(\frac{\hat\Theta^{(0)}}{\sqrt{NT}}\right) \geq c_f\sqrt{c_{NT}^{-1/4}\,\sigma_1\!\left(\frac{\hat\Theta^{(0)}}{\sqrt{NT}}\right)} \right\}$$ where $c_{NT} = \min(\sqrt{N}, \sqrt{T})$ and $c_f = 0.6$ (the default threshold constant, recommended by SWW2025 for balanced panels and shown to work well for mild unbalancedness). ```{r sel_unb, eval=FALSE} # Not run during build (NNR cross-validation takes ~10 s) sel_unb <- ife_select_r_unb(sales ~ price, data = cigar_unb, index = c("state", "year"), verbose = FALSE) cat(sprintf("SVT selects r_hat = %d\n", sel_unb$r_hat)) ``` The function returns (invisibly): - `r_hat`: the selected number of factors (at least 1); - `sv`: the normalised singular values used in the threshold comparison; - `threshold`: the computed threshold value; - `nu_used`: the NNR penalty chosen by cross-validation. ## Analytical bias correction (SWW2025 Theorem 4.2) Even with correct $r$, the unbalanced IFE estimator has asymptotic bias $O((NT)^{-1/2})$ from incidental parameters. Setting `bias_corr = TRUE` applies the SWW2025 analytical correction $\hat\beta^{abc} = \hat\beta - (NT)^{-1/2}\hat W_X^{-1}\hat b$ where $\hat b = \hat b_3 + \hat b_4 + \hat b_5 + \hat b_6$ (strictly exogenous case): $$\hat b_{3k} = \frac{1}{\sqrt{NT}}\sum_{i,t} d_{it}\hat u_{it}^2 \hat\delta_{ki}'[\hat L_{ff}^{-1}]_t\hat\lambda_i$$ $$\hat b_{4k} = \frac{1}{\sqrt{NT}}\sum_{i,t} d_{it}\hat u_{it}^2 \hat\omega_{kt}'[\hat L_{\lambda\lambda}^{-1}]_i\hat F_t$$ $$\hat b_{5k} = \frac{1}{\sqrt{NT}}\sum_{i,t} d_{it}\hat u_{it}^2 \hat\lambda_i'\hat\Xi_{kt}\hat\lambda_i$$ $$\hat b_{6k} = \frac{1}{\sqrt{NT}}\sum_{i,t} d_{it}\hat u_{it}^2 \hat F_t'\hat\Delta_{ki}\hat F_t$$ Here $[\hat L_{ff}^{-1}]_t = (\sum_i d_{it}\hat\lambda_i\hat\lambda_i')^{-1}$, $[\hat L_{\lambda\lambda}^{-1}]_i = (\sum_t d_{it}\hat F_t\hat F_t')^{-1}$, and $\hat\Delta_{ki}$, $\hat\Xi_{kt}$ are $r \times r$ auxiliary matrices defined in SWW2025 (page 21). ### Strictly exogenous regressors (`exog = "strict"`) Use when $\mathbb{E}[u_{it} \mid X, F, \Lambda] = 0$ for all $(i,t)$. This is the default and covers the vast majority of cross-sectional panel applications. ```{r bc_strict} fit_bc_strict <- ife_unbalanced(sales ~ price, data = cigar_unb, index = c("state", "year"), r = 2, se = "standard", bias_corr = TRUE, exog = "strict") cat(sprintf("Raw beta: %.5f\n", fit_bc_strict$coef_raw["price"])) cat(sprintf("Corrected beta: %.5f\n", fit_bc_strict$coef["price"])) cat(sprintf("b_hat (b3+...+b6): %.5f\n", fit_bc_strict$b_hat["price"])) ``` ### Weakly exogenous regressors (`exog = "weak"`) Use when the regressor includes the lagged dependent variable or is otherwise predetermined (correlated with past but not future errors). This adds the dynamic bias term $\hat b_2$: $$\hat b_{2k} = \frac{1}{\sqrt{NT}}\sum_i\sum_{t < s} \Gamma\!\left(\frac{s-t}{L_T}\right) \hat u_{it}\, x_{isk}\,\hat F_t'[\hat L_{\lambda\lambda}^{-1}]_i\hat F_s$$ The total bias correction becomes $\hat b = \hat b_2 + \hat b_3 + \hat b_4 + \hat b_5 + \hat b_6$. ```{r bc_weak, eval=FALSE} # Not run: exog="weak" adds the b2 dynamic term # (requires HAC SE for valid inference in dynamic models) fit_bc_weak <- ife_unbalanced(y ~ lag_y + x, data = df_dynamic, index = c("unit", "time"), r = 2, se = "hac", bias_corr = TRUE, exog = "weak") ``` **Summary of `exog` × `se` combinations:** | Regressor type | Errors | Recommended | Bias term | |---|---|---|---| | Strictly exogenous | i.i.d. | `se="standard"` or `"cluster"` | $b_3$–$b_6$ | | Strictly exogenous | Serially correlated | `se="hac"` | $b_3$–$b_6$ (HAC versions) | | Weakly exogenous (lag $y$) | i.i.d. | `se="cluster"` + `exog="weak"` | $b_2$–$b_6$ | | Weakly exogenous (lag $y$) | Serially correlated | `se="hac"` + `exog="weak"` | $b_2$–$b_6$ (HAC) | ## Balanced vs unbalanced on the same data As a sanity check, `ife_unbalanced()` on a complete (balanced) panel should agree with `ife(force="none")`: ```{r bal_vs_unbal} fit_bal <- ife(sales ~ price, data = cigar, index = c("state", "year"), r = 2, force = "none", se = "standard") fit_unb2 <- ife_unbalanced(sales ~ price, data = cigar, index = c("state", "year"), r = 2, se = "standard") cat(sprintf("ife (balanced): %.6f\n", fit_bal$coef["price"])) cat(sprintf("ife_unbalanced: %.6f\n", fit_unb2$coef["price"])) cat(sprintf("Difference: %.2e\n", abs(fit_bal$coef["price"] - fit_unb2$coef["price"]))) ``` The difference is at most $10^{-4}$ in finite samples (the EM inner loop has a convergence tolerance of $10^{-7}$). --- # Input validation and common errors Both `ife()` and `ife_unbalanced()` perform strict input validation and stop with an informative message if any check fails. | Check | Error message | |---|---| | `data` is not a `data.frame` | `"'data' must be a data.frame"` | | index columns not in `data` | `"'index' columns not found: ..."` | | NA in outcome or covariates | `"Missing values (NA) in variable '...'"` | | Duplicate (unit, time) pairs | `"Duplicate (unit, time) pairs found"` | | `r` < 1 (unbalanced) | `"'r' must be a positive integer"` | | Panel unbalanced in `ife()` | `"Panel is not balanced: ..."` | | `se` not in allowed set | `"'se' must be one of: ..."` | | `exog` not `"strict"`/`"weak"` | `"'exog' must be 'strict' or 'weak'"` | | `L_T < 1` | `"'L_T' must be a positive integer"` | --- # Quick reference: argument tables ## `ife()` — balanced panel | Argument | Default | Description | |---|---|---| | `formula` | — | `outcome ~ covariate1 + ...` | | `data` | — | balanced long-format `data.frame` | | `index` | — | `c("unit_col", "time_col")` | | `r` | `1L` | number of factors | | `force` | `"two-way"` | additive FE: `"none"` / `"unit"` / `"time"` / `"two-way"` | | `se` | `"standard"` | `"standard"` / `"robust"` / `"cluster"` | | `method` | `"static"` | `"static"` / `"dynamic"` (Moon & Weidner 2017) | | `bias_corr` | `FALSE` | apply analytical bias correction | | `M1` | `1L` | MA order for dynamic BC (used when `method="dynamic"`) | | `tol` | `1e-9` | outer-loop convergence tolerance | | `max_iter` | `10000L` | outer-loop iteration cap | ## `ife_unbalanced()` — unbalanced panel | Argument | Default | Description | |---|---|---| | `formula` | — | `outcome ~ covariate1 + ...` | | `data` | — | long-format `data.frame` (may be unbalanced) | | `index` | — | `c("unit_col", "time_col")` | | `r` | `1L` | number of factors | | `se` | `"standard"` | `"standard"` / `"robust"` / `"cluster"` / `"hac"` | | `init` | `"ols"` | initialisation: `"ols"` or `"nnr"` | | `bias_corr` | `FALSE` | apply SWW2025 analytical bias correction | | `exog` | `"strict"` | exogeneity: `"strict"` or `"weak"` | | `L_T` | `NULL` | HAC bandwidth ($\lfloor 2T^{1/5}\rfloor$ if `NULL`) | | `c_f` | `0.6` | SVT threshold constant (used with `ife_select_r_unb()`) | | `nu_NT` | `NULL` | NNR penalty (cross-validated if `NULL`) | | `tol` | `1e-9` | outer-loop convergence tolerance | | `max_iter` | `10000L` | outer-loop iteration cap | | `tol_em` | `1e-7` | EM inner-loop convergence tolerance | | `max_iter_em` | `500L` | EM inner-loop iteration cap | --- # References Bai, J. (2009). Panel data models with interactive fixed effects. *Econometrica*, **77**(4), 1229--1279. Bai, J. and Ng, S. (2002). Determining the number of factors in approximate factor models. *Econometrica*, **70**(1), 191--221. Baltagi, B. H. (1995). *Econometric Analysis of Panel Data*. Wiley. Moon, H. R. and Weidner, M. (2017). Dynamic linear panel regression models with interactive fixed effects. *Econometric Theory*, **33**, 158--195. Su, L., Wang, F. and Wang, Y. (2025). Estimation and inference for interactive fixed effects panel data models with unbalanced panels. SSRN Working Paper 5177283.