Approximate functional differencing
Abstract
EconStor is a publication server for scholarly economic literature, provided as a non-commercial public service by the ZBW.
Full text
Dhaene, Geert; Weidner, Martin Article Approximate functional differencing SERIEs - Journal of the Spanish Economic Association Provided in Cooperation with: Spanish Economic Association Suggested Citation: Dhaene, Geert; Weidner, Martin (2023) : Approximate functional differencing, SERIEs - Journal of the Spanish Economic Association, ISSN 1869-4195, Springer, Heidelberg, Vol. 14, Iss. 3/4, pp. 379-416, https://doi.org/10.1007/s13209-023-00283-1 This Version is available at: https://hdl.handle.net/10419/286582 Standard-Nutzungsbedingungen: Die Dokumente auf EconStor dürfen zu eigenen wissenschaftlichen Zwecken und zum Privatgebrauch gespeichert und kopiert werden. Sie dürfen die Dokumente nicht für öffentliche oder kommerzielle Zwecke vervielfältigen, öffentlich ausstellen, öffentlich zugänglich machen, vertreiben oder anderweitig nutzen. Sofern die Verfasser die Dokumente unter Open-Content-Lizenzen (insbesondere CC-Lizenzen) zur Verfügung gestellt haben sollten, gelten abweichend von diesen Nutzungsbedingungen die in der dort genannten Lizenz gewährten Nutzungsrechte. Terms of use: Documents in EconStor may be saved and copied for your personal and scholarly purposes. You are not to copy documents for public or commercial purposes, to exhibit the documents publicly, to make them publicly available on the internet, or to distribute or otherwise use the documents in public. If the documents have been made available under an Open Content Licence (especially Creative Commons Licences), you may exercise further usage rights as specified in the indicated licence. https://creativecommons.org/licenses/by/4.0/
SERIEs (2023) 14:379–416 https://doi.org/10.1007/s13209-023-00283-1 ORIGINAL ARTICLE Approximate functional differencing Geert Dhaene1·Martin Weidner2 Received: 1 February 2023 / Accepted: 12 May 2023 / Published online: 9 June 2023 © The Author(s) 2023 Abstract Inference on common parameters in panel data models with individual-specific fixed effects is a classic example of Neyman and Scott’s (Econometrica 36:1–32, 1948) incidental parameter problem (IPP). One solution to this IPP is functional differencing (Bonhomme in Econometrica 80(4):1337–1385, 2012), which works when the number of time periods Tis fixed (and may be small), but this solution is not applicable to all panel data models of interest. Another solution, which applies to a larger class of models, is “large-T” bias correction [pioneered by Hahn and Kuersteiner (Econometrica 70(4):1639–1657, 2002) and Hahn and Newey (Econometrica 72(4):1295–1319, 2004)], but this is only guaranteed to work well when Tis sufficiently large. This paper provides a unified approach that connects these two seemingly disparate solutions to the IPP. In doing so, we provide an approximate version of functional differencing, that is, an approximate solution to the IPP that is applicable to a large class of panel data models even when Tis relatively small. Keywords Panel data ·Discrete choice ·Incidental parameters ·Bias correction · Functional differencing JEL Classification C23 We thank Stéphane Bonhomme for useful comments and discussions, and a referee for useful comments. This research was supported by the European Research Council Grant ERC-2018-CoG-819086-PANEDA and the Flemish Research Council Grant G073620N. BMartin Weidner [email protected] Geert Dhaene [email protected] 1KU Leuven, Leuven, Belgium 2University of Oxford, Oxford, UK 123
380 SERIEs (2023) 14:379–416 1 Introduction Panel data offer the potential to account for unobserved heterogeneity, typically through the inclusion of unit-specific parameters; see Arellano (2003) and Arellano and Bonhomme (2011) for reviews. Nonlinear panel data models, however, remain challenging to estimate, precisely because in many models the presence of unit-specific—or “incidental”—parameters makes the maximum likelihood estimator (MLE) of the common parameters inconsistent when the number of observations per unit, T, is finite (Neyman and Scott 1948). The failure of maximum likelihood has prompted two kinds of reactions. One approach is to look for point-identifying moment conditions that are free of incidental parameters. Such moment conditions can come from a conditional or a marginal likelihood (e.g., Rasch 1960; Lancaster 2000), an invariant likelihood (Moreira 2009), an integrated likelihood (Lancaster 2002), functional differencing (Bonhomme 2012), or from some other reasoning to eliminate the incidental parameters, for example, differencing in linear dynamic models (e.g., Arellano and Bond 1991). This approach is usually model-specific and is “fixed-T”, i.e., it seeks consistent estimation when T is fixed (and usually small). However, point-identifying moment conditions for small Tmay not exist because point identification simply may fail; see Honoré and Tamer (2006) and Chamberlain (2010) for examples. The other main approach is motivated by “large-T” arguments and seeks to reduce the large-Tbias of the MLE or of the likelihood function itself or its score function (e.g., Hahn and Kuersteiner 2002; Alvarez and Arellano 2003; Hahn and Newey 2004; Arellano and Bonhomme 2009; Bonhomme and Manresa 2015; Dhaene and Jochmans 2015b; Arellano and Hahn 2016; Fernández-Val and Weidner 2016). This approach is less model-specific and may also be applied to models where point identification fails for small T. The functional differencing method of Bonhomme (2012) provides an algebraic approach to systematically find valid moment conditions in panel models with incidental parameters—if such moment conditions exist. Related ideas are used in Honoré (1992), Hu (2002), Johnson (2004), Kitazawa (2013), Honoré and Weidner (2020), Honoré, Muris and Weidner (2021), and Davezies, D’Haultfoeuille and Mugnier (2022). In this paper, we extend the scope of functional differencing to models where point identification may fail. In such models, exact functional differencing (as in Bonhomme 2012) is not possible, but an approximate version thereof yields moment conditions that are free of incidental parameters and that are approximately valid in the sense that their solution yields a point close to the true common parameter value. Bonhomme’s method relies on the existence of (one or more) zero eigenvalues of a matrix of posterior predictive probabilities (or a posterior predictive density function) defined by the model. Our extension considers the case where all eigenvalues are positive and, therefore, point identification fails, but where some eigenvalues are very close to zero. This occurs as the number of support points of the outcome variable increases. Eigenvalues close to zero then lead to approximate moment conditions obtained as a bias correction of an initially chosen moment condition. The bias correction can be iterated, possibly infinitely many times. In point-identified models, the infinitely iterated bias correction is equivalent to functional differencing. Therefore, approximate functional 123
SERIEs (2023) 14:379–416 381 differencing can be viewed as finite-Tinference in point-identified models, and as a large-Titerative bias correction method in models that are not point-identified. The construction of approximate moment conditions is our main focus. Once such moment conditions are found, estimation follows easily using the (generalized) method of moments, and the discussion of estimation is therefore deferred to later parts of the paper (from Sect.6). We illustrate approximate functional differencing in a probit binary choice model. The implementation, including the iteration, is straightforward, only requiring elementary matrix operations. Indeed, one of the contributions of this paper is to show how to iterate score-based bias correction methods for discrete choice panel data relatively efficiently. After introducing the model setup in Sect. 2, we review the main ideas behind functional differencing in Sect.3. In Sect.4we introduce our novel bias corrections and explain how they relate to functional differencing. Section5examines the eigenvalues of the matrix of posterior predictive probabilities in a numerical example. Section6 briefly discusses estimation. Further numerical illustration of the methods and some Monte Carlo simulation results are presented in Sect.7. Section8discusses some extensions, in particular, a generalization of the estimation method to average effects. Finally, we provide some concluding remarks in Sect. 9. 2 Setup We observe outcomes Yi∈Yand covariates Xi∈Xfor units i=1,...,n.We only consider finite outcome sets Yin this paper, but in principle all our results can be generalized to infinite sets Y. There are also latent variables Ai∈A, which are treated as nuisance parameters. We assume that (Yi,Xi,Ai),i=1,...,n, are independent and identical draws from a distribution with conditional outcome probabilities Pr Yi=yiXi=xi,Ai=αi=fyixi,α i,θ 0,(1) where the function fyixi,α i,θis known (this function specifies “the model”), but the true parameter value θ0∈⊂Rdθis unknown. Our primary goal in this paper is inference on θ0. Let π0(αi|xi)be the true distribution of Aiconditional on Xi=xi. Then, the conditional distribution that can be identified from the data is Pr Yi=yiXi=xi=A fyixi,α i,θ 0π0(αi|xi)dαi.(2) No restrictions are imposed on π0(αi|xi)nor on the marginal distribution of Xi, that is, we have a semi-parametric model with unknown parametric component θ0and unknown nonparametric component π0(αi|xi). The setup just described covers many nonlinear panel data models with fixed effects. There, we observe outcomes Yit and covariates Xit for unit iover time periods t=1,...,T. For static panel models we then set Yi=(Yi1,...,YiT)and Xi=(Xi1,...,XiT), and the model is typically specified as 123
382 SERIEs (2023) 14:379–416 fyixi,α i,θ= T t=1 f∗yit xit,α i,θ, where f∗yit xit,α i,θ 0=Pr Yit =yit Xit =xit,Ai=αi.Here, f∗yit xit, αi,θ)often depends on xit,αi,θonly through a single index x itθ+αi, where θis a regression coefficient vector of the same dimension as xit, and αi∈Ris an individualspecific fixed effect. Of course, θmay also contain additional parameters (e.g., the variance of the error term in a Tobit model). For dynamic nonlinear panel models, we usually have to model the dynamics explicitly. For example, we may include a lagged dependent variable in the model. In that case, assuming that Yit at t=0 is observed, we have Yi=(Yi1,...,YiT)and Xi=(Yi0,Xi1,...,XiT), and the model is usually specified as fyixi,α i,θ= T t=1 f∗yit yi,t−1,xit,α i,θ, where f∗yityi,t−1,xit,α i,θ=Pr Yit =yitYi,t−1=yi,t−1,Xit =xit,Ai=αi. Here, the initial observation Yi0is included in the conditioning variable Xi.Inthis way, the setup in Eqs. (1) and (2) also covers dynamic panel data models. The setup may also be relevant for applications outside of standard panel data, e.g., pseudo-panels, network models, or games. But one typically needs Yito be a vector of more than one outcome to learn anything about θ0since, in most models, the value of αialone can fully fit any possible outcome value if there is only a single outcome per unit (i.e., if the sample is purely cross-sectional). The main insights of our paper are therefore applicable more broadly, but our focus will be on panel data. In particular, the following static binary choice panel data model will be our running example throughout the paper. Example 1A (Static binary choice panel data model) Consider a static panel data model with Yi=(Yi1,...,YiT)and Xi=(Xi1,...,XiT)where the outcomes Yit ∈{0,1} are generated by Yit =1(X it θ0+Ai≥Uit) and the errors Uit are independent of Xiand Ai∈R, and are i.i.d. across iand twith cdf F(u). This implies fyixi,α i,θ= T t=1[1−F(x it θ+αi)]1−yit [F(x it θ+αi)]yit . For the probit model, we have F(u)=(u), where is the standard normal cdf, and for the logistic model we have F(u)=(1+e−u)−1. To make the example even more specific, we consider a single binary covariate Xit ∈{0,1}such that for all i=1,...,nwe have 123
SERIEs (2023) 14:379–416 383 Xit =1(t>T0), for some T0∈{1,...,T−1}, that is, Xit is equal to zero for the initial T0time periods, and is equal to one for the remaining T1=T−T0time periods. Here T0, and therefore Xi, is non-random and constant across i, so we can simply write fyiαi,θinstead of fyixi,α i,θ. The parameter of interest, θ0∈R, is one-dimensional. Example 1B (Example 1A reframed) Consider Example 1A, but denote the binary outcomes now as Y∗ it ∈{0,1}and define the outcome Yifor unit ias the pair Yi=Yi,0,Yi,1:= ⎛ ⎝ T0 t=1 Y∗ it, T t=T0+1 Y∗ it⎞ ⎠∈{0,...,T0}×{0,...,T1}=Y.(3) Here, Yi,0=T t=1Y∗ it (1−Xit)is the number of outcomes for unit ifor which Y∗ it =1 within those time periods that have Xit =0, while Yi,1=T t=1Y∗ it Xit is the number of outcomes with Y∗ it =1 for the time periods with Xit =1. This implies that fyiαi,θ=T0 yi,0[1−F(αi)]T0−yi,0[F(αi)]yi,0T1 yi,1 ×[1−F(θ +αi)]T1−yi,1[F(θ +αi)]yi,1, wherewedropxifrom fyixi,α i,θsince it is non-random and constant across i. The parameter of interest, θ0∈R, is unchanged. From the perspective of parameter estimation, Example 1B is completely equivalent to Example 1A, because Yiin Example 1B is a minimal sufficient statistic for the parameters (θ0,α i)in Example 1A. Nevertheless, the outcome space in Example 1A is larger (|Y|=2T) than the outcome space in Example 1B (|Y|=(T0+1)(T1+1)), and this will make a difference in our discussion of moment conditions in these two examples below. 3 Main idea behind functional differencing We now explain the main idea behind the functional differencing method of Bonhomme (2012). Our presentation is similar to that in Honoré and Weidner (2020). However, our goal here is much closer to that in Bonhomme’s original paper because we want to describe a general estimation method, one that is applicable to a very large class of models, as opposed to obtaining an analytical expression for moment conditions in specific models. 3.1 Exact moment conditions Consider the model described by (1) and (2), where our goal is to estimate θ0. Functional differencing (Bonhomme 2012) aims to find moment functions m(yi,xi,θ) ∈ 123
384 SERIEs (2023) 14:379–416 Rdmsuch that the model implies, for all xiand αi, that Em(Yi,Xi,θ 0)Xi=xi,Ai=αi=0(4) or, equivalently, y∈Y m(y,xi,θ) fyxi,α i,θ=0, since we want (4) to hold for all possible θ0∈. Verifying that m(yi,xi,θ) satisfies this conditional moment condition only requires knowledge of the model fyixi,α i,θ, not of the observed data. Note that m(yi,xi,θ)does not depend on αi, but nevertheless should have zero mean conditional on any realization Ai=αi. This is a strong requirement, and we will get back to this below. Once we have found such valid moment functions m(yi,xi,θ), we can choose an arbitrary (matrix-valued) function g(xi,θ)∈Rdm×dm, and define m(yi,xi,θ):= g(xi,θ)m(yi,xi,θ), which is a vector of dimension dm. By the law of iterated expectations, we then obtain, under weak regularity conditions, the unconditional moment condition E[m(Yi,Xi,θ 0)]=0,(5) which we can use to estimate θ0by the generalized method of moments (GMM, Hansen 1982). The nuisance parameters αido not feature in the GMM estimation at all, that is, functional differencing provides a solution to the incidental parameter problem (Neyman and Scott 1948). Of course, the key condition for consistent GMM estimation is that E[m(Yi,Xi,θ) ]= 0 for any θ= θ0. This identification condition is violated if m(yi,xi,θ) does not depend on θ(a special case of which is m(yi,xi,θ)=0, which is a trivial solution to (4). Hence the moment functions must depend on θto be informative about θ0. Uninformative moment functions in Example 1A To give an example of a moment function that is uninformative about θ0, consider Example 1A.Lettand sbe two time periods where Xit =Xis.LetYi,−(t,s)∈ {0,1}T−2be the outcome vector Yifrom which the outcomes Yit and Yis are dropped. Then, since Xit =Xis, the outcomes Yit and Yis are exchangeable and therefore EYit Yi,−(t,s)=EYis Yi,−(t,s). This implies that for any function g:{0,1}T−2→Rthe moment function m(yi,xi,θ):= (yit −yis)g(yi,−(t,s))(6) 123
SERIEs (2023) 14:379–416 385 satisfies (4). This moment function does not depend on θand is therefore not useful for parameter estimation. (It is useful for model specification testing, but we will not discuss this.) Furthermore, one can show that every moment function m(yi,xi,θ)that satisfies (4) in Example 1A is equal to a corresponding valid moment function in Example 1B plus a linear combination of moment functions of the form (6).1Thus, from the perspective of constructing valid moment functions that are informative about θ0, without loss of generality we can focus on Example 1B instead of Example 1A. Example 1A is useful because it is a completely standard panel model and it gives a simple example of valid moment functions that do not depend on θ. From here onward, however, we will always use Example 1B as our running example. Informative moment functions in Example 1B for logistic errors Consider Example 1B with logistic error distribution, F(u)=(1+e−u)−1. Then, Yi,0+Yi,1is a sufficient statistic for Ai, so the distribution of Yiconditional on Yi,0+ Yi,1does not depend on Ai. It is well known that this implies that the corresponding conditional MLE of θ0is consistent as n→∞, for any fixed T≥2; see, e.g., Rasch (1960), Andersen (1970), and Chamberlain (1980). Here, instead of considering conditional maximum likelihood, we focus purely on the existence of moment conditions. Let ¯y=(¯y0,¯y1)∈{0,...,T0}×{0,...,T1} and ˜y=(˜y0,˜y1)∈{0,...,T0}×{0,...,T1}be two possible realizations of Yisuch that ¯y0+¯y1=˜y0+˜y1and ¯y=˜y. Since Yi,0+Yi,1is a sufficient statistic for Ai,it must be the case that the ratio r(θ) := f¯yαi,θ f˜yαi,θ does not depend on αi. This implies that m(yi,θ):= 1{yi=¯y}−r(θ) 1{yi=˜y}(7) satisfies Em(Yi,θ 0)Ai=αi=0. A short calculation gives r(θ) =T0 ¯y0T1 ¯y1T0 ˜y0−1T1 ˜y1−1 exp[(¯y1−˜y1)θ]. Since we assume ¯y=˜y, the moment function m(yi,θ)indeed depends on θ. Furthermore, m(yi,θ) is strictly monotone in θwhen yi=˜yand constant in θotherwise, 1Let yi=y(y∗ i)be the mapping between an outcome y∗ iin Example 1A and an outcome yiin Example 1B, as defined by (3), and let Y∗(yi)=y∗ i:yi=y(y∗ i)be the set of outcomes y∗ ithat map to yi. Starting from a valid moment function mA(y∗ i,θ)in Example 1A we obtain a valid moment function in Example 1B as mB(yi,θ)=Y∗(yi)−1y∗ i∈Y∗(yi)mA(y∗ i,xi,θ). The null space of this linear mapping mA→ mB is spanned by moment functions of the form (6). This implies that the difference mA(y∗ i,θ)−mB(y(y∗ i), θ) is a linear combination of moment functions of the form (6). 123
386 SERIEs (2023) 14:379–416 and all outcomes are realized with positive probability. Hence E[m(Yi,θ) ]is strictly monotone in θand the condition E[m(Yi,θ 0)]=0 uniquely identifies θ0. The observation that the existence of a sufficient statistic for the nuisance parameter Aiallows for identification and estimation of θ0is quite old (e.g., Rasch 1960). However, the reason that the functional differencing method is truly powerful is that moment functions satisfying (4) and identifying θ0may exist even in models where no sufficient statistic for Aiis available. Examples of this are given by Honoré (1992), Hu (2002), Johnson (2004), Kitazawa (2013), Honoré and Weidner (2020), Honoré, Muris and Weidner (2021), and Davezies, D’Haultfoeuille and Mugnier (2022). Bonhomme (2012) provides a computational method for obtaining moment functions m(yi,xi,θ) such that (4) holds in a large class of models, while Honoré and Weidner (2020) discuss how to obtain explicit algebraic formulas for moment conditions in specific models. Dobronyi, Gu and Kim (2021) show that additional moment inequalities may exist that contain identifying information on θ0that is not contained in the moment equalities. Our example of a moment function in (7) is convenient and easy to understand, but it is not really representative of the potential complexity of more general moment functions. The papers cited in the previous paragraph give a better view of the true capability of the functional differencing method in more challenging settings. 3.2 Approximate moment conditions Functional differencing is a very powerful and useful method. Nevertheless, there are many models to which it is not applicable. The reason is that the condition in Eq. (4)is actually quite strong. It requires us to find a function m(Yi,Xi,θ)that does not depend on Aiat all, but that is supposed to have a conditional mean of zero for any possible realization of Ai. In most standard panel data models Aitakes values in R(Aican also be a vector), implying that (4) imposes an infinite number of linear restrictions. It is therefore perhaps unsurprising that there are many panel data models for which (4) has no non-trivial solution at all. In Example 1B we have shown the existence of valid moment functions for the logit model, but it turns out that no valid moment function exists for the probit model when θ0= 0 (we have verified this non-existence numerically for many values of Tand T0). Instead of trying to find moment functions satisfying (4), and hence (5), exactly, we argue that it can also be fruitful to search for moment functions that satisfy these conditions only approximately, i.e., E[m(Yi,Xi,θ 0)]≈0.(8) For a given model fyixi,α i,θwe might not be able to find an exact solution to (5), but we might be able to find a very good approximate solution. Examples of approximate moment conditions are provided by the “large-T” panel data literature, which considers asymptotic sequences where also T→∞(jointly with n→∞). To illustrate the insights of this literature, let αi(θ) be the MLE of αi obtained from maximizing fYiXi,α i,θover αi∈A, and let ψ(yi,xi,α i,θ) be a moment function that satisfies E[ψ(Yi,Xi,Ai,θ 0)]=0 in model (1), e.g., 123
SERIEs (2023) 14:379–416 393 Remark 2 If the set Ais finite with cardinality nA=|A|, then by construction rank[Q(x,θ)]≤nA. Thus, whenever nA<nY,Q(x,θ)has nY−nAzero eigenvalues, implying that exact moment functions, free of α, are available. Notice, however, that this assumes not only that αtakes on only a finite number of values, but also that these values are known (they constitute the known set A). By contrast, the literature on discretizing heterogeneity in panel data (e.g., Bonhomme and Manresa 2015; Su, Shi and Phillips 2016; Bonhomme, Lamadon and Manresa 2022) usually considers the support points of Ato be unknown. For our purposes, the fact that rank[Q(x,θ)]≤nA matters only in our numerical implementation, where the rank of Q(x,θ) might be truncated my the discretization of the set A. 5 Eigenvalues of Q(x,Â): numerical example Lemma 1guarantees that all eigenvalues of the matrix Q(x,θ)lie in the interval [0,1], and Lemma 2shows that exact moment conditions that are free of the incidental parameter Aare only available if Q(x,θ) has a zero eigenvalue. However, even in models where Q(x,θ) does not have a zero eigenvalue, we suggest that calculating the eigenvalues of Q(x,θ)is generally informative about whether moment conditions exist that are approximately free of the incidental parameters. This is because in typical applications we expect that the distinction between a zero eigenvalue and a very small eigenvalue of Q(x,θ) should be practically irrelevant, that is, as long as Q(x,θ) has one or more eigenvalues that are very close to zero, then very good approximate moment conditions in the sense of (8) should exist. It is difficult to make a general statement about how small an eigenvalue of Q(x,θ) needs to be to qualify as sufficiently small. However, in a typical model with a sufficiently large number nYof outcomes (which for discrete choice panel data usually requires only moderately large T) one will often have multiple eigenvalues of Q(x,θ) that are so small (say smaller than 10−5) that there is little doubt that they can be considered equal to zero for practical purposes. To illustrate this, consider Example 1B with normally distributed errors, F(u)= (u), even values of T, and T0=T1=T/2, which implies nY=(1+T/2)2.Forthe prior distribution of Awe choose the standard normal distribution, πprior(α) =φ(α).2 We then calculate the eigenvalues of the nY×nYmatrix Q(θ) for θ=1(thereare no longer covariates xin this example as they are assigned non-random values). For T=2wehavenY=4, and the four eigenvalues of Q(1)are λ1=1, λ2=0.47463, λ3=0.10727, and λ4=0.00016. For T=4 and T=6wehavenY=9 and nY=16, respectively, and the corresponding eigenvalues of Q(1)are plotted in Figs.1 and 2. Figure3plots only the smallest eigenvalues of Q(1)for T=2,4,...,20. From these figures, we see that for T≥4 the smallest eigenvalue of Q(1)is less than 10−9, which we argue can be considered equal to zero for practical purposes. In 2In fact, for our numerical implementation, we discretize the standard normal prior by choosing 1000 grid points αj=−1(j/1001),j=1,...,1000, and we implement a prior that gives equal probability to each of these grid points. The approximation bias that results from this discretization is negligible for our purposes, as long as nYis much smaller than 1000. 123
394 SERIEs (2023) 14:379–416 123456789 0 0.2 0.4 0.6 0.8 1 j eigenvalue λj 123456789 10−10 10−8 10−6 10−4 10−2 100 j eigenvalue λj Fig. 1 The eigenvalues λj(θ) of the matrix Q(θ) in Example 1B are plotted for the case θ=1, T=4, T0=T1=2, and where both the error distribution and the prior distribution of Aare standard normal. The left and right plots show the same eigenvalues, just with a different scaling of the y-axis 051015 0 0.2 0.4 0.6 0.8 1 j eigenvalue λj 051015 10−20 10−15 10−10 10−5 100 j eigenvalue λj Fig. 2 Same eigenvalue plot as in Fig.1,butforT=6andT0=T1=3 Fig. 3 For the same setting as in Fig.1, but for different values of T(with T0=T1=T/2), we plot only the smallest eigenvalue of Q(θ) for θ=1. Notice that the smallest eigenvalue is never zero, that is, Q(1)has full rank for all values of Tconsidered 2468101214161820 10−163 10−123 10−83 10−43 10−3 T smallest eigenvalue Figs.1and 2we see that the largest eigenvalue, λ1, is equal to one (because Q(1)is a stochastic matrix), but then the eigenvalues λjdecay exponentially fast as jincreases.3 If we were to replace the standard normal distribution of the errors by the standardized logistic cdf F(u)=(1+e−πu/√3)−1(normalized to have variance one), 3The eigenvalues of Q(1)presented in this section were obtained using Mathematica with a numerical precision of 1000 digits. 123
SERIEs (2023) 14:379–416 395 then the left-hand side (non-logarithmic) plots in Figs.1and 2would look almost identical, but there would be (T/2)2eigenvalues exactly equal to zero. These zero eigenvalues for the logit model are due to the existence of a sufficient statistic for A and, correspondingly, the existence of exact moment functions (associated with the left null-space of Q(θ)), as discussed in Sect.3.1 above. Given that the change from standardized logistic errors to standard normal errors is a relatively minor modification of the model, it is not surprising that we see many eigenvalues close to zero in Figs. 1 and 2. Figure 3shows that the smallest eigenvalue of Q(1)in this example also decays exponentially fast as Tincreases.4However, for none of the values of Tthat we considered here, did we find an exact zero eigenvalue for the static binary choice probit model. We conjecture that this is true for all T≥2,butwehavenoproof. 5 This example illustrates that eigenvalues of Q(x,θ) very close to zero but not exactly zero may exist in interesting models. When aiming to estimate the parameter θ0in a particular model of the form (2), our first recommendation is to calculate the eigenvalues of Q(x,θ)for some representative values of θand xto see if some of them are zero or close to zero. If some are equal to zero, then exact functional differencing (Bonhomme 2012) is applicable. If some are very close to zero, then approximate moment functions (as in (8)) are available. The eigenvalues of Q(x,θ) are useful to examine whether exact or approximate moment functions for θare available in a given model. However, as explained in Sect.3.1, the corresponding moment functions also have to depend on θto be useful for parameter estimation. For example, the matrix Q(1)in Example 1A has exactly the same non-zero eigenvalues as the matrix Q(1)in Example 1B that we just discussed, but in addition, it has a zero eigenvalue with multiplicity equal to 2T−(T0+1)(T1+1), corresponding to the uninformative moment functions in equation (6). As a diagnostic tool, it can also be useful to calculate the matrix Q(x)for the model f(y|x,θ i,α i), which has no common parameters and where both θiand αiare individual-specific fixed effects (this requires choosing a prior for θas well, which may have finite support to keep the computation simple). Every zero eigenvalue of that matrix Q(x) then corresponds to a moment function, for that value of x, that does not depend on θ (within the range of the chosen prior for θ). The existence of uninformative moment functions (6) in Example 1A, for example, can be detected in this way. 4Presumably this finding and the fast shrinkage of the identified sets of common parameters (Honoré and Tamer 2006) and average effects (Chernozhukov, Fernández-Val, Hahn and Newey 2013) are manifestations of the same phenomenon. 5One needs to be careful with such conclusions for all T. For example, we also experimented with another error distribution. If, in Example 1B with θ=1andT0=T1=T/2, one chooses the error distribution F(u)to be the Laplace distribution with mean zero and scale one, then numerically we found that for any choice of prior the matrix Q(1)does not have a zero eigenvalue for T=2andT=4, but it does for T=6. So it is not impossible that something similar could happen for the probit model for sufficiently large T, although we do not expect it. 123
396 SERIEs (2023) 14:379–416 6 Estimation Suppose we have chosen a prior distribution πprior, an initial score function s(y,x,θ), and an order of bias correction, q. This gives the bias-corrected score function s(q)(y,x,θ) as the moment function m(y,x,θ) for which the approximate moment condition (8) is assumed to hold. For simplicity, suppose that ds=dθ, so that the number of moment conditions equals the number of common parameters we want to estimate. We can then define a pseudo-true value θ∗∈as the solution of E[m(Yi,Xi,θ ∗)]=0.(21) The corresponding method of moments estimator θsatisfies 1 nn i=1m(Yi,Xi, θ) = 0.Under appropriate regularity conditions, including existence and uniqueness of θ∗, we then have, as n→∞, √n( θ−θ∗)d →N(0,V∗), with asymptotic variance given by V∗=[G ∗]−1Var [m(Yi,Xi,θ ∗)]G−1 ∗,G∗=E∇θm(Yi,Xi,θ ∗).(22) In Sect.7we will report the bias θ∗−θ0and the asymptotic variance V∗for different choices of moment functions in Example 1B. Note that reporting the bias θ∗−θ0of the parameter estimates is more informative than reporting the bias E[m(Yi,Xi,θ 0)] of the moment condition, in particular since the moment condition can be rescaled by an arbitrary factor. How should qbe chosen? If Q(x,θ)is singular, then a natural choice is q=∞, as described in Sect.4.3, because this delivers an exactly unbiased moment condition. If Q(x,θ) is nonsingular but some of its eigenvalues are small, then, recalling our discussion in Sect.5, our general recommendation is to choose relatively large values of q. The larger the chosen q, the more we rely on the smallest eigenvalues of Q(x,θ), because contributions to s(q)(y,x,θ)from larger eigenvalues of Q(x,θ) are downweighted more heavily as qincreases. If none of the eigenvalues Q(x,θ) is close to zero, then there are no moment conditions that hold approximately in the sense of (8). Yet, even then, setting q>0 is likely to improve on q=0, even though the remaining bias will still be non-negligible in general. Whatever the eigenvalues of Q(x,θ)are (and, indeed, whether Q(x,θ)is singular or not), qis a tuning parameter and a principled way to choose qwould be to optimize some criterion, for example, the (estimated) mean squared error of θ. We leave this for further study. In our numerical illustrations in Sect.7we just consider finite values of qup to q=1000, and q=∞. 7 Asymptotic and finite-sample properties In this section, we report on asymptotic and finite-sample properties of θfor different choices of moment functions in the model of Example 1B with standard normal errors (i.e., the panel probit model with a single, binary regressor and fixed effects) and a 123
SERIEs (2023) 14:379–416 397 variation thereof, the model of Example 1A with a continuous regressor. Throughout, we set θ0=1, we use a standard normal prior (i.e., πprior(α) =φ(α),asinFigs.1 and 2), we choose the integrated score (12) as the initial score function, and we vary q, the number of iterations of the bias correction procedure. We first present results on asymptotic and finite-sample biases and variances of θ for three cases where Tis relatively small: Example 1B with T0=T1=T/2 and T∈{4,6}(Case 1); Example 1B with T0=1,T1=T−1 and T∈{4,10}(Case 2); Example 1A with a continuous regressor Xit ∼N(0.5,0.25)and T∈{4,6} (Case 3). In all three cases we set the true distribution of Aequal to N(1,1)(i.e., π0(α) =φ(α −1)); note that this implies that πprior is rather different from π0. Then, in the setup of Example 1B with standard normal errors, we numerically explore Conjecture 1 by examining the asymptotic bias, θ∗−θ0,forTup to 512, q up to 3, and various choices of π0as detailed below. 7.1 Case 1: binary regressor, T0=T1=T/2 Table 1reports θ∗−θ0and V∗for the case where T0=T1=T/2 and T∈{4,6}.The uncorrected estimate of θ0(q=0)has a large positive bias, 0.5050 when T=4 and 0.4056 when T=6. Bias correction (q>0) reduces the bias considerably, though non-monotonically in q. The least bias is attained at q=∞, where the bias is very small: −0.52 ×10−4when T=4 and −0.25 ×10−10 when T=6. Our calculation of V∗shows that there is, overall, a bias-variance trade-off in this example. When T=4, V∗slightly decreases as we move from q=0toq=1, but for q≥1 we see that V∗increases in q; when T=6, V∗increases in qthroughout. Strikingly, V∗at q=∞is much larger than, say, at q=1000.6 Table 1also reports, for a cross-sectional sample size of n=1000, the approximate RMSE of θand the approximate coverage rate of the 95% confidence interval with bounds θ±(n V∗)1/2−1(0.975), where V∗is the empirical analog to V∗. The approximate RMSE is calculated as RMSE =(V∗/n+(θ∗−θ0)2)1/2and the approximate coverage rate as CI0.95 =Pr[|Z|≤−1(0.975)]where Z∼ N((θ∗−θ0)/(V∗/n)1/2,1). Of course, RMSE and CI0.95 follow mechanically from θ∗,V∗,nand heavily depend on the chosen n. For our choice of n=1000, when T=4, the RMSE is minimized at q=2, but this is a rather fortuitous consequence of the bias having a local minimum (in magnitude) at q=2. In practice, when Tis very small we would not recommend choosing qless than 10, say, because otherwise, the remaining bias is often non-negligible. From q=10 or 20 onward, the bias and confidence interval coverage rates are reasonably good. On the other hand, we would also not recommend choosing qto be very large (including q=∞), one reason being asymptotic variance inflation. We also conducted a real Monte Carlo simulation, under the exact same setup as described (and with n=1000 in particular). Table 2gives the results, based on 1000 Monte Carlo replications. The column “bias” is E( θ−θ0)(estimated by Monte Carlo), 6The limit limq→∞ θ∗(q)can be obtained by solving E[S(θ∗)UnY(θ∗)[U−1(θ∗)]nYδ(Y)]=0forθ∗, where UnY(θ) is the submatrix of U(θ) whose columns are the right-eigenvectors of Q(θ) corresponding to λnY(θ), the smallest eigenvalue of Q(θ),and[U−1(θ)]nYis the submatrix of [U−1(θ)]whose rows are the corresponding left-eigenvectors. 123
398 SERIEs (2023) 14:379–416 and the second column is n×var( θ) (with var( θ) estimated by Monte Carlo), to be compared with V∗in Table 1. The columns RMSE and CI0.95 are the finite-nRMSE and coverage rate. All the results are close to those in Table 1, confirming that large-n asymptotics provide a good approximation to the finite-ndistribution of θ. Note that Table 2does not report simulation results for q=∞. This is because in some Monte Carlo runs, in particular for T=6, it turned out to be too difficult to numerically distinguish between the smallest and the second smallest eigenvalue of Q(x,θ)and, therefore, to reliably select the eigenvector associated with the smallest eigenvalue. This is another reason not to recommend choosing q=∞. 7.2 Case 2: binary regressor, T0=1,T1=T−1 There is nothing special about the case T0=T1=T/2, which we just discussed. Any other T0≥1 and T1=T−T0≥1 lead to qualitatively similar results. We illustrate this for the case T0=1 and T1=T−1. Table 3is similar to Table 1and reports results for (T0,T1)=(1,T−1)with T∈{4,10}. With T0=1 fixed, we find that the bias for q=0 is nearly constant in T(and large), while the bias of the bias-corrected estimates (q>0)decreases in T. Again, the bias is not monotonic in q(it changes sign) and it becomes very small as qbecomes sufficiently large, albeit more slowly than in the case T0=T1=T/2. Table 4presents the corresponding simulations, showing that the finite-sample results are, again, very close to asymptotic results reported in Table 3. 7.3 Case 3: continuous regressor Here we illustrate approximate functional differencing in a panel probit model with a single continuous regressor and fixed effects. Apart from the continuity of the regressor, the setup is identical to that in Cases 1 and 2 above. Specifically, we consider Example 1A with Uit ∼N(0,1),θ0=1, πprior(α) =φ(α),π0(α) =φ(α −1), and the integrated score as initial score. We set Xit ∼N(0.5,0.25), so that Xit has the same mean and variance (across t) as the binary regressor in Cases 1 and 2. For T∈{4,6} and n=1000, we generated a single data set Xit (t=1,...,T;i=1,...,n)to form X={X1,...,Xn}, so the results are to be understood with reference to this X. Table 5presents the asymptotic biases and variances for qup to 1000. (We do not consider q=∞here because, even though θ∗and θremain well-defined in the limit q→∞, the limiting values are generically determined by a single Xi∈X, that is, estimation would be based on a single observation i,forwhichQ(Xi,θ)has the smallest eigenvalue within the sample.) The results are similar to those in Table 1, albeit the asymptotic biases and variances are somewhat larger (for all q). Table 6 presents the corresponding simulation results, which are in line with those in Table 5. 7.4 Numerical calculations related to Conjecture 1 Table 7reports bias calculations for larger values of Tthat support our conjecture about the rate of the bias as Tgrows. The model is as in Example 1B with standard normal 123
SERIEs (2023) 14:379–416 399 Table 1 Asymptotic biases and variances, and approximate RMSEs and coverage rates (n=1000) qT 0=T1=2,T=4T0=T1=3,T=6 θ∗−θ0V∗RMSE CI0.95 θ∗−θ0V∗RMSE CI0.95 0 0.5050 3.3313 0.5083 0.0000 0.4056 2.3063 0.4084 0.0000 1 0.1525 3.3116 0.1630 0.2452 0.0787 2.3477 0.0924 0.6315 2−0.0039 3.4940 0.0592 0.9495 −0.0172 2.5114 0.0530 0.9364 3−0.0513 3.6583 0.0793 0.8644 −0.0321 2.6289 0.0605 0.9041 4−0.0577 3.8016 0.0844 0.8453 −0.0281 2.7157 0.0592 0.9160 5−0.0516 3.9307 0.0812 0.8694 −0.0221 2.7789 0.0572 0.9296 6−0.0433 4.0438 0.0769 0.8955 −0.0175 2.8231 0.0559 0.9375 7−0.0358 4.1383 0.0736 0.9139 −0.0144 2.8532 0.0553 0.9416 8−0.0297 4.2145 0.0714 0.9256 −0.0124 2.8735 0.0550 0.9438 9−0.0252 4.2745 0.0701 0.9328 −0.0111 2.8873 0.0549 0.9451 10 −0.0218 4.3210 0.0692 0.9373 −0.0102 2.8969 0.0548 0.9459 20 −0.0124 4.4722 0.0680 0.9460 −0.0070 2.9262 0.0545 0.9481 40 −0.0104 4.5149 0.0680 0.9473 −0.0042 2.9365 0.0544 0.9493 60 −0.0091 4.5311 0.0679 0.9479 −0.0031 2.9392 0.0543 0.9496 80 −0.0080 4.5415 0.0679 0.9484 −0.0026 2.9404 0.0543 0.9497 100 −0.0071 4.5495 0.0678 0.9487 −0.0024 2.9411 0.0543 0.9498 200 −0.0046 4.5708 0.0678 0.9495 −0.0020 2.9435 0.0543 0.9498 400 −0.0033 4.5793 0.0678 0.9497 −0.0017 2.9464 0.0543 0.9499 600 −0.0031 4.5819 0.0678 0.9498 −0.0015 2.9487 0.0543 0.9499 800 −0.0030 4.5849 0.0678 0.9498 −0.0014 2.9508 0.0543 0.9499 1000 −0.0030 4.5881 0.0678 0.9498 −0.0012 2.9528 0.0544 0.9499 ∞−0.0452 19.2259 0.1387 0.9500 −0.01025 13.9013 0.1179 0.9500 The model is as in Example 1B with standard normal errors, πprior(α) =φ(α),π0(α) =φ(α −1),andθ0=1. The pseudo-true value θ∗is based on the integrated score. Notation: 0r=00 ...0 rzeros , e.g., −0.0452 =−0.000052 123
400 SERIEs (2023) 14:379–416 Table 2 Simulation results for n=1000. Setup as in Table 1. Results based on 1000 Monte Carlo replications qT 0=T1=2,T=4T0=T1=3,T=6 Bias n×var RMSE CI0.95 Bias n×var RMSE CI0.95 00.5067 3.5931 0.5102 0.0000 0.4076 2.3947 0.4105 0.000 10.1543 3.5243 0.1653 0.2260 0.0807 2.4378 0.0946 0.614 2−0.0020 3.6881 0.0608 0.9400 −0.0151 2.6022 0.0532 0.936 3−0.0493 3.8578 0.0793 0.8530 −0.0299 2.7221 0.0601 0.902 4−0.0556 4.0160 0.0843 0.8300 −0.0259 2.8112 0.0590 0.912 5−0.0494 4.1629 0.0813 0.8590 −0.0198 2.8760 0.0572 0.924 6−0.0409 4.2930 0.0773 0.8830 −0.0151 2.9209 0.0561 0.937 7−0.0333 4.4024 0.0742 0.9080 −0.0120 2.9512 0.0556 0.946 8−0.0272 4.4908 0.0723 0.9190 −0.0100 2.9712 0.0554 0.952 9−0.0225 4.5604 0.0712 0.9270 −0.0087 2.9845 0.0553 0.953 10 −0.0191 4.6142 0.0706 0.9340 −0.0078 2.9934 0.0553 0.951 20 −0.0096 4.7858 0.0698 0.9380 −0.0046 3.0148 0.0551 0.955 40 −0.0075 4.8264 0.0699 0.9390 −0.0018 3.0159 0.0549 0.956 60 −0.0062 4.8383 0.0698 0.9370 −0.0007 3.0150 0.0549 0.957 80 −0.0051 4.8446 0.0698 0.9380 −0.0003 3.0146 0.0549 0.957 100 −0.0043 4.8489 0.0698 0.9380 −0.0001 3.0145 0.0549 0.958 200 −0.0018 4.8588 0.0697 0.9420 0.0003 3.0150 0.0549 0.957 400 −0.0005 4.8603 0.0697 0.9440 0.0006 3.0151 0.0549 0.957 600 −0.0003 4.8611 0.0697 0.9440 0.0008 3.0152 0.0549 0.957 800 −0.0003 4.8630 0.0697 0.9450 0.0009 3.0153 0.0549 0.957 1000 −0.0002 4.8652 0.0698 0.9450 0.0010 3.0156 0.0549 0.959 123
SERIEs (2023) 14:379–416 401 Table 3 Asymptotic biases and variances, and approximate RMSEs and coverage rates (n=1000) qT 0=1,T1=3,T=4T0=1,T1=9,T=10 θ∗−θ0V∗RMSE CI0.95 θ∗−θ0V∗RMSE CI0.95 00.6704 2.7135 0.6725 0.0000 0.6721 1.4919 0.6732 0.0000 10.3131 3.0434 0.3179 0.0001 0.2545 2.1031 0.2586 0.0002 20.0758 3.6051 0.0967 0.7564 0.0567 2.7400 0.0771 0.8087 3−0.0252 3.9220 0.0675 0.9312 0.0026 2.9696 0.0546 0.9497 4−0.0576 4.0648 0.0859 0.8525 −0.0122 3.0664 0.0567 0.9444 5−0.0641 4.1438 0.0908 0.8310 −0.0172 3.1265 0.0585 0.9391 6−0.0623 4.2037 0.0899 0.8395 −0.0192 3.1701 0.0595 0.9366 7−0.0583 4.2566 0.0875 0.8545 −0.0200 3.2030 0.0600 0.9356 8−0.0544 4.3047 0.0852 0.8683 −0.0202 3.2285 0.0603 0.9354 9−0.0510 4.3484 0.0833 0.8793 −0.0201 3.2485 0.0604 0.9356 10 −0.0481 4.3877 0.0819 0.8877 −0.0199 3.2646 0.0605 0.9360 20 −0.0354 4.6294 0.0767 0.9184 −0.0164 3.3403 0.0601 0.9408 40 −0.0281 4.8096 0.0748 0.9310 −0.0130 3.3925 0.0597 0.9442 60 −0.0250 4.8776 0.0742 0.9351 −0.0114 3.4177 0.0596 0.9456 80 −0.0231 4.9177 0.0738 0.9375 −0.0104 3.4340 0.0595 0.9464 100 −0.0217 4.9491 0.0736 0.9391 −0.0097 3.4462 0.0595 0.9469 200 −0.0173 5.0518 0.0732 0.9432 −0.0078 3.4839 0.0595 0.9480 400 −0.0147 5.1243 0.0731 0.9452 −0.0065 3.5246 0.0597 0.9486 600 −0.0141 5.1505 0.0731 0.9455 −0.0058 3.5524 0.0599 0.9489 800 −0.0139 5.1726 0.0732 0.9457 −0.0054 3.5744 0.0600 0.9491 1000 −0.0137 5.1968 0.0734 0.9459 −0.0051 3.5935 0.0602 0.9492 ∞−0.0378 15.7761 0.1256 0.9500 −0.0620 35.4401 0.1883 0.9500 The model is as in Example 1B with standard normal errors, πprior(α) =φ(α),π0(α) =φ(α −1),andθ0=1. The pseudo-true value θ∗is based on the integrated score. Notation: 0r=00 ...0 rzeros 123
402 SERIEs (2023) 14:379–416 Table 4 Simulation results for n=1000 qT 0=1,T1=3,T=4T0=1,T1=9,T=10 Bias n×var RMSE CI0.95 Bias n×var RMSE CI0.95 00.6721 2.6261 0.6740 0.0000 0.6740 1.4847 0.6751 0.0000 10.3147 3.1535 0.3197 0.0000 0.2559 2.1953 0.2602 0.0000 20.0774 3.8340 0.0991 0.7500 0.0576 2.8768 0.0787 0.7950 3−0.0239 4.1970 0.0690 0.9160 0.0034 3.1192 0.0560 0.9470 4−0.0562 4.3578 0.0867 0.8480 −0.0115 3.2218 0.0579 0.9420 5−0.0627 4.4455 0.0915 0.8320 −0.0166 3.2857 0.0597 0.9360 6−0.0608 4.5111 0.0906 0.8360 −0.0186 3.3324 0.0606 0.9330 7−0.0568 4.5686 0.0883 0.8550 −0.0194 3.3679 0.0612 0.9310 8−0.0528 4.6209 0.0861 0.8630 −0.0196 3.3955 0.0615 0.9310 9−0.0493 4.6682 0.0843 0.8720 −0.0195 3.4172 0.0616 0.9320 10 −0.0464 4.7110 0.0829 0.8810 −0.0192 3.4347 0.0617 0.9320 20 −0.0334 4.9732 0.0780 0.9050 −0.0157 3.5184 0.0613 0.9370 40 −0.0258 5.1642 0.0764 0.9130 −0.0122 3.5780 0.0610 0.9410 60 −0.0226 5.2313 0.0758 0.9190 −0.0105 3.6080 0.0610 0.9420 80 −0.0206 5.2686 0.0754 0.9210 −0.0094 3.6283 0.0610 0.9430 100 −0.0190 5.2971 0.0752 0.9260 −0.0087 3.6438 0.0610 0.9430 200 −0.0145 5.3904 0.0748 0.9350 −0.0066 3.6935 0.0611 0.9420 400 −0.0117 5.4577 0.0748 0.9330 −0.0051 3.7479 0.0614 0.9430 600 −0.0111 5.4847 0.0749 0.9350 −0.0042 3.7836 0.0617 0.9440 800 −0.0108 5.5093 0.0750 0.9350 −0.0037 3.8109 0.0618 0.9480 1000 −0.0106 5.5370 0.0752 0.9340 −0.0033 3.8338 0.0620 0.9470 Setup as in Table 3. Results based on 1000 Monte Carlo replications 123
SERIEs (2023) 14:379–416 409 Using the posterior distribution in (11), a natural baseline estimating function (q=0) is w(0)(y,x,θ):= A μ(x,α,θ)π post(α |y,x,θ)dα. The corresponding estimator μ(0)of μ0can again be motivated by “large-T” panel data considerations, where, under regularity conditions, the posterior distribution concentrates around the true value Aas T→∞. Let W(x,θ)be the nY-vector with entries w(0)(y(k),x,θ),k=1,...,nY. Then, the analog of the limiting estimating function in (20), corresponding to q→∞,for average effects is w(∞)(y,x,θ):= W(x,θ)Q†(x,θ)δ(y) =W(x,θ) h(∞)[Q(x,θ)]δ(y), h(∞)(λ) := λ−1for λ>0, 0forλ=0, (24) where Q†(x,θ) is a pseudo-inverse of Q(x,θ), and the application of a function h(∞):[0,1]→Rto the matrix Q(x,θ) was defined in equation (18). The motivation for choosing w(∞)(y,x,θ) in this way is that it gives an unbiased estimator of the average effect (i.e., μ(∞) ∗=μ0) whenever we can write μ(x,α,θ)=y∈Yν(y,x,θ)fyx,α,θfor some function ν(y,x,θ).7Of course, average effects with this form of μ(α, x,θ)are a very special case, but they are usually the only cases for which we can expect unbiased estimation of the average effect to be feasible (for fixed T); see also Aguirregabiria and Carro (2021). Notice that we do not assume here that μ(α, x,θ)is of this form, it is just used to motivate (24). As we have seen before, the non-zero eigenvalues of Q(x,θ) can be very small, which implies that the pseudo-inverse Q†(x,θ) can have very large elements. The corresponding estimator μ(∞)based on (24) therefore typically has a very large variance and we do not recommend this estimator in practice. Instead, to balance the bias-variance trade-off of the average effect estimator, some regularization of the pseudo-inverse of Q(x,θ) in (24) is required. There are various ways to implement regularization, in the same way that there are various ways to implement approximate functional differencing (see Sect.8.2). 7This is because in that special case we have W(x,θ) =N(x,θ)Q(x,θ),whereN(x,θ) is the nY-vector with entries ν(y(k),x,θ), and therefore Ew(∞)(Y,X,θ 0)X=x,A=α= N(x,θ 0)Q(x,θ 0)Q†(x,θ 0)Eδ(Y)X=x,A=α=N(x,θ 0)Eδ(Y)X=x,A=α= μ(x,α,θ 0). 123
410 SERIEs (2023) 14:379–416 Here, regularization means that we want to find functions hq(λ) that approximate the inverse function 1/λ well for large values of λ∈[0,1], but that deviate from 1/λ for values of λclose to zero to avoid divergence.8This gives,9for q∈{0,1,2,...}, hq(λ) = q r=0 (1−λ)r=1−(1−λ)q+1 λfor λ>0, q+1forλ=0.(25) The corresponding estimating function that regularizes w(∞)(y,x,θ) is therefore given by w(q)(y,x,θ):= W(x,θ) hq[Q(x,θ)]δ(y), hq[Q(x,θ)]= q r=0InY−Q(x,θ) r. (26) This is a polynomial in Q(x,θ), as was the case for s(q)(y,x,θ). Choosing a value of q that is not too large therefore ensures that the variance of the corresponding estimator μ(q)remains reasonably small (for fixed q), because we don’t need the pseudo-inverse of Q(x,θ). Note also that w(q)(y,x,θ) and the corresponding estimators μ(q)have a largeTbias-correction interpretation very similar to s(q)(y,x,θ). For example, we have h1(λ) =2−λ, and therefore w(1)(y,x,θ)=2w(0)(y,x,θ)−W(x,θ) Q(x,θ)δ(y). We conjecture that the estimator of μ0corresponding to only W(x,θ)Q(x,θ)δ(y) has twice the leading order 1/Tasymptotic bias of the estimator μ(0)corresponding to w(0)(y,x,θ), that is, w(1)(y,x,θ)is exactly the jackknife linear combination that eliminates the large-Tleading order bias in μ(0); see Dhaene and Jochmans (2015b). Appropriate iterations of this jackknife bias correction also give the estimating functions w(q)(y,x,θ)for q>1. We are not considering average effects further here. But we found it noteworthy that there is a formalism for average effect calculation that closely mirrors the development of approximate functional differencing for the estimation of θ0introduced 8In previous sections, the functions hq(λ) =(1−λ)qwere polynomial approximations of (rescaled versions of) the function h∞(λ) =1{λ=0}. The regularization that is analogous to s(q)(y,x,θ)in (17) isgivenbyaq-th order Taylor expansion of the function 1/λ around λ=1. 9Here,weusetheconventionthat0 0=1, which also implies that InY−Q(x,θ) 0=InYeven though Q(x,θ) has an eigenvalue equal to one. Also, there is some ambiguity in what value we should assign to hq(λ) for λ=0. We choose hq(0)=q+1 because it results in the simple polynomial expression (26)for hq[Q(x,θ)], which is convenient since hq[Q(x,θ)]can be evaluated without ever calculating the eigenvalues and eigenvectors of Q(x,θ). However, if we want to obtain w(∞)(y,x,θ)in (24) as the limit of w(q)(y,x,θ)as q→∞, then we should assign hq(0)=0forλ=0, but this would deviate from the polynomial expression. 123
SERIEs (2023) 14:379–416 411 above. However, this does not imply that we expect the results for average-effect estimation to be necessarily similar to those for the estimation of the common parameters θ0. In particular, for small values of T, the identified set for the average effects in discrete-choice panel data models tends to be much larger than the identified set of the common parameters (see, e.g., Chernozhukov, Fernández-Val, Hahn and Newey 2013; Davezies, D’Haultfoeuille and Laage 2021; Liu, Poirier and Shiu 2021; Pakel and Weidner 2021). Therefore we expect larger values of Tto be required for the point estimators μ(q)to perform well, and we also expect the bias-variance trade-off in the choice of qto be quite different. For a closely related discussion see Bonhomme and Davezies (2017), and also the section on “Average marginal effects” in the 2010 working paper version of Bonhomme (2012). 9 Conclusions We have linked the large-Tpanel data literature with the functional differencing method through a bias correction that converges to functional differencing when iterated. Our numerical illustrations show that in models where exact functional differencing is not possible, one may still apply it approximately to obtain estimates that can be essentially unbiased, even when the number of time periods Tis small. The key element in our construction is the nY×nYmatrix Q(x,θ). The eigenvalues of this matrix are informative about whether (approximate) functional differencing is applicable in a given model. The matrix Q(x,θ) also features prominently in our bias-corrected score functions in (17) and in our regularized estimating functions for average effects in (26). We have assumed a discrete outcome space with a finite number of elements nY. When the outcome space is infinite, the matrix Q(x,θ)has to be replaced by the corresponding operator. The goal of this paper was primarily to introduce and illustrate an approximate version of functional differencing. Future work is needed to better understand the properties of the method and to explore its usefulness in empirical work, both for the estimation of common parameters, which was our primary focus, and for the estimation of average effects, briefly introduced in Sect.8.3. Data availability Data sharing not applicable to this article as no datasets were generated or analysed during the current study. Declarations Conflict of interest We have no financial or non-financial interests to disclose. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. 123
412 SERIEs (2023) 14:379–416 A Proofs Proof of Lemma 1Define Q(y|y,x,θ)=Afyx,α,θf(y|x,α,θ)π prior(α |x)dα pprior(y|x,θ) 1/2pprior(y|x,θ) 1/2 and let Q(x,θ)be the nY×nYmatrix with elements Qk,(x,θ)=Q(y(k)|y(),x,θ). Also define the nY×nYdiagonal matrix Pprior(x,θ)=diag pprior(y(k)|x,θ) k=1,...,nY. From (11) and (14) we obtain Q(y|y,x,θ)=Afyx,α,θf(y|x,α,θ)π prior(α |x)dα pprior(y|x,θ) =pprior(y|x,θ) 1/2Q(y|y,x,θ) pprior(y|x,θ) −1/2, which in matrix notation is Q(x,θ)=Pprior(x,θ) 1/2Q(x,θ) Pprior(x,θ) −1/2. This shows that the matrices Q(x,θ)and Q(x,θ)are similar and therefore have the same eigenvalues.10 The matrix Q(x,θ)is symmetric and positive semi-definite (by construction), which implies that all its eigenvalues (and therefore all eigenvalues of Q(x,θ)) are non-negative real numbers. Furthermore, Q(x,θ) is diagonalizable because it is symmetric. Hence Q(x,θ)is also diagonalizable, because it is similar to Q(x,θ).11 In addition, Q(x,θ)is a stochastic matrix (by construction), which implies that its spectral radius is equal to one, that is, Q(x,θ)cannot have any eigenvalue larger than one. We thus conclude that all eigenvalues of Q(x,θ)lie in the interval [0,1]. The following lemma is useful for the proof of Lemma 2, which we present afterward. Lemma 3 Let the assumptions of Lemma 2hold. Let w(y,x,θ 0)∈Rbe such that y∈Y w(y,x,θ 0)Qyy,x,θ 0=0,for all y∈Y. 10 Two matrices Aand Bare similar if B=P−1AP for some nonsingular matrix P. Similar matrices have the same eigenvalues. 11 A matrix is diagonalizable if and only if it is similar to a diagonal matrix. Since Q(x,θ)is similar to a diagonal matrix, and Q(x,θ) is similar to Q(x,θ), it must also be the case that Q(x,θ)is similar to a diagonal matrix. 123
SERIEs (2023) 14:379–416 413 Then y∈Y w(y,x,θ 0)fyx,α,θ 0=0,for all α∈A. Proof of Lemma 3The nY×nYdiagonal matrix Pprior(x,θ 0)was defined in the Proof of Lemma 1. In addition, let F(x,α,θ 0)and W(x,θ 0)be the nY-vectors with elements f(y(k)|x,α,θ 0)and w(y(k),x,θ 0), respectively, for k=1,...,nY. Then Q(x,θ 0)=A F(x,α,θ 0)F(x,α,θ 0)π prior(α |x)dαP−1 prior(x,θ 0), (27) and the condition on w(y,x,θ 0)in the lemma can be written as W(x,θ 0)Q(x,θ 0)=0. Plugging in the expression for Q(x,θ 0)in (27) and multiplying with Pprior(x,θ 0) W(x,θ 0)from the right gives A W(x,θ 0)F(x,α,θ 0)F(x,α,θ 0)W(x,θ 0)π prior(α |x)dα=0. Since W(x,θ 0)F(x,α,θ 0)F(x,α,θ 0)W(x,θ 0)≥0 and πprior(α |x)>0 we conclude that W(x,θ 0)F(x,α,θ 0)F(x,α,θ 0)W(x,θ 0)=0,(28) for almost all values α, except possibly for a set of values αthat has measure zero under πprior(α |x). However, since fyx,α,θ 0is assumed to be continuous in α, we conclude that (28) must hold for all α∈A, since any violation on a set of measure zero would require a discontinuity in α. Finally, (28) also implies that W(x,θ 0)F(x,α,θ 0)=0, for all α∈A. This is what we wanted to show, just written in vector notation. Proof of Lemma 2# part (i): Let U0(x,θ 0)be the submatrix of U(x,θ 0)that only contains those columns that are the right-eigenvectors of Q(x,θ 0)corresponding to the eigenvalues λj(x,θ 0)=0. We then have Q(x,θ 0)U0(x,θ 0)=0.Similarly, let [U−1(x,θ 0)]0be the submatrix of U−1(x,θ 0)that only contains the rows that are the left-eigenvectors of Q(x,θ 0)corresponding to the eigenvalues λj(x,θ 0)=0. We then have [U−1(x,θ 0)]0Q(x,θ 0)=0, 123
414 SERIEs (2023) 14:379–416 and according to Lemma 3this implies [U−1(x,θ 0)]0F(x,α,θ 0)=0.(29) Next, by using the definition of s∞(y,x,θ 0)and h∞[Q(x,θ 0)]in the main text we find s∞(y,x,θ 0)=S(x,θ 0)h∞[Q(x,θ 0)]δ(y) =S(x,θ 0)U(x,θ 0)diag 1λj(x,θ 0)=0j=1,...,nY!U−1(x,θ 0)δ(y) =S(x,θ 0)U0(x,θ 0)[U−1(x,θ 0)]0δ(y), and therefore Es∞(Y,X,θ 0)X=x,A=α=S(x,θ 0)U0(x,θ 0)[U−1(x,θ 0)]0F(x,α,θ 0)=0, where in the last step we used (29). # part (ii): Let m(x,θ 0)be the nY-vector with elements m(y(k),x,θ 0),k= 1,...,nY. Then, Em(Y,X,θ 0)X=x,A=α=0 can be written in vector notation as m(x,θ 0)F(x,α,θ 0)=0.(30) From the expression of Q(x,θ 0)in (27) we see that this implies m(x,θ 0)Q(x,θ 0)= 0, that is, if (30) holds for all α∈A, then Q(x,θ 0)has a zero eigenvalue with corresponding left-eigenvector m(x,θ 0). This is the “if” part of the statement in part (ii) of the lemma. Conversely, if Q(x,θ 0)has a zero eigenvalue, then let m(x,θ 0)be a corresponding left-eigenvector. We then have m(x,θ 0)Q(x,θ 0)=0. According to Lemma 3this implies that (30) holds, or equivalently that Em(Y,X,θ 0)X=x,A=α=0. We have thus also shown the “only if” part of the statement in part (ii) of the lemma. # part (iii): Let m(y,x,θ 0)∈Rbe such that Em(Y,X,θ 0)X=x,A=α=0. We choose s(y,x,θ 0)=m(y,x,θ 0). Using the definition of s(1)(y,x,θ)in (16) we then find s(1)(y,x,θ 0)=m(y,x,θ 0), and therefore also s(q)(y,x,θ 0)=m(y,x,θ 0), for all q∈{1,2,...}. We therefore also find s∞(y,x,θ 0)=limq→∞ s(q)(y,x,θ 0)= m(y,x,θ 0), which is what we wanted to show. 123
SERIEs (2023) 14:379–416 415 References Aguirregabiria V, Carro JM (2021) Identification of average marginal effects in fixed effects dynamic discrete choice models. arXiv preprint arXiv:2107.06141 Alvarez J, Arellano M (2003) The time series and cross-section asymptotics of dynamic panel data estimators. Econometrica 71(4):1121–1159 Andersen EB (1970) Asymptotic properties of conditional maximum-likelihood estimators. J R Stat Soc Ser B (Methodol) 32(2):283–301 Arellano M (2003) Discrete choices with panel data. Investig Econ 27(3):423–458 Arellano M, Bond S (1991) Some tests of specification for panel data: Monte Carlo evidence and an application to employment equations. Rev Econ Stud 58(2):277–297 Arellano M, Bonhomme S (2009) Robust priors in nonlinear panel data models. Econometrica 77(2):489– 536 Arellano M, Bonhomme S (2011) Nonlinear panel data analysis. Annu Rev Econ 3(1):395–424 Arellano M, Bonhomme S (2012) Identifying distributional characteristics in random coefficients panel data models. Rev Econ Stud 79(3):987–1020 Arellano M, Hahn J (2007) Understanding bias in nonlinear panel models: some recent developments. Econom Soc Monogr 43:381 Arellano M, Hahn J (2016) A likelihood-based approximate solution to the incidental parameter problem in dynamic nonlinear models with multiple effects. Glob Econ Rev 45(3):251–274 Bonhomme S (2012) Functional differencing. Econometrica 80(4):1337–1385 Bonhomme S, Manresa E (2015) Grouped patterns of heterogeneity in panel data. Econometrica 83(3):1147– 1184 Bonhomme S, Weidner M (2022) Minimizing sensitivity to model misspecification. Quant Econ 13(3):907– 954 Bonhomme S, Lamadon T, Manresa E (2022) Discretizing unobserved heterogeneity. Econometrica 90(2):625–643 Bonhomme S, Davezies L (2017) Panel data, inverse problems, and the estimation of policy parameters Chamberlain G (1980) Analysis of covariance with qualitative data. Rev Econ Stud 47(1):225–238 Chamberlain G (2010) Binary response models for panel data: identification and information. Econometrica 78(1):159–168 Chernozhukov V, Fernández-Val I, Hahn J, Newey W (2013) Average and quantile effects in nonseparable panel models. Econometrica 81(2):535–580 Davezies L, D’Haultfoeuille X, Laage L (2021). Identification and estimation of average marginal effects in fixed effect logit models. arXiv preprint arXiv:2105.00879 Davezies L, D’Haultfoeuille X, Mugnier M (2022) Fixed effects binary choice models with three or more periods. Quant Econ (forthcoming) Dhaene G, Jochmans K (2015) Split-panel Jackknife estimation of fixed-effect models. Rev Econ Stud 82(3):991–1030 Dhaene G, Jochmans K (2015a) Profile-score adjustments for incidental-parameter problems. Working Paper, Sciences Po, Paris Google Scholar Article Location Dobronyi C, Gu J, Kim KI (2021) Identification of dynamic panel logit models with fixed effects. arXiv preprint arXiv:2104.04590 Fernández-Val I, Weidner M (2016) Individual and time effects in nonlinear panel models with large n, t. J Econom 192(1):291–312 Hahn J, Kuersteiner G (2002) Asymptotically unbiased inference for a dynamic panel model with fixed effects when both n and t are large. Econometrica 70(4):1639–1657 Hahn J, Newey W (2004) Jackknife and analytical bias reduction for nonlinear panel models. Econometrica 72(4):1295–1319 Hansen LP (1982) Large sample properties of generalized method of moments estimators. Econom J Econom Soc 45:1029–1054 Honoré BE, Muris C, Weidner M (2021) Dynamic ordered panel logit models. arXiv preprint arXiv:2107.03253 Honoré BE, Weidner M (2020) Moment conditions for dynamic panel logit models with fixed effects. arXiv preprint arXiv:2005.05942 Honoré BE (1992) Trimmed LAD and least squares estimation of truncated and censored regression models with fixed effects. Econometrica 60:533–565 123
416 SERIEs (2023) 14:379–416 Honoré BE, Tamer ET (2006) Bounds on parameters in panel dynamic discrete choice models. Econometrica 74(3):611–629 Horn RA, Johnson CR (1994) Topics in matrix analysis Hu L (2002) Estimation of a censored dynamic panel data model. Econometrica 70(6):2499–2517 Johnson EG (2004) Identification in discrete choice models with fixed effects. Working paper, Bureau of Labor Statistics. CiteSeer Kitazawa Y (2013) Exploration of dynamic fixed effects logit models from a traditional angle. Technical report, No. 60, Kyushu Sangyo University Faculty of Economics Lancaster T (2000) The incidental parameter problem since 1948. J Econom 95(2):391–413 Lancaster T (2002) Orthogonal parameters and panel data. Rev Econ Stud 69:647–66 Liu L, Poirier A, Shiu J-L (2021) Identification and estimation of average partial effects in semiparametric binary response panel models. Working paper Moreira MJ (2009) A maximum likelihood method for the incidental parameter problem. Ann Stat 37(6A):3660–3696 Neyman J, Scott EL (1948) Consistent estimates based on partially consistent observations. Econometrica 36:1–32 Pakel C, Weidner M (2021) Bounds on average effects in discrete choice panel data models. Technical report, Working paper Rasch G (1960) Studies in mathematical psychology: I. Nielsen & Lydiche, Probabilistic models for some intelligence and attainment tests Su L, Shi Z, Phillips PC (2016) Identifying latent structures in panel data. Econometrica 84(6):2215–2264 Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. 123