Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
128,332 characters · 20 sections · 122 citation commands
Bias and Consistency in Three-way Gravity Models
\vskip -.5cm We study the incidental parameter problem for the “three-way” Poisson {Pseudo-Maximum Likelihood} (“PPML”) estimator recently recommended for identifying the effects of trade policies {and in other panel data gravity settings}. Despite the number and variety of fixed effects involved, we confirm PPML is consistent for fixed $T$ and we show it is in fact the only estimator among a wide range of PML gravity estimators that is generally consistent in this context when $T$ is fixed. At the same time, asymptotic confidence intervals in fixed-$T$ panels are not correctly centered at the true point estimates, and cluster-robust variance estimates used to construct standard errors are generally biased as well. We characterize each of these biases analytically and show both numerically and empirically that they are salient even for real-data settings with a large number of countries. We also offer practical remedies that can be used to obtain more reliable inferences of the effects of trade policies and other time-varying gravity variables, {which we make available via an accompanying Stata package called \href{https://github.com/tomzylkin/ppml_fe_bias}{\tt{ppml_fe_bias}}}.
\enlargethispage{\baselineskip}
JEL Classification Codes: C13; C50; F10 \\ Keywords: Structural Gravity; Trade Agreements; Asymptotic Bias Correction \setcounter{page}{0} \thispagestyle{empty}
Despite intense and longstanding empirical interest, the effects of bilateral trade agreements on trade are still considered highly difficult to assess.
As emphasized in a recent practitioner's guide put out by the WTO (yotov_advanced_2016), many current estimates in the literature suffer from easily identifiable sources of bias (or “estimation challenges”). This is not for a lack of awareness. Papers showing leading causes of bias in the gravity equation are often among the most widely celebrated and cited in the trade field, if not in all of Economics.\footnote{For some context, if we start citation counts in 2003, anderson_gravity_2003 and santos_silva_log_2006 are, respectively, the most cited articles in the American Economic Review and the Review of Economics and Statistics. Paling only slightly in this exclusive company, baier_free_2007 is the 4th most-cited article in the Journal of International Economics, having gathered “only” 2,500 citations. Readers familiar with these other papers will also likely be familiar with helpman_estimating_2008's work on the selection process underlying zero trade flows, an issue we do not take up here.} In particular, it is now generally accepted that trade flows across different partners are interdependent via {the network structure of trade} (the main contribution of anderson_gravity_2003), that log-transforming the dependent variable is not innocuous (as argued by santos_silva_log_2006), and\textemdash most relevant to the context of trade agreements\textemdash that earlier, puzzlingly small estimates of the effects of free trade agreements were almost certainly biased downwards by treating them as exogenous (baier_free_2007).
As a consequence\textemdash and aided by some recent computational developments\textemdash researchers seeking to identify the effects of trade agreements have naturally moved towards more advanced estimation strategies that take on board all of the above concerns.\footnote{larch2019currency, ppmlhdfe, and stammann2017fast describe algorithms that enable fast estimation of the three-way models considered here. } In particular, a “three-way” fixed effects Poisson {Pseudo-Maximum Likelihood} (“FE-PPML”) estimator with time-varying exporter and importer fixed effects to account for {network dependence} and time-invariant exporter-importer (“pair”) fixed effects to address endogeneity has recently emerged as a logical workhorse method for empirical trade policy analysis.\footnote{Pair fixed effects are of course no substitute for good instruments. However, instruments for trade policy changes which are also exogenous to trade are understandably hard to come by. As discussed in head_gravity_2014's essential handbook chapter on gravity estimation, pair fixed effects have the advantage that the effects of trade agreements and other trade policies are identified from time-variation in trade within pairs. Causal interpretations follow if standard “parallel trend” assumptions are satisfied.} {It also has clear potential application to the study of network data more generally, such as data on urban commuting or migration (e.g., brinkman2019freeway; rothenberg2020; allen2018border; {beverelli2019migration)}.
However, one reason why {some} researchers may hesitate in embracing this estimator is the current lack of clarity regarding how the three fixed effects in the model may bias estimation, especially in the standard “fixed $T$” case where the number of time periods is small. Even though FE-PPML estimates can be shown to be asymptotically unbiased with a single fixed effect (a well-known result) as well as in a two-way setting where both dimensions of the panel become large (fernandez-val_individual_2016), the latter result does not come strictly as a generalization of the former one, leaving it potentially unclear whether a three-way model with a fixed time dimension should be expected to inherit the nice asymptotic properties of these other models.
Accordingly, the question we investigate in this paper is the extent to which the three-way FE-PPML estimator is affected by incidental parameter problems (IPPs). As is well known in both statistics and econometrics (neyman_consistent_1948,lancaster_orthogonal_2002), IPPs arise when estimation noise from estimates of fixed effects and other “incidental parameters” contaminates the scores of the main parameters of interest, inducing bias. In the worst case, this bias renders the estimates inconsistent, making estimation inadvisable. As we will show, while inconsistency is actually not a problem for three-way FE-PPML, both the estimated coefficients and standard errors are affected by meaningful biases due to IPPs that researchers should be aware of.
To state our main results more precisely, in gravity settings where the number of countries ($N$) goes to infinity and $T$ is small, we find the following: \vskip -4em
{To first explain our consistency results}, our basic strategy involves using the first-order conditions of FE-PPML to “{profile out}” (solve for) the pair fixed effect terms from the first-order conditions of the other parameters. Notably, this allows us to re-express the three-way gravity model as a two-way model {in which the only remaining incidental parameters are the exporter-time and importer-time fixed effects.} Three-way FE-PPML is therefore consistent in {fixed}-$T$ settings for largely the same reasons the two-way {models} considered in fernandez-val_individual_2016 are consistent, and we provide suitably modified versions of the regularity conditions and consistency results established by fernandez-val_individual_2016 for the simpler two-way case.
At the same time, it does not also follow that fernandez-val_individual_2016's earlier results for the asymptotic unbiased-ness of the two-way FE-PPML estimator similarly carry over to the three-way case when $T$ is fixed. {The key is that the resulting two-way estimator that is obtained after profiling out the fixed effects has its own special properties with respect to IPPs. When $T$ is fixed, the estimation noise in the remaining exporter-time and importer-time fixed effects induces an asymptotic bias of order $1/N$ as $N$ grows large, a result that is broadly consistent with most of the two-way settings studied in fernandez-val_individual_2016.} However, when instead both $N$ and $T$ grow large at the same rate, the estimator turns out to be unbiased asymptotically, analogous to what fernandez-val_individual_2016 found in the two-way FE-PPML case.
{The reason why asymptotic bias is a concern is that the asymptotic standard deviation is itself of order {$1/(N\sqrt{T})$}. Thus, when $T$ is {fixed}, both the bias in point estimates will be of comparable magnitude to their standard errors as $N \rightarrow \infty$, causing the asymptotic distribution of estimates to be incorrectly centered as discussed above. In practice, this is a less severe problem than inconsistency, but it does mean that standard hypothesis tests for assessing statistical significance are not reliable. One of the objectives of this paper will be to adapt some of the leading remedies from the recent literature on “large $T$” IPPs (see, e.g., arellano_understanding_2007) in order to re-center the asymptotic distribution of estimates and thereby restore asymptotically valid inferences.\footnote{The new literature on “large $T$” asymptotic bias in nonlinear FE {models} has emerged as a recent response to the well-known “fixed $T$” consistency problem first described in neyman_consistent_1948. Examples include phillips_linear_1999, hahn_asymptotically_2002, lancaster_orthogonal_2002, woutersen_robustness_2002, alvarez_time_2003, carro_estimating_2007, arellano_robust_2009, fernandez-val_bias_2011, and kato_asymptotics_2012. Unlike in most other settings explored in this literature, the panel estimator we consider is consistent regardless of $T$.}}
The bias in the estimated standard errors is similar to one that has been found in two-way gravity {settings} by several recent studies (egger_glm_2015,jochmans_two-way_2016,pfaffermayr2019gravity,pfaffermayr2021confidence). Intuitively, because the origin-time and destination-time fixed effects in the {model} each converge to their true values at a rate of only $1/\sqrt{N}$ (not $1/N$), the cluster-robust sandwich estimator for the variance has a leading bias of order $1/N$ (not $1/N^{2}$), and standard errors in turn have a bias of order $1/\sqrt{N}$. This latter type of bias is related to the general result that standard “heteroskedasticity-robust” variance estimators are downward-biased in small samples (see, e.g., mackinnon1985some,imbens2016robust), including for PML estimators (kauermann2001note), but is more severe in this setting due to an IPP. We should therefore be concerned that estimated confidence intervals may be too narrow in addition to being off-center.
For the bias in point estimates, we construct two-way analytical and jackknife bias corrections inspired by the corrections proposed in Fern\'andez-Val and Weidner fernandez-val_individual_2016,ARE. For the bias in standard errors, we show how kauermann2001note's method for correcting the PML sandwich estimator may be adapted to the case of a conditional estimator with multi-way fixed effects and cluster-robust standard errors. Our simulations confirm that these methods are usually effective at improving inferences. The jackknife correction reduces more of the bias in point estimates than the analytical correction in smaller samples, but the analytical correction does a better job at improving coverage, especially when also paired with corrected standard errors.
\enlargethispage{1em} For our empirical applications, we {first} estimate the average effects of a free trade agreement (FTA) on trade for a range of different industries using what would typically be considered a large trade data set, with 167 countries and 5 time periods. The biases we uncover vary in size across the different industries, but are generally large enough to indicate that our bias corrections should be worthwhile in most three-way gravity settings. For aggregate trade data (which yields results that are fairly representative), the estimated coefficient for FTA has an implied downward bias about 15%-22% of the estimated standard error, and the implied downward bias in the standard error itself is about 11% of the original standard error. As a means of further demonstration, we also apply our corrections to replication data from several recent papers that have used three-way gravity {models}. This latter exercise reveals several instances in which our methods make a material difference for assessing statistical significance. It also highlights the possibility that the bias in standard errors can sometimes be severe, as much as 40% or more in some cases.
Aside from fernandez-val_individual_2016's work on two-way nonlinear {models}, pesaran_estimation_2006, bai_panel_2009, hahn_reducing_2006, and moon_dynamic_2017 have each conducted similar analyses for two-way linear {models} with interacted individual and time fixed effects. Turning to three-way models, hinzetal have recently developed bias corrections for dynamic three-way probit and logit {models} based on asymptotics suggested by ARE where all three panel dimensions grow at the same rate. Though widely applicable, this approach is not appropriate for our setting because of the different role played by the time dimension when the estimator is FE-PPML.\footnote{Also related are the GMM-based differencing strategies for two-way FE {models} proposed by charbonneau_multiple_2012 and jochmans_two-way_2016. These strategies rely on differencing the data in such as way that the resulting GMM moments do not depend on any of the incidental parameters. In principle, these methods could be extended to allow for differencing across a time dimension as well in a three-way panel.} In the network context, graham2017econometric, dzemski2018empirical, and chen2014nonlinear have studied asymptotic bias in network {models} with node-specific (possibly sender- and receiver-specific) fixed effects. chen2014nonlinear's analysis is espeically notable in that they allow these node-specific effects to be vectors rather than scalars, similar to the exporter-time and importer-time fixed effects that feature in gravity {models}. Our bias expansions substantially differ from those of chen2014nonlinear because the equivalent outcome variable in our setting (trade flows observed over time for a given pair) is also a vector rather than a scalar and because we work with a conditional moment {model} where the distribution of the outcome may be misspecified.
In what follows, Section (ref) first provides a discussion of why IPPs are a concern for gravity {models} and of the no-IPP properties of FE-PPML. Section (ref) then establishes bias and consistency results for the three-way gravity {model} specifically and discusses how to implement bias corrections. Sections (ref) and (ref) respectively present simulation evidence and empirical applications. Section (ref) concludes, and an Appendix adds further simulation results and technical details, including proofs.
Gravity models are now routinely estimated using FE-PPML with multiple sets of fixed effects. As we discuss in this section, these practices follow naturally from the gravity model's theoretical microfoundations but are not without need for further scrutiny. In particular, because the underlying model is nonlinear, it is important to clarify that, while PPML is known to be free from incidental parameter bias in some special cases, it is by no means immune to IPPs in general. It will also be useful for us to provide some general discussion of IPPs and the different ways in which they may manifest.
As documented in head_gravity_2014, the emergence of rich and varied theoretical foundations for the gravity equation has fueled a “fixed effects revolution” in the gravity literature over the last two decades. As such, we find it useful to briefly describe a simple trade model and discuss how it may be used to motivate an estimating equation with either two-way or three-way fixed effects.
To establish some notation we will use throughout the paper, we will consider a world with $N$ countries and we will let $i$ and $j$ respectively be indices for exporter and importer. For now, we will focus on deriving a two-way gravity model where the two fixed effects account for each country's multilateral resistance. Later, we will add a time dimension and a third fixed effect that absorbs all time-invariant components of trade costs.
To add some theoretical structure, suppose that trade flows are given by the following gravity equation:
Here, $y_{i}:=\sum_{j}y_{ij}$ and $y_{j}:=\sum_{i}y_{ij}$ are the market sizes of the two countries, $\tau_{ij}\ge1$ is a bilateral trade cost, $\theta>0$ is the trade elasticity, and $\Pi_{i}$ and $P_{j}$ respectively are the outward and inward multilateral resistances from anderson_gravity_2003, which capture how bilateral trade flows depend on each country's opportunities for trade with third countries. More formally, these latter terms are derived from the following two relationships that are inherent to all general equilibrium gravity models:
As shown, these terms respectively aggregate the exporter's ability to export goods to more desirable import markets and the importer's ability to import from more capable exporters.\footnote{This presentation of the gravity model readily conforms to the trade models used in eaton_technology_2002 or anderson_gravity_2003, though the interpretation of $\theta$ differs across the two models. With some minor modifications, this setup can also be made compatible with any of the theoretical gravity models considered in head_gravity_2014 or costinot_trade_2014. To be clear, our econometric results do not require any particular microfoundation for the gravity equation.}
For estimation, it is typical to parameterize the trade cost $\tau_{ij}$ as depending exponentially on some variables of interest, i.e.,
where $x_{ij}$ are the components of trade costs whose effects we wish to estimate. Because not all trade costs are reflected in $x_{ij}$, we also allow for “unobserved” trade costs via the {idiosyncratic} trade cost term $\omega_{ij}$. Combining (ref) with (ref) then delivers the following estimating equation:
where $\alpha_{i}=\ln(y_{i}/\Pi_{i}^{-\theta})$ and $\gamma_{j}=\ln(y_{j}/P_{j}^{-\theta})$ are origin and destination fixed effects that absorb market sizes and multilateral resistances and $\omega_{ij}$ now provides a multiplicative error term.\footnote{Alternatively, it is sometimes common to write trade costs as a log-linear function, i.e, $\ln\tau_{ij}=x_{ij}\beta+e_{ij}$, with $e_{ij}$ now reflecting unobserved (log) trade costs. Interestingly, these two ways of specifying the error term do not necessarily have equivalent implications for estimation. In the log-linear formulation, if the log-error term $e_{ij}$ is assumed to be heteroskedastic with mean zero, Jensen's inequality implies that estimation in levels will be biased and log-OLS will be consistent.} When we introduce the three-way model, all of the terms shown in (ref) will have a further subscript for time, and the the unobserved trade cost will have a time-invariant component that will motivate the use of an added $ij$ fixed effect.
To motivate the arc of the rest of the paper, several points stand out from the estimation suggested by (ref). First, the implied moment condition for estimation is
{which follows after imposing $\mathbb{E}(\omega_{ij}|x_{ij},\!\alpha_i,\!\gamma_j)=1.$}\footnote{Note that consistent estimation of $\beta$ actually does not require $\mathbb{E}(\omega_{ij}|\cdot)=1$ in this case but rather $\mathbb{E}(\omega_{ij}|\cdot)=\widetilde{\omega}_{i}\widetilde{\omega}_{j}$, where $\widetilde{\omega}_{i}$ and $\widetilde{\omega}_{j}$ could be country-specific components of unobserved trade costs that would be absorbed by the fixed effects. For the three-way model, one requires $\mathbb{E}(\omega_{ijt}|\cdot)=\widetilde{\omega}_{it}\widetilde{\omega}_{jt}\widetilde{\omega}_{ij}$.} As discussed in santos_silva_log_2006, consistent estimation of the trade cost parameters in $\beta$ therefore generally requires a nonlinear model. {Second, because unobserved trade costs enter the country-specific terms $\alpha_{i}$ and $\gamma_{j}$ through the system of multilateral resistances, we treat $\alpha_{i}$ and $\gamma_{j}$ as unknown parameters that will be noisily estimated, raising concerns about a possible IPP}.\footnote{ {As demonstrated in pfaffermayr2021confidence, if we assume (i) there are no unobserved trade costs (such that $\omega_{ij}$ does not enter the system in (ref)) and (ii) the aggregate quantities $y_i$ and $y_j$ are perfectly observed, then it is better to regard $\alpha_i$ and $\gamma_j$ as reflecting constraints rather than as incidental parameters. Since fally_structural_2015 shows that constrained PPML and FE-PPML produce the same estimates in this context, there is no concern about IPPs if these assumptions are met.}} As we go on to discuss, the FE-PPML estimator that is most often used in this context has some special robustness against IPPs, but this robustness does not hold for FE-PPML in general, especially once we deviate from the two-way gravity setting implied by (ref).
In the context of fixed effects models, IPPs occur when the estimation noise in the fixed effects contaminates the scores of the other parameters being estimated, inducing a bias. This bias can manifest in a variety of different ways; thus it is useful to provide a generic characterization that can illustrate the different possibilities that may arise. To that end, let $n$ be the total number of observations and let $p$ be the total number of parameters being estimated, inclusive of any fixed effects. As described in ARE, what we need to be concerned with the number of observations that are available to estimate each fixed effect, i.e., $\ensuremath{n/p}$. More precisely, when appropriate regularity conditions are satisfied, the estimated $\widehat{\beta}$ may generally be thought of as having the following bias and standard deviation:
where $b\in\mathbb{R}$ and $c>0$ are constants that depend on the model being estimated. {As this presentation emphasizes, {all} estimators in nonlinear settings are generally biased in small samples, but this is not the same thing as saying that the bias always poses a problem for inferring statistical significance. As we will discuss, in larger samples, what matters is whether the bias disappears faster than the standard error as $n\rightarrow \infty.$\footnote{In the panel data literature, these results for the bias and standard deviation are usually derived not for $\widehat{\beta}$ directly, but for the asymptotic distribution of $\widehat{\beta}$ (because there are cases where $\widehat{\beta}$ may not have a first or second moment, but nevertheless has a well-defined limiting distribution with finite moments). We ignore this distinction for our heuristic discussion here.}
To provide a simple taxonomy of the cases that can arise, consider first the standard textbook treatment of maximum likelihood estimation, where we usually have $p$ fixed while $n\rightarrow\infty$. In this case, the bias in $\widehat{\beta}$ becomes asymptotically negligible as compared to the standard deviation, which crucially means that estimated confidence intervals can be expected to be centered at the truth when the data becomes sufficiently large. By contrast, in the classical IPP of neyman_consistent_1948, the number of parameters grows at the same rate as the number of observations, implying that the bias does not converge to zero asymptotically. In that case, the fixed effect estimator is inconsistent.}
The gravity model with two-way fixed effects then serves to illustrate a third possibility that will also be applicable to the results that follow for three-way gravity models. In the two-way gravity setting, $p$ is on the order of $2N$, where $N$ is the number of countries, and $n$ is on the order of $N^{2}$. Consequently, the bias and standard deviation of $\widehat{\beta}$ are given by
In this case, as $N\rightarrow\infty$, the estimated $\widehat{\beta}$ is consistent (both the standard deviation and bias converge to zero), but we also have \[ \lim_{N\rightarrow\infty}\;\frac{{\rm bias}(\widehat{\beta})}{{\rm std}(\widehat{\beta})}=\frac{2\,b}{c}. \] {As we discuss below, the two-way PPML gravity {estimator} is a special case where we actually have that $b=0$. However, for other two-way gravity {estimators} (such as two-way Gamma PML for example), we generally have that $b\neq0$, meaning the bias will not disappear relative to the standard error as $N\rightarrow\infty$.} Compared to neyman_consistent_1948, the IPP these estimators suffer from is not an inconsistency problem but rather an asymptotic bias problem, whereby the slow convergence of the fixed effects causes the {asymptotic distribution} for $\widehat{\beta}$ to be incorrectly centered as it converges to the truth.\footnote{Consistency here follows from how the number of fixed effects grows only with the square root of the sample size, as discussed in egger_trade_2011. The bias in the asymptotic distribution for two-way {models} was proven by fernandez-val_individual_2016, discussed below.} This version of the IPP is more benign, but ignoring the bias will nonetheless result in invalid inferences and test results. The “large $T$” panel data literature therefore discusses various methods for bias correction of $\widehat{\beta}$ that restore asymptotically valid inference. Importantly, the degree to which inferences are biased depends on the bias constant $b$, which {cannot be} known beforehand without applying such a correction.
{To provide a more visual illustration of these ideas, Figure (ref) presents simulation results for the three cases we have just discussed: inconsistency (top-left), no asymptotic bias (top-right), and asymptotic bias (bottom-left). Since asymptotic bias will ultimately be our focus, it is worth noting from the figure how the estimates are consistent in this case\textemdash the distribution will collapse to the true value as $N\rightarrow\infty$\textemdash but confidence bounds based on these estimates will clearly be inappropriate.}
Our discussion of IPPs thus far has been for the generic estimation of a nonlinear model. However, the PPML estimator that is most commonly used to estimate gravity models actually behaves very differently than other estimators in this context. In the classic panel data setting with “one way” fixed effects, for example, FE-PPML has the very special property that the IPP bias constant $b$ turns out to be zero, meaning that it is {asymptotically unbiased} in situations where other estimators tend to be inconsistent. As discussed in wooldridge_distribution-free_1999, the reason behind this result is that {the same estimator can be obtained from a multinomial model} that does not depend on the fixed effects.\footnote{The earliest references to present versions of this result include andersen1970asymptotic, palmgren_fisher_1981, and hausman_econometric_1984. wooldridge_distribution-free_1999's contribution is to show that FE-PPML is consistent even when the assumed distribution of the data is misspecified. Our Lemma (ref) in the Appendix clarifies that FE-PPML is relatively unique in this regard versus similar models.}
This special property of FE-PPML has important implications for estimating gravity models as well. For two-way gravity settings, the asymptotic bias of $\widehat{\beta}$ when $N\rightarrow\infty$ was worked out in fernandez-val_individual_2016. They show that the two IPP contributions from $\alpha_{i}$ and $\gamma_{j}$ “decouple” asymptotically, such that the overall bias can be decomposed as the sum of two bias terms that would be expected in a one-way setting, i.e., $b_{(\alpha)}/N+b_{(\gamma)}/N$, {where the $N$'s come from the number of observations associated with each fixed effect.} Because $b=0$ for the one-way FE-PPML case, we also have $b_{(\alpha)}=b_{(\gamma)}=0$ in the two-way setting; that is, two-way FE-PPML gravity estimates for $\widehat{\beta}$ are asymptotically unbiased just as one-way FE-PPML estimates are.\footnote{Note that Theorem 4.1 in fernandez-val_individual_2016 is written for the correctly specified case, where $y_{ij}$ is actually Poisson distributed. However, Remark 3 in the paper gives the extension to conditional moment models, where for the FE-PPML case only the moment condition in (ref) needs to hold. Their paper considers standard panel models, as opposed to trade models, but the only technical difference is that $y_{ij}$ is often not observed for the trade model when $i\neq j$. This missing diagonal has no meaningful effect on any of the results we discuss.}
Taken together, these results might create the impression that FE-PPML is generally immune to IPPs, regardless of what fixed effects are included in the model. Thus, it is important to clarify that a key feature of the two-way gravity model is that both fixed effect dimensions grow only with the square root of the panel size, such that the estimation noise in the estimated $\widehat{\alpha}_{i}$'s and $\widehat{\gamma}_{j}$'s disappears asymptotically. {As we discuss in the Appendix, if we instead consider a model where both fixed effects grow with $n$ rather than with its square root, the IPPs associated with each fixed effect do not decouple from one another, and FE-PPML in this case is actually inconsistent.}
To synthesize these points, FE-PPML has a very special property\textemdash one can condition out one of the fixed effects\textemdash but this property has an important limitation\textemdash the resulting multinomial {model} does not inherit the same no-bias properties as the original PPML estimator with respect to any further fixed effects. Both of these results will be fleshed out in more detail in the following section when we recast the three-way gravity model as a {two-way multinomial model} in order to obtain an appropriate expression for the bias. Doing so will also allow us to highlight another reason why FE-PPML is not immune to IPPs: even for the two-way gravity model, while the $\alpha_{i}$ and $\gamma_{j}$ parameters do not induce an IPP bias in $\widehat{\beta}$, they nonetheless have implications for the estimated variance that are not innocuous; we thus will devote attention to this issue as well.
To recap the sequence of results just described, we know that FE-PPML estimates with one fixed effect do not suffer from an IPP. We also know that FE-PPML may have an IPP in models with more than one fixed effect, but it is both consistent and asymptotically unbiased in two-way gravity settings where neither fixed effect dimension grows at the same rate as the size of the panel. As we will now show, each of these earlier results will be useful for understanding the more complex case of a three-way gravity model that adds a time dimension and a third set of fixed effects to the above two-way model. We also describe a series of bias corrections for the three-way model, including for the possible downward bias of the estimated standard errors.
{To formally introduce the three-way model, we add an explicit time subscript $t\in\{1,\ldots,T\}$ to $y_{ij}$, $x_{ij}$, and $\omega_{ij}$ from the prior model and add a “country-pair”-specific fixed effect $\eta_{ij}$, such that trade costs are now given by $\tau_{ijt}^{-\theta}=e^{x_{ijt}'\beta+\eta_{ij}}\omega_{ijt}$. All other elements in the original trade model likewise acquire a time subscript, meaning that $\alpha_{it}=\ln{y_{it} / \Pi_{it}^{-\theta}}$ and $\gamma_{jt}=\ln{y_{jt}/P_{jt}^{-\theta}}$ also must be indexed by $t$.} The model now reads
where the three fixed effects now respectively index exporter-time, importer-time, and country-pair.\footnote{{Note that a multiplicative error term is not necessary to deliver the moment condition in (ref). {We could instead have an additive error term $\varepsilon_{ijt}=y_{ijt}-\lambda_{ijt}$. In this case}, it would be more natural to think of it as coming from measurement error. In addition, note that we assume the true model is as written in (ref) and assume away, e.g., any unobserved heterogeneity in $\beta$. Allowing for this type of heterogeneity is an important extension for future work to address.}} The unobserved trade cost $\omega_{ijt}\ge0$ continues to serve as an {error term}, such that $y_{ijt}=\lambda_{ijt}\omega_{ijt}\ge0$. {We thus allow for zero trade flows}. For the asymptotics using the three-way model, we consider $T$ fixed, while $N\rightarrow\infty$. The FE-PPML estimator maximizes
over $\beta$, $\alpha$, $\gamma$ and $\eta$.
{Our strategy for showing the consistency of this estimator will capitalize on the special properties of FE-PPML discussed in the previous section. In particular, we will exploit the fact} that not all of the fixed effect dimensions grow at the same rate as $N$ increases. The numbers of exporter-time and importer-time fixed effects each increase with $N$ (as before), but the dimension of {the pair fixed effect} $\eta$ increases with $N^2$, since adding another country to the data adds another $N-1$ {pairs} to the estimation. It therefore makes sense to first “profile out” (i.e., solve for) $\eta$ {so that we may deal with the remaining two fixed effects in turn}. For given values of $\beta$, $\alpha$, $\gamma$, {maximizing over $\eta$ gives us}
We therefore have
with {
thus leaving us with the likelihood of a multinomial model where the only incidental parameters are $\alpha_{it}$ and $\gamma_{jt}$. } {Using (ref), one can easily verify that there is no bias in the score of the profile log-likelihood $ \ell_{ij}(\beta, \alpha_{it}, \gamma_{jt}) $ when evaluated at the true parameters $\beta^0$, $\alpha_{it}^0$, and $\gamma_{jt}^0$.} {The reason for this is exactly the same as for the classic panel data setting discussed above.} {Furthermore, the remaining fixed effects $\alpha_{it}$ and $\gamma_{jt}$ grow only with the square root of the sample size as $N\rightarrow\infty$, implying that they are consistently estimated.} {This in turn leads us to the following result}:
{{{Intuitively, this result follows because of how the special properties of FE-PPML allow us to rewrite the three-way gravity model as a two-way model without introducing a $1/T$ bias. The form of the bias in $\widehat{\beta}$ is therefore the same as in (ref), such that three-way FE-PPML is consistent as $N\rightarrow\infty$ largely for the same reason two-way FE-PPML and other two-way PML gravity estimators are generally consistent. However, in the context of three-way estimators, we can also state a stronger result that applies more narrowly to FE-PPML in particular:
{In other words, three-way FE-PPML is unique among three-way PML estimators} in that its consistency does not require strong assumptions about the conditional variance of $y_{ijt}$. To draw an appropriate contrast, it is possible to obtain a closed form solution for the pair fixed effect $\widehat{\eta}_{ij}$ so long as $g(\widehat{\lambda}_{ijt})$ is of the form $g(\widehat{\lambda}_{ijt})=\widehat{\lambda}_{ijt}^q$, where $q$ can be any real number. Notably, this latter class of estimators not only includes FE-PPML (for which $q=0$), but also includes other popular gravity estimators such as Gamma PML ($q=-1$) and Gaussian PML ($q=1$). However, as we discuss in the Appendix, these other estimators are only consistent if the conditional variance is proportional to ${\lambda}_{ijt}^{1-q}$, in which case they inherit the properties of their associated MLE estimators.}
Because three-way FE-PPML inherits the consistency properties of the two-way estimator, one might expect that it also inherits its {“no asymptotic bias” properties} as well. However, this is where the limitations of PPML's no-IPP properties become apparent. While the profile log-likelihood in (ref) is now of a similar form to the two-way {models} considered in fernandez-val_individual_2016, notice that it no longer resembles the original FE-PPML log-likelihood. Their no-bias result for two-way FE-PPML therefore does not carry over to {the three-way model}, and it is possible to show that FE-PPML {estimates have} an asymptotic bias in this setting.
{Before proceeding, it is helpful to first revisit the intuition established in Section (ref) that shapes how we expect the bias to behave. As we have discussed, in a model with $p$ parameters and $n$ observations, the bias should be proportional to $p/n$, whereas the standard error should vary with $1/\sqrt{n}$. After profiling out $\eta$, the resulting two-way model has $\sim2NT$ parameters vs. $\sim N^{2}T$ observations. If $T$ is held fixed, we would expect the bias and the standard error to decrease at the same rate ($1/N$), raising concerns about a possible asymptotic bias problem and guiding us as to its form. Readers should keep this intuition in mind in reading through the technical details that follow.}
{To illustrate more precisely where the bias comes from}, {it is necessary to examine how the estimated fixed effects enter the score for ${\beta}$ using a Taylor expansion. To that end, first let $\phi:={\rm vec}(\alpha,\gamma)$ be a vector that collects all of the exporter-time and importer-time fixed effects, such that we can rewrite $\ell_{ij}$ slightly as $\ell_{ij}=\ell_{ij}(\beta,\phi)$. We can then similarly define the function $\widehat{\phi}(\beta)$ as collecting the estimated fixed effects $\widehat{\alpha}$ and $\widehat{\gamma}$ as functions of $\beta$. Next, we construct a second-order expansion of the score for ${\beta}$ around the true {set of fixed effects} $\phi^{0}$ and evaluated at the true parameter $\beta^{0}$:
}This expression is near-identical to a similar expansion that appears in fernandez-val_individual_2016\textemdash differing mainly in that $\ell_{ij}$ is a vector rather than a scalar\textemdash and communicates the same essential insights: because the latter two terms in (ref) are generally not equal to zero, the score for ${\beta}$ is biased, with the bias depending on the interaction between the higher-order partial derivatives of $\ell_{ij}$ and the estimation errors in $\widehat{\alpha}_{i}$ and $\widehat{\gamma}_{j}$ as well as their variances and covariances.
Demonstrating the bias in {this particular setting} then requires that we introduce some additional notation, mainly to provide some shorthand for the higher-order partial derivatives of $\ell_{ij}$ that appear in (ref). To do so, we first find it convenient to let $\vartheta_{ijt}:=\lambda_{ijt}/\sum_{\tau}\lambda_{ij\tau}$. {We then define the $T\times1$ “score” vector $S_{ij}$, the $T \times T$ “Hessian” matrix $H_{ij}$ and the $T\times T\times T$ cubic tensor $G_{ij}$ (the “third partial”), with their respective elements given by
where it should be understood that all of these terms are evaluated at the true values for all parameters.} Explicit formulas for $G_{ij,tsr}$ are provided in the Appendix.
{The value of defining these objects is that they allow us to easily form terms identified by (ref) as being important for the bias of the score. For example, $S_{ij}$ allows us to obtain $\partial\ell_{ij}/\partial\beta^{k}={x_{ij,k}^{\prime}S_{ij}}$. Likewise, we also have that $\partial^{2}\ell_{ij}/\partial\alpha_{i}\partial\beta^{k}=\partial^{2}\ell_{ij}/\partial\gamma_{j}\partial\beta^{k}=-H_{ij}x_{ij,k}$ and that \[ \frac{\partial^{3}\ell_{ij}}{\partial\alpha_{i}\partial\alpha_{i}^{\prime}\partial\beta^{k}}=\frac{\partial^{3}\ell_{ij}}{\partial\alpha_{i}\partial\gamma_{j}^{\prime}\partial\beta^{k}}=\frac{\partial^{3}\ell_{ij}}{\partial\gamma_{j}\partial\alpha_{i}^{\prime}\partial\beta^{k}}=\frac{\partial^{3}\ell_{ij}}{\partial\gamma_{j}\partial\gamma_{j}^{\prime}\partial\beta^{k}}=G_{ij}x_{ij,k}, \] where we use the convention that $G_{ij}x_{ij,k}$ is a $T\times T$ matrix with elements $[G_{ij}x_{ij,k}]_{st}=\sum_{r}G_{ijrst}x_{ijr,k}$.} We also find it useful to define the expected Hessian $\bar{H}_{ij}=\mathbb{E} (H_{ij} \, |\,x_{ij} )$ and, similarly, the expected third partial $\bar{G}_{ij}=\mathbb{E} (G_{ij} \,|\,x_{ij} )$.\footnote{ {Because $\bar{H}_{ij}$ is only positive semi-definite (not positive definite), we use a Moore-Penrose pseudoinverse whenever the analysis requires we work with an inverse of $\bar{H}_{ij}$. Specifically, we have that $\bar{H}_{ij}\,\iota_{T}=0$, where $\iota_{T}=(1,\ldots,1)^{\prime}$ is a T-vector of ones. Thus, $\bar{H}_{ij}$ is only of rank $T-1$ rather than of rank $T$.} } Finally, we define the $K$-vector $\widetilde{x}_{ij}$ as an appropriate two-way within-transformation of $x_{ij}$ {that purges it of the fixed effects}; see the Appendix for details.}
{Obtaining a tractable expression for the bias then involves following the logic of (ref) and plugging in the just-defined objects $S_{ij}$, $H_{ij}$, $G_{ij}$, and $\widetilde{x}_{ij}$ where appropriate. Before doing so, we invoke the assumption that observations are serially correlated within pairs but independent across pairs, as is commonly assumed in the literature (see yotov_advanced_2016.) This assumption turns out to cause the IPPs associated with $\alpha_i$ and $\gamma_j$ to “decouple”,\footnote{ In particular, all elements of the cross-partial objects $\mathbb{E}[\partial^{2}\ell_{ij}/\partial\alpha_i \partial\gamma_{j}]$, $\mathbb{E}[\partial^{3}\ell_{ij}/\partial\alpha_i \alpha_i^{\prime} \partial\gamma_{j}]$, etc. can be shown to be asymptotically small for $N\rightarrow \infty$. Thus, in what follows, $B_N$ reflects the contribution of the $\alpha_i$ parameters to the bias and $D_N$ reflects the contribution of the $\gamma_j$ parameters. {As we discuss in the Appendix, relaxing this assumption can change the expression of the bias.}} leading to the following proposition:}
The above proposition establishes the asymptotic distribution of the three-way gravity estimator as $N\rightarrow \infty$, including the asymptotic bias $(N-1)^{-1}W_{N}^{-1}(B_{N}+D_{N})$. {Intuitively, this bias can be decomposed as the product of the {inverse expected Hessian with respect to $\beta$} (i.e.\ $W_{N}^{-1}$), the rate of asymptotic convergence (essentially $1/N$), and the bias of the score from (ref), which here is given by the combined term $B_{N}+D_{N}$. $B_N$ reflects the contribution to the bias from the noise in $\widehat \alpha_i$, whereas $D_N$ reflects the contribution from the noise in $\widehat \gamma_j$. {The first terms in both $B_N$ and $D_N$ come from the second term in (ref), reflecting the estimation error in the estimated fixed effects, and the second terms in $B_N$ and $D_N$ echo the third term in (ref), reflecting their variance.}}
Thus, in the end, the bias reduces to essentially the same simple formula we gave in Section (ref) for the bias of the two-way gravity model, i.e., \[ \frac{1}{N-1}\,b_{(\alpha)}+\frac{1}{N-1}\,b_{(\gamma)}, \] where $b_{(\alpha)}=W_{N}^{-1} B_{N}$ and $b_{(\gamma)}=W_{N}^{-1} D_{N}$ are constants that do not vary with $N$. Importantly, and unlike in the two-way FE-PPML setting, the three-way model does not give us the no-bias result that $B_{N}=D_{N}=0$, as the following discussion helps to illustrate.
Admittedly, the complexity of the objects that appear in Proposition (ref) may make it difficult to appreciate the general point that the three-way estimator is not {asymptotically unbiased}. One way to make these details more transparent is to focus our attention on the simplest possible panel model where $T=2$. The convenient thing about this simplified setting is that the likelihood function $\ell_{ij}$ can be reduced to just a scalar: $\ell_{ij}=y_{ij1}\log\vartheta_{ij1}+y_{ij2}\log\left(1-\vartheta_{ij1}\right)$, where now {
with $\Delta x_{ij}=x_{ij1}-x_{ij2}$, $\Delta \alpha_{i}= \alpha_{i1}- \alpha_{i2}$, and $\Delta \gamma_{j}= \gamma_{j1}- \gamma_{j2}$. {These normalizations allow us to express all of the objects that appear in Proposition (ref) as also just scalars, and we can therefore easily derive the following result:} }
{Two points then stand out based on the above expression: (i) all bias terms in $B_{N}^{k}$ and $D_{N}^{k}$ generally depend on the distribution of $y_{ij}$ and do not depend on it in the same way; (ii) none of these terms generally equals 0. The first of these two observations can be seen from how the bias depends on the expected second moments of $y_{ij}$ (e.g., $\mathbb{E}(y_{ij1}^{2})$, $\mathbb{E}(y_{ij1}y_{ij2})$, etc.), marking an important difference from the models that were considered in fernandez-val_individual_2016.\footnote{The specific examples they use are the Poisson model, which is unbiased, and the probit model, which requires the distribution of $y_{ij}$ to be correctly specified. They also provide a bias expansion for “conditional moment” models that allow the distribution of $y_{ij}$ to be misspecified. Beyond this theoretical discussion, bias corrections for misspecified models have yet to receive much attention, however.} Among other things, the difficulty associated with estimating these second moments means that analytical bias corrections may not necessarily offer superior performance to distribution-free method such as the jackknife. The second observation mainly follows from the first. It can also be shown that $\sum_{j\neq i}\bar{G}_{ij}\Delta\widetilde{x}_{ij}\neq0$, which ensures that the second term can never be zero.}
While Proposition (ref) only focuses on asymptotics where $N\rightarrow\infty$, the three-way gravity panel also features a time dimension ($T$), and it is interesting to wonder how the above results may depend on changes in $T$. {As we show in the Appendix}, { large $T$ only makes a difference for the asymptotic order of the bias of $\widehat \beta$ if there is only weak time dependence between observations belonging to the same country pair, in the sense described by hansen2007asymptotic.\footnote{By “weak” time dependence, we mean that any such dependence dissipates as the temporal distance between observations increases. Alternatively, if observations are correlated regardless of how far apart they are in time, the standard error is always of order $1/N$ (see hansen2007asymptotic), and the same will also be true for the asymptotic bias. The latter is arguably a less natural assumption in this context, however.} We will henceforth assume any time dependence is weak. The following remark then describes some additional asymptotic results for when $T$ is large. }
To elaborate further, letting $T\rightarrow\infty$ is obviously not sufficient for either $\alpha$ or $\gamma$ to be consistently estimated and does not solve the IPP, as stated in part (i). However, as part (ii) tells us, $T$ still plays an interesting role in conditioning the bias when both $N$ and $T$ jointly become large. Intuitively, because $W_{N}^{-1}$ is of order $1/T$ as $T\rightarrow\infty$, {whereas $B_{N}$ and $D_{N}$ are both of order 1}, the bias in $\widehat{\beta}$ effectively vanishes at a rate of $1/(NT)$ as both $N,T\rightarrow \infty$, such that it disappears asymptotically in relation to the order-$1/(N\sqrt{T})$ standard error. However, since $T$ is usually small relative to $N$ in this context, it remains to be seen whether these asymptotic results carry over to practical settings.
Of course, even if the point estimates are correctly centered, inferences will still be unreliable if the estimates of the variance used to construct confidence intervals are not themselves unbiased. For PPML, confidence intervals are typically obtained using a “sandwich” estimator for the variance that accounts for the possible misspecification of the model. However, as shown by kauermann2001note, the PPML sandwich estimator is generally downward-biased in finite samples. Furthermore, for gravity models (both two-way and three-way), the bias in the sandwich estimator can itself be formalized as a kind of IPP.\footnote{This type of IPP has similar origins to the one described in verdier2018, who considers a dyadic linear model with two-way FEs and sparse matching between the two panel dimensions.}
To illustrate the bias of the sandwich estimator in our three-way setting, recall that we can express the variance of $\widehat{\beta}$ as ${\rm Var}(\widehat{\beta}-\beta)=N^{-1}(N-1)^{-1}W_{N}^{-1}\Omega_{N}W_{N}^{-1}$. As is also true for the linear model (cf., mackinnon1985some,imbens2016robust), {the bias arises because plugin estimates for the “meat” of the sandwich $\Omega_{N}$ depend on the estimated score variance $\mathbb{E}(\widehat{S}_{ij}\widehat{S}_{ij}^{\prime})$ rather than on the true variance $\mathbb{E}(S_{ij}S_{ij}^{\prime})$.} Even though $\mathbb{E}(\widehat{S}_{ij}\widehat{S}_{ij}^{\prime})$ is a consistent estimate for $\mathbb{E}(S_{ij}S_{ij}^{\prime})$, it will generally be downward-biased in finite samples. Notably, this bias may be especially slow to vanish for models with gravity-like fixed effects.
To see this, we follow the same approach as kauermann2001note. Specifically, we use the special case where $\mathbb{E}(S_{ij}S_{ij}^{\prime})=\kappa\bar{H}_{ij}$ (such that $\Omega_{N}=\kappa W_{N}$, meaning PPML is correctly specified) to demonstrate that $\mathbb{E}(\widehat{S}_{ij}\widehat{S}_{ij}^{\prime})$ generally has a downward bias. Under this assumption, it is possible to show that the expected outer product of the fitted score $\mathbb{E}(\widehat{S}_{ij}\widehat{S}_{ij}^{\prime})$ has a first-order bias of}
where $W_{N}^{(\phi)}:=\mathbb{E}_{N}[-\partial^{2}\ell_{ij}/\partial\phi\partial\phi^{\prime}]$ captures the expected Hessian of the concentrated likelihood with respect to $\alpha$ and $\gamma$ and where $d_{ij}$ is a $T\times dim(\phi)$ matrix of dummies such that each row satisfies $d_{ijt}\phi=\alpha_{it}+\gamma_{jt}$.\footnote{A detailed derivation of (ref) is provided in the Appendix.}
The two terms on the right-hand side of (ref) are both negative definite, implying that the sandwich estimator is generally downward-biased\textemdash and definitively so if the model is correctly specified. {Since we work with cluster-robust standard errors, a relevant comparison to draw here is with cameron2008bootstrap, who have previously shown that the cluster-robust sandwich estimator has a downward bias that depends on the number of clusters. In our setting, the standard cameron2008bootstrap bias is reflected in the first term on the righthand-side of (ref), which captures how the the bias depends on the variance of $\widehat{\beta}$. The second term, which arises because of an IPP, captures how much of the bias is due to the variance in the estimated origin-time and destination-time fixed effects in $\widehat{\phi}$. The former term decreases with $1/N^{2}$---i.e., with the number of pairs/clusters---but the latter term only decreases with $1/N$, since increasing $N$ by $1$ only adds $1$ additional observation of each origin-time and destination-time fixed effect.\footnote{ The $[N (N-1)]^{-1} W_{N}^{(\phi)-1}$ matrix that appears in the second term is the inverse Hessian with respect to the fixed effects and thus reflects their variance. Because adding a new country only adds one new observation of each fixed effect, the diagonal elements of this matrix decrease with only $1/N$ as $N$ increases, despite how the formula is written. pfaffermayr2019gravity makes a similar point about the order of the bias of the standard errors for the two-way FE-PPML estimator, albeit using a slightly different analysis.}}
All together, this analysis implies that the estimated standard error for $\widehat{\beta}$ will exhibit a bias that only disappears at the relatively slow rate of $1/\sqrt{N}$. We should therefore be concerned that asymptotic confidence intervals for $\widehat{\beta}$ may exhibit inadequate coverage even in moderately large samples, similar to what has been found for the two-way FE-PPML estimator in recent simulation studies by egger_glm_2015, jochmans_two-way_2016, and pfaffermayr2019gravity. Indeed, the bias approximation we have derived in (ref) can be readily adapted to the two-way setting or even to more general settings with $k$-way fixed effects.
We now present two methods for correcting the bias in estimates: a jackknife method based on the split-panel jackknife of dhaene2015split and an analytical correction based on the expansion shown in Proposition (ref). We also provide an analytical correction for the downward bias in standard errors.
The advantage of the jackknife correction is that it does not require explicit estimation of the bias yet still has a simple and powerful applicability. To see this, note first that the asymptotic bias we characterize can be written as \[ \frac{1}{N}B^{\beta}+ o_p(N^{-1}), \label{eq:bias} \] where $B^{\beta}$ is a combined term that captures any suspected asymptotic bias contributions of order $1/N$. The specific jacknife we will apply for our current purposes is a split-panel jackknife based on dhaene2015split. As in dhaene2015split, we want to divide the overall data set into subpanels of roughly even size and then estimate $\widehat{\beta}_{(p)}$ for each subpanel $p$. Given the gravity structure of the model, we first divide the set of countries into evenly-sized groups $a$ and $b$. We then consider 4 subpanels of the form “$(a,b)$”, where “$(a,b)$” denotes a subpanel where exporters from group $a$ are matched with importers from group $b$. The other three subpanels are $(a,a)$, $(b,a)$, and $(b,b)$. For randomly-generated data, we can define $a$ and $b$ based on their ordering in the data (i.e., $a:={i:i\le N/2}$; $b:={i:i> N/2}$). For actual data, it would be more sensible to draw these subpanels randomly and repeatedly.\footnote{This is just one possible way to construct a jackknife correction for two-way panels. We have also experimented with splitting the panel one dimension at a time as in fernandez-val_individual_2016, but we find the present method performs significantly better at reducing the bias.}
The split-panel jackknife estimator for $\beta$, $\widetilde{\beta}_{N}^J$, is then defined as
This correction works to reduce the bias because, {so long as the distribution of $y_{ij}$ and $x_{ij}$ is homogeneous across both the $i$ and $j$ dimensions of the panel},\footnote{ { By “homogeneity” we mean that the vector $(y_{ij}, x_{ij}, \alpha_i, \gamma_j)$ is identically distributed across both $i$ and $j$, which is the appropriate translation of Assumption 4.3 in fernandez-val_individual_2016 to our setting. This does allow $(y_{ij}, x_{ij})$ to be heterogeneously distributed {conditional on the fixed effects} in the sense that the fixed effects themselves introduce heterogeneity into the model. Nonetheless, this is a strong assumption. One of the main advantages of the analytical bias correction is that it does not require such assumptions}. } each $\widehat{\beta}_{(p)}$ has a leading bias term equal to $2B^{\beta}/N$. The average $\widehat{\beta}_{(p)}$ across these four subpanels thus also has a leading bias of $2B^{\beta}/N$ and any terms depending on $B^{\beta}/N$ cancel out of (ref). Thus, the bias-corrected estimate $\widetilde{\beta}_{N}^{J}$ only has a bias of order $o_p(N^{-1})$, which is obtained by combining the second-order bias from $\widehat{\beta}$ with that of the average subpanel estimate. This latter bias can be shown to be larger than the original second-order bias in (ref), but the overall bias should still be smaller because of the elimination of the leading bias term.
{ Our analytical correction for the bias is based on the bias expression in Proposition (ref). In the Appendix, we show how appropriate sample analogs $\widehat{W}_{N}$, $\widehat{B}_{N}$, $\widehat{D}_{N}$ of the expressions for $W_N$, $B_N$, $D_N$ can be formed. The resulting bias-corrected estimate is then given by $$\widehat \beta - (N-1)^{-1}\widehat{W}_{N}^{-1}(\widehat{B}_{N}+\widehat{D}_{N}) . $$ It is possible to show that these plug-in corrections lead to estimates that are asymptotically unbiased as $N\rightarrow\infty$.} Still, for finite samples, it is evident that the bias in some of these plug-in objects could cause the analytical bias correction to itself exhibit some bias. For this reason, it is not obvious a priori whether the analytical correction will outperform the jackknife at reducing the bias in $\widehat{\beta}$. One clear advantage the analytical correction has over the jackknife is that {it does not require any homogeneity restrictions on the distribution of $y_{ij}$ and $x_{ij}$} in order to be valid.
Under the assumption of clustered errors within pairs, a natural correction for the variance estimate is available based on (ref). Specifically, let \[ \widehat{\Omega}^{U}\!:=\!\frac{1}{N\!\left(N\!-\!1\right)}\!\sum_{i,j}\widehat{\widetilde{x}}_{ij}\!\!\left[\mathbf{I}_{T}-\frac{1}{N\!\left(N\!-\!1\right)}\bar{H}_{ij}\widehat{\widetilde{x}}_{ij}\widehat{W}_{N}^{-1}\widehat{\widetilde{x}}^{\prime}\!-\frac{1}{N\!\left(N\!-\!1\right)}\bar{H}_{ij}d_{ij}\widehat{W}_{N}^{(\phi)-1}d_{ij}^{\prime}\right]^{-1}\!\!\!\widehat{S}_{ij}\widehat{S}_{ij}^{\prime}\widehat{\widetilde{x}}_{ij}, \] where $\mathbf{I}_{T}$ is a $T\times T$ identity matrix and $\widehat{W}_{N}^{(\phi)}$ is a plugin estimate for $W_{N}^{(\phi)}$. The corrected variance estimate is then given by
The logic of this adjusted variance estimate follows directly from Kauermann and Carroll (2001): if the PPML estimator is correctly specified (such that $E(S_{ij}S_{ij}^{\prime})=\kappa\bar{H}_{ij}$), then $\widehat{V}^{U}$ can be shown to eliminate the first-order bias in $\widehat{V}(\widehat{\beta}-\beta^{0})$ shown in (ref). It is not generally unbiased otherwise, but it is plausible that it should eliminate a significant portion of any downward bias under other variance assumptions as well.
{
As we have noted, implementations of our analytical corrections are available via our Stata command {\tt ppml_fe_bias}. Since this command can be applied to data sets that do not conform exactly to our theoretical framework, we provide here some brief comments on its applicability in general cases that may not be explicitly covered in the notation and formulas used above. For example, it is worth pointing out that our corrections may be used with data sets that have missing values. They also continue to apply if the data includes the diagonal ($i=j$) terms versus treating them as missing or inapplicable. Similarly, we can allow for the data to have unequal numbers of exporters and importers. In the latter case, we need the numbers of exporters and importers to grow at the same rate asymptotically for our results to apply.
One practicality that is not covered in our current framework is the case of “four way” gravity models that have an added index for industry. Though a full analytical characterization of this type of model is left for future work, our Appendix describes how a modified version of the heuristic from ARE may be used to assess the order of the bias and derive a suitable jackknife correction. Interestingly, this discussion reveals that the asymptotic bias problem may be more severe for four-way models than for three-way models. The heuristic we propose may be used to assess IPPs in other, non-trade settings as well.}
For our simulation analysis, we assume the following: (i) the data generating process (DGP) for the dependent variable is of the form $y_{ijt}=\lambda_{ijt}\omega_{ijt}$, where $\omega_{ijt}$ is a log-normal disturbance with mean $1$ and variance $\sigma_{ijt}^{2}$. (ii) $\beta=1$. (iii) The model-relevant fixed effects $\alpha$, $\gamma$, and $\eta$ are each $\sim{\cal N}(0,1/16)$. (iv) $x_{ijt}=x_{ijt-1}/2+\alpha+\gamma+\nu_{ijt}$, where $\nu_{ijt}\sim{\cal N}(0,1/16)$.\footnote{These assumptions on $\alpha$, $\gamma$, $\eta$, $x_{ijt}$, and $\nu_{ijt}$ are taken from fernandez-val_individual_2016. Notice that $x_{ijt}$ is strictly exogenous with respect to $\omega_{ijt}$ conditional on $\alpha$, $\gamma$, and $\eta$.} (v) Taking our cue from santos_silva_log_2006, we consider 4 different assumptions about the {disturbance} $\omega_{ijt}$:
where we also allow for serial correlation within pairs by imposing
such that the degree of correlation weakens for observations further apart in time.\footnote{The $0.3$ that appears here serves as a quasi-correlation parameter. Replacing $0.3$ with $1$ would be analogous to assuming disturbances are perfectly correlated within pairs. Replacing it with $0$ removes any serial correlation. Choosing other values for this parameter produces similar results.}
The relevance of these various assumptions to commonly used error distributions is best described by considering the conditional variance ${\rm Var}(y_{ijt}|x_{it},\alpha,\gamma,\eta)$. For example, DGP I assumes that the conditional variance is constant, as in a Gaussian process with i.i.d disturbances. In DGP II, the conditional variance equals the conditional mean, as in a Poisson distribution. DGP III\textemdash which we will also refer to as “log-homoskedastic”\textemdash is the unique case highlighted in santos_silva_log_2006 where the assumption that the conditional variance is proportional to the square of the conditional mean leads to a homoskedastic error when the model is estimated in logs using a linear model. Finally, DGP IV provides a “quadratic” error distribution that mixes DGP II and DGP III and also allows for overdispersion that depends on $x_{ijt}$. As shown by santos_silva_log_2006, this type of DGP tends to induce a relatively large degree of bias.
Tables (ref) and (ref) present simulation evidence comparing the uncorrected three-way FE-PPML estimator with results computed using the analytical and jackknife corrections described in Section (ref). As in the prior simulations, we again compute results for a variety of different panel sizes\textemdash in this case for $N=20,50,100$ and $T=2,5,10$.\footnote{Note that the trade literature currently recommends using wide intervals of 4-5 years between time periods so as to allow trade flows time to adjust to changes in trade costs (see cheng_controlling_2005.) Thus, for practical purposes, $T=10$ may be thought of as a relatively “long” panel in this context that might span 40+ years. IPPs do not necessarily vanish for larger values of $T$, as discussed further below.} In order to validate our analytical predictions regarding these estimates, we compute the average bias of each estimator, the ratios of the average bias to the average standard error and of the average standard error to the standard deviation of the simulated estimates, and the probability that the estimated $95\%$ confidence interval covers the true estimate of $\beta=1$. In particular, we expect that the bias in $\widehat{\beta}$ should be decreasing in either $N$ or $T$ but should remain large relative to the estimated standard error and induce inadequate coverage for small $T$. We are also interested in whether the usual cluster-robust standard errors accurately reflect the true dispersion of estimates. Results for DGPs I and II are shown in Table (ref), whereas Table (ref) shows results for DGPs III and IV.
The results in both tables collectively confirm the presence of bias and the viability of the analytical and jackknife bias corrections. The average bias is generally larger for DGPs I and IV than II and III. As expected, it generally falls with both $N$ and $T$ across all the different DGPs, though only weakly so for DGP III (the log-homoskedastic case), which generally only has a small bias.\footnote{Numerically, we have found that the two terms that appear in both $B_N$ and $D_N$ in Proposition (ref) tend to have opposite signs when the DGP is log-homoskedastic, thereby mitigating one another.} To use DGP II\textemdash the Poisson case, where PPML should otherwise be an optimal estimator\textemdash as a representative example, we see that the average bias falls from $3.840\%$ for the smallest sample where $N=20$, $T=2$ to a low of $0.249\%$ at the other extreme where $N=100,$ $T=10$. For DGP IV, the least favorable of these cases, the average bias ranges from $-6.1364\%$ down to $-1.878\%$. On the whole, these results support our main theoretical findings that $\beta$ should be consistently estimated even for {fixed} $T$ but has an asymptotic bias that depends on the number of countries and on the number of time periods.
Interestingly, while the average bias almost always decreases with $T$, the ratio of the bias to standard error usually does not, seemingly contrary to the expectations laid out in Remark (ref). {Evidently, increasing $T$ does not automatically reduce the bias at a rate of $1/T$. As we discuss in more detail below, researchers should thus be careful to note that the implications of Remark (ref) do not necessarily apply to settings with small $T$ or even moderately large $T$.} Furthermore, the estimated cluster-robust standard errors themselves clearly exhibit a bias in all cases as well. Even when $N=100$, SE/SD ratios are uniformly below 1; generally they are closer to $0.9$ or $0.95$, and for DGPs I and IV, they are often closer to $0.85$ or even $0.8$. Because of these biases, the simulated FE-PPML coverage ratios are unsurprisingly below the $0.95$ we would expect for an unbiased estimator.
Bias corrections to the point estimates do help with addressing some, but not all, of these issues. The jackknife generally performs more reliably than the analytical correction at reducing the average bias when compared across all values of $N$ and $T$; notice how, for the Poisson case, for example, the average bias left by the jackknife correction is never greater than $0.25\%$ in absolute magnitude, whereas the analytical-corrected estimates still have average biases ranging between $0.01\%$ and $1.12\%$. However, when $N=100$, the analytical correction begins to closely match the jackknife, especially when $T=10$. All the same, both corrections generally have a positive effect, and the better across-the-board bias-reduction performance of the jackknife comes at the important cost of a relatively large increase in the variance. Thus, the analytical correction generally performs as well as or better than the jackknife in terms of improving coverage even in the smaller samples. Neither correction is sufficient to bring coverage ratios to $0.95$, however, {though corrected Gaussian-DGP estimates and Poisson-DGP estimates both exceed $0.94$ using the analytical correction when $N=100$ and $T=10$, with the latter reaching $0.948$.}
Table (ref) then evaluates the efficacy of our bias correction for the estimated variance. Keeping in mind that this correction is calibrated for the case of a correctly specified variance (which corresponds to DGP II), we would naturally expect that the effect of this correction should vary depending on the conditional distribution of the data. In that light, it is encouraging that we observe positive effects across all cases. The best results by far are for the DGPs I, II, and III, where combining the analytical bias correction for the point estimates with the correction for the variance yields coverage ratios that range between $0.925$ and $0.952$ when $N$ is either 50 or 100 and generally get closer to the the target value of $0.95$ as either $N$ or $T$ increases. These corrections lead to dramatic improvements in coverage for DGP IV as well, but there the remaining biases in both the point estimate and the standard error remain large even for $N=100$ and $T=10$.
To summarize, these simulations suggest that combining an analytical bias correction for $\widehat{\beta}$ with a further correction for the variance based on (ref) should be a reliable way of reducing bias and improving coverage. At the same time, it should be noted that neither should be expected to offer a complete bias removal. For smaller samples, if reducing bias on average is heavily favored, {and if the distribution of $y_{ij}$ and $x_{ij}$ can be reasonably assumed to be homogeneous across $i$ and $j$}, then the split-panel jackknife method might be preferable to the analytical correction method.
{Based on our Remark (ref), one might expect that increasing the size of the time dimension should reduce the asymptotic bias in relation to the standard error. However, in the range of $T$ values we used in Tables (ref)-(ref), this is not what we observe. The question thus arises: can we say if there exists a “large enough” value of $T$ beyond which researchers may feel relatively secure about IPPs?
Figure (ref) addresses this question by presenting simulated bias/SE and coverage ratios for a wider range of $T$ values, spanning from $2$ to $100$. $N$ is fixed at 100, and the data is otherwise generated the same way as before. If we focus just on the first two DGPs, the Gaussian and Poisson cases, we do indeed observe steady improvements in both ratios as $T$ increases, though coverage fails to hit $0.95$ in either case. However, for both DGP III (log-homoskedastic) and DGP IV (quadratic), we actually observe bias/SE ratios getting worse as $T$ approaches 100. In the case of DGP IV, coverage actually gets worse as well.\footnote{ {For scale reasons, coverage results for DGP IV are not shown. For $T=100$, we find that coverage is only $0.52$ in this case. As shown in the left-hand panel of Figure (ref), the reason is because the bias tends to decrease more slowly than the standard error as $T$ increases while $N$ is fixed under this DGP.}} The main takeaway is that gravity panels with seemingly large time spans are not necessarily immune to IPPs. To reconcile these findings with our theory, note that Remark (ref) only says that three-way PPML estimates become asymptotically unbiased as both $N$ and $T$ become large simultaneously. In further simulations, we have confirmed that both the bias/SE ratio and coverage improve across all DGPs when we compare, e.g., $N=T=200$ with $N=T=100.$}
For our {main} empirical application, we estimate the average effects of an FTA for a variety of different industries using a panel with a relatively large number of countries. The value of this exercise is that we expect that trade flows could be distributed very differently across different industries. Based on our results so far, this should lead to a range of differerent bias behaviors in the data.
Our trade data is from the BACI database of gaulier_baci:_2010, from which we extract data on trade flows between 167 countries for the years 1995, 2000, 2005, 2010, and 2015. Countries are chosen so that the same 167 countries always appear as both exporters and importers in every period; hence, the data readily maps to the setting just described with $N=167$ and $T=5$. We combine this trade data with data on FTAs from the NSF-Kellogg database maintained by Scott Baier and Jeff Bergstrand, which we crosscheck against data from the WTO in order to incorporate agreements from more recent years.\footnote{This database is available for download on Jeff Bergstrand's website: \url{https://www3.nd.edu/ jbergstr/}. The most recent version runs from 1950-2012. The additional data from the WTO is needed to capture agreements that entered into force between 2012 and 2015. } The specification we estimate is
where $y_{ijt}$ is trade flows (measured in current USD), $FTA_{ijt}$ is a $0/1$ dummy for whether or not $i$ and $j$ have an FTA at time $t$, and $\omega_{ijt}$ is an error term. As we have noted, estimation of specifications such as (ref) via PPML has become an increasingly standard method for estimating the effects of FTAs and other trade policies and is currently recommended as such by the WTO (see yotov_advanced_2016.)
Table (ref) presents results from FE-PPML estimation of (ref), including results obtained using our bias corrections. The estimation is applied separately to 28 2 digit ISIC (rev. 3) industries as well as to aggregate trade. While we do indeed see a range of different biases in the industry-level estimates, the results for aggregate trade flows, shown in the bottom row of Table (ref), are fairly representative. To provide some basic interpretation, the coefficient on $FTA_{ijt}$ for aggregate trade is initially estimated to be $0.082$, which equates to an $e^{0.082}-1=8.5\%$ average “partial” effect of an FTA on trade.\footnote{The term “partial effect” is conventionally used to distinguish this type of estimate from the “general equilibrium” effects of an FTA, which would typically be calculated by solving a general equilibrium trade model where prices, incomes, and output levels (which are otherwise absorbed by the $\alpha_{it}$ and $\gamma_{jt}$ fixed effects) are allowed to evolve endogenously in response to the FTA. In the context of such models, $\beta$ can usually be interpreted as capturing the average effect of an FTA on bilateral trade frictions specifically, holding fixed all other determinants of trade.} The estimated standard error is $0.027$, implying that this effect is statistically different from zero at the $p<0.01$ significance level. Our bias-corrected estimates do not paint an altogether different picture, but do highlight the potential for meaningful refinement. Both the analytical and jackknife bias corrections for $\beta$ suggest a downward bias of $0.04$-$0.06$, or about $15\%$-$22\%$ of the estimated standard error. As our bias-corrected standard errors show (in the last column of Table (ref)), the initially estimated standard error itself has an implied downward bias of $11\%$ (i.e., $0.027$ versus $0.030$).
Turning to the industry-level estimates, the analytical bias correction more often than not indicates a downward bias ranging between $5\%$-$20\%$ of the estimated standard error, though exceptions are present on both sides of this range. Estimates for the Chemical and Furniture industries appear to be unbiased, for example, and some (such as Tobacco) are associated with an upward bias. On the other end of the spectrum, implied downward biases can also be larger than $20\%$ of the standard error, as is seen for Petroleum ($47\%$), Fabricated Metal Products ($31\%$), Electrical Equipment ($27\%$), and Agriculture ($22\%$). The biases implied by the jackknife are often even larger (see Fabricated Metal Products, for example), consistent with what we found in our simulations for smaller panel sizes. One possible interpretation is that the jackknife-corrected estimates are giving us a less conservative alternative to the analytical corrections in these cases. Indeed, the general correspondence between the two sets of results adds validity to both methods. However, as we have noted, these jackknife estimates could also be reflecting {non-homogeneity in the data} and/or the higher variance introduced by the jackknife. Implied biases in the standard error, meanwhile, tend to range between $10\%$-$20\%$ of the original standard error, again with some exceptions.
{To further illustrate the types of results that can occur, we also obtain replication data for several recent articles that have used three-way gravity models and re-examine their findings using our bias corrections. The results, reported in Table (ref), help to demonstrate how these corrections can matter for assessing statistical significance. Instances where conventional significance levels are affected include the coefficients for $\text{EIA}\times\text{CONTIG}$ and $\text{EIA}\times\text{LEGAL}$ from baier2018heterogeneous and the coefficient for log approval rating from rose2019soft.\footnote{Note that baier2018heterogeneous and baier_economic_2014 use first-differenced OLS with added pair time trends. We estimate the three-way FE-PPML equivalents.} The implied bias to standard error ratios are sometimes as large as 40%-45%, as occurs for the $\text{EIA}\times\text{LANG}$ coefficient from baier2018heterogeneous and the Total EIA Effect from bergstrand_economic_2015. However, the largest effects are actually for the standard error, which is downward-biased by more than 40% in several cases (the Total FTA effect from baier2019widely, for example). Notably, there are meaningful differences found even for the data used in larch2019currency, a large data set with 213 countries and 66 time periods. These results reinforce our earlier finding that bias corrections may be useful even in settings with seemingly large $N$ and $T$}.
{
Thanks to recent methodological and computational advances, nonlinear models with three-way fixed effects have become increasingly popular for investigating the effects of trade policies on trade flows. However, the asymptotic and finite-sample properties of three-way fixed effects estimators have not been rigorously studied, especially with regards to potential IPPs. The performance of the FE-PPML estimator in particular is of natural interest in this context, both because FE-PPML is known to be relatively robust to IPPs as well as because it is likely to be a researcher's first choice for estimating three-way gravity models. Our results regarding the consistency of PPML in this setting reflect these unique properties of PPML and support its current status as a workhorse estimator for estimating the effects of trade polices.
Given the consistency of PPML in this setting, and given the nice IPP-robustness properties of PPML in general, it may come as a surprise that three-way PPML estimates nonetheless suffer from an asymptotic bias that affects the validity of inferences. In theory, the bias should become less of a problem when the country and time dimensions are both large, but our experiments with the time dimension indicate the bias can be of comparable magnitude to the standard error even in ostensibly large trade data sets. Typical cluster-robust estimates of the standard error are also biased, implying estimated confidence intervals not only off-center but also too narrow.
These issues are not so severe that they leave researchers in the wilderness, but we do recommend taking advantage of the corrective measures we have described. In particular, we find that analytical bias corrections based on Taylor expansions to both the point estimates and standard errors generally lead to improved inferences when applied simultaneously. We caution that we have not found these corrections to be a panacea, however, and several avenues remain open for future work. For example, confidence interval estimates could be adjusted further to account for the uncertainty in the estimated variance\textemdash kauermann2001note describe such a correction for the PPML case. A quasi-differencing approach similar to jochmans_two-way_2016 could provide another angle of attack, and a recent contribution by pfaffermayr2021confidence suggests that jackknife and bootstrap confidence interval methods hold promise as well. Turning to broader applications, the essential dyadic structure of our bias corrections could be easily {adapted} to network models that study changes in network behavior over time, including settings that involve studying the number of interactions between network members. }
{.67em}