EconBase
← Back to paper

Binary Response Models for Heterogeneous Panel Data with Interactive Fixed Effects

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.

73,832 characters · 11 sections · 125 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

\newtheorem{corollary}{Corollary} \newtheorem{definition}{Definition} \newtheorem{lemma}{Lemma} \newtheorem{proposition}{Proposition} \newtheorem{remark}{Remark} \newtheorem{theorem}{Theorem} \newtheorem{assumption}{Assumption}

\numberwithin{corollary}{section} \numberwithin{definition}{section} \numberwithin{equation}{section} \numberwithin{lemma}{section} \numberwithin{proposition}{section} \numberwithin{remark}{section} \numberwithin{theorem}{section}

\allowdisplaybreaks[4]

titlepage{ \begin{center} { \bf Binary Response Models for Heterogeneous Panel Data with Interactive Fixed Effects \begingroup \footnote{ \begin{itemize} • Correspondence: Bin Peng, Department of Econometrics and Business Statistics, Monash University, Caulfield East, VIC 3145, Australia. Email: [email removed] \end{itemize} } \addtocounter{footnote}{-1} \endgroup } {\sc Jiti Gao$^\sharp$ and Fei Liu$^{\ast}$ and Bin Peng$^{\sharp}$ and Yayi Yan$^{\sharp}$} $^\sharp$Monash University, Australia and $^{\ast}$Nankai University, China \today \begin{abstract} In this paper, we investigate binary response models for heterogeneous panel data with interactive fixed effects by allowing both the cross-sectional dimension and the temporal dimension to diverge. From a practical point of view, the proposed framework can be applied to predict the probability of corporate failure, conduct credit rating analysis, etc. Theoretically and methodologically, we establish a link between a maximum likelihood estimation and a least squares approach, provide a simple information criterion to detect the number of factors, and achieve the asymptotic distributions accordingly. In addition, we conduct intensive simulations to examine the theoretical findings. In the empirical study, we focus on the sign prediction of stock returns, and then use the results of sign forecast to conduct portfolio analysis. \end{abstract} \end{center} {\em Keywords}: Binary Response, Heterogeneous Panel, Interactive Fixed Effects, Portfolio Analysis {\em JEL classification}: C18, C23, G11 }

Introduction

Varieties of binary response panel data models have been proposed and studied over the past a couple of decades, and earlier developments date back at least to Chamberlain and the references therein. The challenges in the previous studies often arise due to the identification issues caused by short time periods of data and non-closed form estimators (e.g., Manski1987, Chamberlain2010; among others). With the rise and availability of big and rich datasets, recent studies on binary response panel data models gradually shift focuses to the cases where both the cross-sectional dimension and the temporal dimension are allowed to diverge. An excellent review is given in FW2018. Recently, an important strand of the literature is devoted to binary response models with interactive fixed effects (e.g., BL2017, Wang2019, Chen2020). Within these studies, the central questions are a) how to estimate the coefficients together with the factors and the factor loadings? and b) to achieve the optimal efficiency in a), how to determine the number of factors?

For linear additive models, the aforementioned questions have been addressed well in the literature by utilizing different techniques. For example, using principal component analysis (PCA) and random matrix theory, BN2002, Onatski, LamYao2012 and AhnHorenstein2013 are able to detect the number of factors for large panel data using information criteria or eigenanalysis; Pesaran2006 introduces a common correlated effects (CCE) estimator that takes an advantage of a factor structure involving both dependent and independent variables; Bai2009 and Moon develop alternative methods to estimate the coefficients together with the factors and the factor loadings; LCL2020 and HJPS respectively use the maximum likelihood method and the classifier-Lasso method to consider the cases with heterogeneous coefficients; and so forth.

However, for non-linear panel data models, especially for binary response panel data models involving interactive fixed effects, there is limited progress, which is mainly due to the fact that a variety of tools adopted for linear additive models may no longer be directly applicable and useful. Below, we comment on the relevant literature. In BL2017, the authors extend the CCE approach to a framework with binary responses, in which a key step is to estimate unobservable factors from the regressors. As a consequence, the approach requires an explicit structure of the regressors, and the usual limitation of a CCE type estimator occurs, e.g., the number of unobservable factors cannot be larger than the number of regressors (cf., BL2017). Wang2019 and Chen2020 propose similar solutions to binary response panel data models, and the main difference is that the former does not include any regressors in the model. Thereby, we may regard Wang2019 and Chen2020 as the binary response counterparts of BN2002 and Bai2009, respectively. Recently, AB2020, AL2020 and CDG2020 bring attention to additive panel data models with interactive fixed effects, in which closed form estimators are less obvious. Specifically, quantile regressions are investigated in all three papers, of which only AL2020 include regressors, while AB2020 propose a Bayesian approach.

In view of the aforementioned literature, we specifically consider a binary response panel data model with interactive fixed effects by incorporating heterogeneous coefficients. From a modelling perspective, it is similar to BL2017, but we require less structure on the regressors, which allows us to avoid adding any restriction between the number of regressors and the number of unobservable factors. Our investigation establishes a link between a maximum likelihood estimation and a nonlinear least squares approach. As a consequence, some tools adopted for linear additive models immediately become applicable. For example, the identification restrictions provided in Bai2009 and Moon are readily to be applied to the binary response models with very minor modifications. Meanwhile, we are able to estimate the unknown heterogeneous parameters as well as the unobservable factors and factor loadings, and establish the asymptotic distributions accordingly. Our approach may be considered as the binary response counterpart of that considered in BN2013. In addition, we propose a simple information criterion to detect the number of factors. Last but not least, we conduct intensive numerical studies to examine the theoretical findings, and demonstrate the practical relevance.

From a practical perspective, the proposed framework can be relevant and applicable to the following fields for instance. Predicting the probability of corporate failure has gained its attention since the seminal work of Altman. Along this line of research, our paper provides a more generalized framework to extend those panel data driven studies (e.g., CCL2014). Similarly, our model and estimation method can be applied to panel data based credit rating analysis (e.g., JJW2015). In the empirical study of the paper, we pay particular attention to portfolio analysis. More often than not, in order to calculate the optimal weights assigned to each stock, one adopts all stocks to construct a covariance matrix (e.g., CLL2019, ELW2019), which then naturally falls into the category of high dimensional covariance matrix estimation. Thereby, to boost the estimation accuracy, we see the increasing popularity of methods using rank reduction (PX2019), penalization (CLL2019), or both (FLM13) among others. It is worth emphasizing that by default the aforementioned techniques eventually include all stocks in practice, though some of them may have relatively small weights compared to the others. As pointed out in CD2006 and Nyberg2011, the sign of stock market returns may be predictable even if the returns themselves are not predictable. Also, CD2006 mention that “As volatility moves, so too does the probability of a positive return: the higher the volatility, the lower the probability of a positive return". In connection with the fact that a primary goal of portfolio analysis is to minimize the volatility (ELW2019), a binary response panel data model with interactive fixed effects naturally marries the above studies by modelling the probabilities of positive returns, so we can drop those having low probabilities.

In summary, the main contributions of the paper are as follows. (i). We consider a binary panel data model with both heterogeneous coefficients and interactive fixed effects, and establish a link between the maximum likelihood estimation and the nonlinear least squares estimation. As a consequence, the traditional type of identification conditions (such as those in Bai2009 and Moon) for linear panel data models with interactive fixed effects is readily to be applied with very minor modification. (ii). We provide a simply information criterion to select the number of factors. (iii). In addition to extensive simulation studies, we use the newly established model and approach to bridge two strands of studies on stock returns, so that better performance can be achieved for portfolio analysis.

The structure of this paper is as follows. Section (ref) proposes the model, develops the methodology, and then establishes the asymptotic results. Section (ref) provides intensive simulations to examine the finite-sample performance of the theoretical findings. Section (ref) considers an empirical portfolio analysis. Section (ref) concludes. Appendix (ref) sketches the outline of the theoretical development, and comments on the bias correction and the average partial effects. For the sake of space, we only provide the proofs of some selected main results in this appendix. The omitted proofs and the preliminary lemmas are given in the online supplementary Appendix (ref) of the paper.

Before proceeding further, we introduce some mathematical symbols that will be used repeatedly throughout the paper. $\|\cdot\|$ denotes the Euclidean norm of a vector or the Frobenius norm of a matrix; $O(1)$ always stands for a finite positive constant, and may be different at each appearance; $\to_P$ and $\to_D$ stand for convergence in probability and convergence in distribution respectively; $\Pr(A\, | \,B)$ represents the probability of the event $A$ occurring conditional on the event $B$; $E(Y\, | \,X)$ denotes the expectation of the variable $Y$ conditional on the variable $X$; for a matrix $W$ with full column rank, let $M_W = I -P_W$ with $P_W =W(W'W)^{-1}W'$; for a square matrix $W$, $\rho_{\max}(W)$ stands for its largest eigenvalue; $ a\wedge b = \min\{a,b\} $ and $a\vee b= \max\{a,b\} $; for a square matrix $A$, $\rho_{\max}(A)$ returns the maximum eigenvalue.

The Model and Methodology

In this section, we present the model and the methodology with an algorithm for numerical implementation in practice, and establish the associated asymptotic results. Specifically, we provide the basic setup in Section (ref), and present a numerical estimation procedure for practical implementation; Section (ref) summaries the relevant asymptotic results; and Section (ref) considers the selection of the number of factors.

The Setup

The model we consider is a binary response panel data model with interactive fixed effects of the form:

eqnarray[eqnarray omitted — 173 chars of source]

where $i=1,\ldots,N$ and $t=1,\ldots,T$. In the model (ref), we observe the binary dependent variable $y_{it}$ and the $d_\beta\times 1$ explanatory variables $x_{it}$ with $d_\beta$ being finite. For ease of notation, we suppose that $\{\varepsilon_{it}\}$ is an array of identically distributed random errors in $(i,t)$ with the respective probability density function (PDF) and the cumulative distribution function (CDF) being known. Specifically, we denote the PDF and the CDF as $g_\varepsilon(\cdot)$ and $G_\varepsilon(\cdot)$. Both the factor loading $\gamma_{0i}$ and the factor $f_{0t}$ are $d_f\times 1$, where $d_f$ is finite. For the time being, we assume that $d_f$ is known, and we will come back to its estimation with the corresponding asymptotic result in Section (ref) later. In what follows, we are interested in recovering

eqnarray[eqnarray omitted — 171 chars of source]

For notational simplicity, we let $\theta_{0i}=(\beta_{0i}', \gamma_{0i}')'$, and $\Theta_0=(B_0,\Gamma_0)=(\theta_{01},\ldots,\theta_{0N})'$ throughout the paper.

remarkWe now comment on why heteroskedasticity is ruled out in the above setting. When heteroskedasticity occurs, (say, $\varepsilon_{it}\sim N(0 ,\sigma_i^2)$), we can always rewrite the model as follows. \begin{eqnarray} y_{it} = \left\{ \begin{array}{cc} 1, & x_{it}'\beta_{0i}^*+\gamma_{0i}^{*\prime}f_{0t} -\varepsilon_{it}^*\ge 0, \\ 0, & otherwise, \end{array}\right. \end{eqnarray} where $\beta_{0i}^* =\frac{\beta_{0i}}{\sigma_i}$, $\gamma_{0i}^{*} =\frac{\gamma_{0i}}{\sigma_i}$, and $\varepsilon_{it}^* =\frac{\varepsilon_{it}}{\sigma_i} $. The transferred model (ref) indicates that we can only estimate the true parameters up to unknown constants $\sigma_i$'s unless some further restrictions are imposed.

We now start presenting our estimators for (ref). Simple algebra shows that

eqnarray[eqnarray omitted — 257 chars of source]

which immediately yields $E[y_{it} \, | \, x_{it},\gamma_{0i}, f_{0t}] = G_\varepsilon(x_{it}'\beta_{0i}+\gamma_{0i}'f_{0t}).$ Thus, the likelihood function is specified as follows:

eqnarray[eqnarray omitted — 224 chars of source]

where $B =(\beta_1,\ldots, \beta_N)'$, $F = (f_1,\ldots, f_T)'$ and $\Gamma = (\gamma_1,\ldots, \gamma_N)'$ are $N\times d_\beta$, $T\times d_{f}$ and $N\times d_{f}$ matrices respectively. Moreover, for the purpose of identification, we require

eqnarray[eqnarray omitted — 102 chars of source]

Finally, the log-likelihood function is defined below:

eqnarray[eqnarray omitted — 241 chars of source]

where $F \in \mathbf{F} $. The estimators are thus given by

eqnarray[eqnarray omitted — 140 chars of source]

where $\widehat{B} =(\widehat{\beta}_1,\ldots, \widehat{\beta}_N)'$, $\widehat{F} = (\widehat{f}_1,\ldots,\widehat{f}_T)'$, and $\widehat{\Gamma} = (\widehat{\gamma}_1,\ldots, \widehat{\gamma}_N)'$.

In what follows, the main goals are studying the asymptotic behaviours of the estimators of (ref). We sketch the strategy of our asymptotic analysis. Before examining any estimator, we first use the Taylor expansion to investigate the log-likelihood function, which allows us to bridge the maximum likelihood estimation of (ref) and a nonlinear least squares approach in Lemma (ref) below. As a consequence, the identification restrictions provided in Bai2009 and Moon are readily to be applied with very minor modifications. On this point, some examples are provided in Section (ref) for the purpose of demonstration. Afterwards, the rates of converges and the asymptotic distributions can be established accordingly.

Up to this point, it is worth commenting on the practical implementation of (ref). A variety of algorithms have been proposed for the panel data models with interactive fixed effects, especially for the cases where the nonclosed form estimators are involved, e.g., Ando, Wang2019, Chen2020, AL2020, just to name a few. Our numerical implementation largely follows these methods.

enumerate• Initial $\widehat{F}^{(0)}$ by using some random number generators, e.g., the standard normal distribution for each element of $\widehat{F}^{(0)}$. Implement the SVD decomposition on $\widehat{F}^{(0)}$, and update $\widehat{F}^{(0)}$ to ensure $\frac{1}{T}\widehat{F}^{(0)'} \widehat{F}^{(0)}=I$. • For Step $j\ (\ge 1)$, obtain $\widehat{B}^{(j)} =(\widehat{\beta}_{1}^{(j)},\ldots, \widehat{\beta}_{N}^{(j)})'$ and $\widehat{\Gamma}^{(j)} =(\widehat{\gamma}_{1}^{(j)},\ldots, \widehat{\gamma}_{N}^{(j)})'$ by maximizing \begin{eqnarray*} \log L^{(j)} (\beta_i,\gamma_i) &=& \sum_{t=1}^T\Big\{(1-y_{it})\log\left[1-G_\varepsilon(x_{it}'\beta_i+\gamma_{i}'\widehat{f}_{t}^{(j-1)})\right] \\ &&+y_{it} \log G_\varepsilon(x_{it}'\beta_i+\gamma_{i}'\widehat{f}_{t}^{(j-1)}) \Big\} \end{eqnarray*} over all $i\ge 1$, where $\widehat{F}^{(j-1)} =(\widehat{f}_1^{(j-1)},\ldots, \widehat{f}_T^{(j-1)})'$ is obtained from the Step $j-1$. Then obtain $\widehat{F}^{(j)} =(\widehat{f}_{1}^{(j)},\ldots, \widehat{f}_{N}^{(j)})'$ by maximizing \begin{eqnarray*} \log L^{(j)} (f_t) &=& \sum_{i=1}^N\Big\{(1-y_{it})\log\left[1-G_\varepsilon(x_{it}'\widehat{\beta}_i^{(j)}+\widehat{\gamma}_{i}^{(j)\prime}f_{t})\right] \\ &&+y_{it} \log G_\varepsilon(x_{it}'\widehat{\beta}_i^{(j)}+\widehat{\gamma}_{i}^{(j)\prime}f_{t}) \Big\} \end{eqnarray*} over all $t\ge 1$. Finally, implement the SVD decomposition on $\widehat{F}^{(j)}$, and update $\widehat{F}^{(j)}$ to ensure $\frac{1}{T}\widehat{F}^{(j)'} \widehat{F}^{(j)}=I$. • Stop after reaching certain criterion (say, $\frac{1}{\sqrt{N}}\|\widehat{B}^{(j)}-\widehat{B}^{(j-1)}\|\le \epsilon$ in which $\epsilon$ is a sufficiently small number).

Compared to BL2017, the above procedure is indeed more time-consuming practically. By utilizing the structure of the regressors, the CCE-type of approach is much appreciated for the computational efficiency. Thus, there is a trade-off between the flexibility of the model and the computational efficiency. In the literature, a few studies aim to justify the aforementioned algorithms theoretically. The early work probably dates back to Chen2014, and recently, Liu2020 and JYGH2020 further provide theoretical evidences for a semiparametric model and a parametric model respectively. In this paper, we do not pursue any theoretical results along this line of research, as it may lead to a different paper. Finally, in practice, one may follow Chen2014 to run the algorithm for several initial values and choose the solution that yields the highest value of the log-likelihood.

The Asymptotic Results

The following conditions are necessary before we present the asymptotic results in this subsection.

assumption• \begin{enumerate}[leftmargin=*] • Suppose that $\{\varepsilon_{it}\}$ is an array of identically distributed random variables in $(i,t)$ with known PDF and CDF as $g_\varepsilon(\cdot)$ and $G_\varepsilon(\cdot)$, respectively. • There exists a set $\Xi_{NT}=[\Xi_{NT}^l, \Xi_{NT}^u]$ such that all $z_{it}^0$'s belong to $\Xi_{NT}$ with probability approaching to 1, and $0<G_\varepsilon(\Xi_{NT}^l)<G_\varepsilon(\Xi_{NT}^u)<1$, where $z_{it}^0 =x_{it}'\beta_{0i} +\gamma_{0i}'f_{0t}$. • Let $\frac{1}{NT}\sum_{i,j=1}^N\sum_{t,s=1}^T E[ |E[e_{it}e_{js}\, | \, w_{it}^0, w_{js}^0]| ]=O(1)$, where $w_{it}^0=(x_{it}, \gamma_{0i}, f_{0t})$ and $e_{it} =\frac{1-y_{it}}{1- G_\varepsilon(z_{it}^0)}- \frac{y_{it}}{G_\varepsilon(z_{it}^0)}$. \end{enumerate}

Assumption (ref).1 requires the identical distribution, which is conventional in the literature on nonlinear models (e.g., Assumption 1 of Chen2020). Also, Remark (ref) partially explains why this condition is reasonable from the perspective of identification.

For Assumption (ref).2, the range of $\Xi_{NT}$ varies respect to the distribution considered. If a distribution is defined on $\mathbb{R}$ (say a normal distribution), we may allow the lower and upper bounds of $\Xi_{NT}$ to diverge to $\pm \infty$ respectively; if an exponential distribution is considered, we may let the lower bound converge to 0 and let the upper bound to diverge; etc. The design of $\Xi_{NT}$ is not new in the literature. For instance, both chenxh2015 and LTG2016 use a similar technique to accommodate some unbounded supports of the regressors.

The current form of Assumption (ref).3 allows for a certain type of weak cross-sectional dependence and time series correlation. Note that $E[e_{it}\, |\, w_{it}^0] =0$ by construction, so one can regard $e_{it}$ as a newly created residual term with mean 0. Assumption (ref).3 essentially imposes a restriction on the second moment of $e_{it}$, which can be verified for example for the independent and identically distributed (i.i.d.) case or for the case involving certain mixing conditions (e.g., Assumption (ref).3 below).

lemmaUnder Assumption (ref), as $(N,T)\to (\infty,\infty)$, \begin{eqnarray*} \frac{1}{NT}\sum_{i=1}^N\sum_{t=1}^T \left[G_\varepsilon(\widehat{z}_{it} ) - G_\varepsilon(z_{it}^0 )\right]^2 =O_P\left(\frac{1}{\sqrt{NT}}\right), \end{eqnarray*} where $z_{it}^0$ is defined in Assumption (ref), $\widehat{z}_{it} = x_{it}'\widehat{\beta}_i + \widehat{\gamma}_i'\widehat{f}_t$, and $\widehat{\beta}_i$, $\widehat{\gamma}_i$, and $\widehat{f}_t$ are defined in (ref).

With very limited restrictions, Lemma (ref) ensures the overall validity of the approach considered in this paper. It is noteworthy that Lemma (ref) still holds even for a model without regressors:

eqnarray[eqnarray omitted — 159 chars of source]

which is one of the models studied in Wang2019. Certainly, the maximum likelihood function should be adjusted in a very obvious manner. We conjecture that Lemma (ref) can help simplify the asymptotic development (and maybe assumptions) of Wang2019.

In addition, Lemma (ref) infers that a maximum likelihood estimation can reduce to a nonlinear least squares estimation more or less. We provide a few examples below.

Example 1: Interestingly, if $\varepsilon_{it}$ follows a uniform distribution, then the expression presented by Lemma (ref) completely possesses a form of the least squares approach:

eqnarray[eqnarray omitted — 205 chars of source]

As a result, most of the arguments made for the term $\widetilde{S}_{NT}(\beta, F)$ on pages 1264-1265 of Bai2009 will apply. We refer the interested readers to detailed discussions therein.

Example 2: If $\Xi_{NT}$ of Assumption (ref) is a compact set with fixed boundaries and $\inf_{z\in \Xi_{NT}} g_\varepsilon (z)\ge c_0 >0$, Lemma (ref) immediately yields that

eqnarray[eqnarray omitted — 264 chars of source]

in which the right hand side again reduces a term identical to $\widetilde{S}_{NT}(\beta, F)$ of Bai2009 ignoring the constant $c_0$. One in fact can allow $\Xi_{NT}=[\Xi_{NT}^l, \Xi_{NT}^u]$ to have diverging boundaries, and allow $g_\varepsilon(\cdot)$ to be a density function of normal distribution. If that is the case, we need an extra condition (i.e., Assumption (ref).1 below) to further regulate $\Xi_{NT}$ in order to achieve asymptotic consistency for all estimators of (ref). Simple algebra shows that under the condition $|\Xi_{NT}^l|+|\Xi_{NT}^u|= o_P(\sqrt{\log (NT)})$, Assumption (ref).1 is fulfilled, so the newly proposed model with the estimation procedure still holds.

remarkWe now discuss the selection of the number of factors before proceeding further. The topic has been studied in a variety of papers for the parametric linear models (e.g., BN2002, Onatski, LamYao2012, AhnHorenstein2013; and references therein). Among the results and discussions available in the literature, an important one is that many asymptotic results still hold true (but losing some efficiency), if the number of factors is over-specified when conducting the regression (e.g., FLM13, Moon). In this study, we find that such an argument is also valid for the binary response panel data model of (ref) under moderate restrictions. For simplicity, consider the Example 2 above. It is clear that the right hand side of (ref) completely reduces to a parametric model, so the arguments made for the parametric linear models apply immediately.

We are now ready to investigate the estimators of (ref), and emphasize again that we temporarily assume the number of factors is known, and work on the selection of the number of factors in Section (ref). The next assumption is essential in order to achieve consistency for each estimator of (ref).

assumption• \begin{enumerate}[leftmargin=*] • Suppose that there exists a sequence $\{a_{NT}\}$ satisfying that $\inf_{w\in\Xi_{NT}}g_\varepsilon(w)\ge a_{NT} >0$ and $a_{NT}\sqrt{NT}\to \infty$, in which $a_{NT}$ may converge to 0 or be constant. • Let $Z_i = X_i(\beta_{0i}-\beta_i ) $, and suppose that the following limits exist given $\| B-B_0\|/\sqrt{N}\le C$, where $C$ is a sufficiently large positive constant. \begin{enumerate} • $\sup_{B} \left|\frac{1}{NT}\sum_{i=1}^N (Z_i\otimes Z_i -E[Z_i\otimes Z_i]) \right| =o_P(1)$; • $\sup_{B} \left| \frac{1}{N\sqrt{T}}\sum_{i=1}^N( \gamma_i\otimes Z_i - E[\gamma_i\otimes Z_i]) \right| =o_P(1)$; • $\frac{1}{N}\Gamma_{0}'\Gamma_{0}\to_P\Sigma_\gamma$ and $\frac{1}{T}F_0'F_0\to_P\Sigma_f$. \end{enumerate} • Suppose that $0<\inf_{ F\in \mathbf{F} } \Omega (F )$, where $\Omega(F )=\operatorname*{\normalfont\textrm{diag}}\{\frac{1}{T}\Omega_{1T}(F ),\ldots, \frac{1}{T}\Omega_{NT}(F ) \}$, and \begin{eqnarray*} \Omega_{iT}(F) &=& E[X_i' M_F X_i\, | \, F ]-E[\gamma_{0i}\otimes (M_F X_i )\, | \, F ]' (\Sigma_\gamma\otimes I_T)^{-1}E[\gamma_{0i}\otimes (M_F X_i )\, | \, F ] . \end{eqnarray*} \end{enumerate}

Assumption (ref).1 bounds the PDF from below, which is similar to the treatments of Hansen2008 and LLL2012 among others. The two examples above should have clearly explained why such a condition is needed. The first two results of Assumption (ref).2 require some uniform consistency, which is due to the fact that the coefficients are indexed by $i$. For the homogeneous coefficient model, such conditions can be completely removed. These two conditions are similar to Assumption G of LCL2020, wherein heterogeneous coefficients are considered as well. Assumption (ref).3 is the identification restriction accounting for the heterogeneous setting of the model which is in the same spirit as Assumption D of Ando and the condition imposed on the term $Q_1(F_1)$ of HJPS. The reason why Assumption (ref).3 is necessary should be crystal clear in view of the two examples under Lemma (ref).

With Assumption (ref) in hand, we summarize the asymptotic consistency in the following lemma.

lemmaUnder Assumptions (ref) and (ref), as $(N,T)\to (\infty,\infty)$, \begin{enumerate} • $\frac{1}{N}\|\widehat{B}-B_0\|^2 =o_P(1)$, • $\frac{1}{NT} \|\widehat{F}\widehat{\Gamma}'-F_0\Gamma_{0}'\|^2 =o_P(1)$, • $\|P_{\widehat{F}} -P_{F_0} \| = o_P(1)$. \end{enumerate}

In Lemma (ref), the first two results guarantee the consistency of the estimators in (ref), while the third result infers that the space spanned by the columns of $F_0$ can be recovered consistently.

To carry on our analysis further, more structures and notations are needed:

eqnarray[eqnarray omitted — 707 chars of source]

where $u_{it}^0=(x_{it}',f_{0t}')'$. Moreover, let $\Omega_u=\operatorname*{\normalfont\textrm{diag}}\{\Sigma_{u,1},\cdots,\Sigma_{u,N}\}$, $\Omega_\gamma=\operatorname*{\normalfont\textrm{diag}}\{\Sigma_{\gamma,1},\cdots,\Sigma_{\gamma,T}\}$, and $\Omega_{u\gamma} =\{\Omega_{u\gamma, i t}\}_{N(d_\beta+d_f)\times Td_f}$. For the sake of space, we explain the necessity of these notations in Appendix (ref).

The following set of conditions are necessary to derive the rates of convergence.

assumption• \begin{enumerate}[leftmargin=*] • Let $ \max_{i\ge 1,t\ge 1}E\| x_{it}\|^{4+\delta}<\infty$, $\max_{i\ge 1}E\| \gamma_{0i}\|^{4+\delta}<\infty$, and $\max_{t\ge 1}E\| f_{0t}\|^{4+\delta}<\infty$, for some constant $\delta\geq2$. Also, suppose that $\max_{i\ge 1,t\ge 1}\| x_{it}\|=O_P(\log (N T))$, $\max_{i\ge 1}\| \gamma_{0i}\|=O_P(\log N)$ and $\max_{t\ge 1}\| f_{0t}\|=O_P(\log T)$. Let $F_0 \in \mathbf{F}$, and suppose that there exists $\delta^\ast\in(0,\delta)$ such that $\frac{T}{N^{1+\delta^\ast/4}}\rightarrow0$ and $\frac{N}{T^{1+\delta^\ast/4}}\rightarrow0$. • $g_\varepsilon(w)$ is twice differentiable on $\Xi_{NT}$, and $\sup_{w\in \Xi_{NT}} (|l^{(2)}_{it}(w)|+|l^{(3)}_{it}(w)|)<\infty$ uniformly in $(i,t)$, where $l^{(k)}_{it}(w)$ is the $k^{th}$ derivative of $l_{it}(w)$, Moreover, there exist $\rho_1,\rho_2<0$ such that \begin{eqnarray*} \sup_{w\in \Xi_{NT}} \rho_{\max}\left(\frac{1}{T}\sum_{t=1}^T l^{(2)}_{it}(w)u_{it}^0u_{it}^{0'}\right)\leq \rho_1, \, \sup_{w\in \Xi_{NT}} \rho_{\max}\left(\frac{1}{N}\sum_{i=1}^N l^{(2)}_{it}(w)\gamma_{0i}\gamma_{0i}^\top\right)\leq \rho_2, \end{eqnarray*} with probability one, where $\rho_{\max}(\cdot)$ has been defined in the last paragraph of Section (ref). • $\{\varepsilon_{it}\}$ is independent of $\{(x_{it}, \gamma_{0i}, f_{0t}): i\geq 1, \, t\geq 1\}$. Let $\{\varepsilon_{it}, x_{it},f_{0t}\}$ be strictly stationary and $\alpha$-mixing across $t$, and let $\alpha_{ij}(|t-s|)$ represent the $\alpha$-mixing coefficient. Moreover, assume that $\sum_{i,j=1}^N\sum_{t=1}^\infty(\alpha_{ij}(t))^{\delta/(4+\delta)}=O(N)$, $\sum_{i,j=1}^N (\alpha_{ij}(0))^{\delta/(4+\delta)}=O(N)$, and $\max_{i\geq 1} \sum_{t=1}^\infty (\alpha_{ii}(t))^{\delta/(4+\delta)}=O(1)$. • Suppose that $\frac{a^2_{NT}\sqrt{NT}}{[\log (NT)]^2}\rightarrow \infty$, $\rho_{\max}(\Omega_{u})<\infty$, $\rho_{\max}(\Omega_{\gamma})<\infty$, $\rho_{\max}\left(\frac{1}{NT}\Omega_{u\gamma}\Omega_{\gamma}^{-1}\Omega_{u\gamma}'\right)<\infty$, and $\rho_{\max}\left(\frac{1}{NT}\Omega_{u\gamma}'\Omega_{u}^{-1}\Omega_{u\gamma}\right)<\infty$. \end{enumerate}

Assumption (ref).1 imposes moments restrictions on $x_{it}$, $\gamma_{0i}$ and $f_{0t}$, which are standard in the literature. In addition, the rates $\log (NT)$, $\log N$ and $\log T$ can be easily fulfilled when measuring the maximum value of a series of random observations. See Assumption A7 of CHL2012 for example. We restrict the divergence of $N$ and $T$ by $\frac{T}{N^{1+\delta^\ast/4}}\rightarrow0$ and $\frac{N}{T^{1+\delta^\ast/4}}\rightarrow0$ which can be satisfied in many cases such as $N/T\rightarrow c$ where $c$ is a constant. This condition is used to establish the uniform convergence of the estimators (i.e., the first two results of Lemma (ref) below). The condition $F_0 \in \mathbf{F}$ is for the purpose of identification only. For instance, given $F_0$, we can always find a rotation matrix $W$ such that

eqnarray[eqnarray omitted — 47 chars of source]

As a consequence, we can write

eqnarray[eqnarray omitted — 69 chars of source]

which infers that instead of having $\gamma_{0i} $ and $f_{0t} $ as true parameters, we in fact use $ \gamma_{0i}'W^{-1}$ and $ Wf_{0t}$ as true parameters under Assumption (ref).1. For linear models (e.g., Bai2009 among others), a rotational matrix usually kicks in through the PCA procedure. However, for nonlinear models without closed-form estimators as in this paper, achieving an analytic form of such a rotation matrix seems to be impossible to the best of our knowledge. Similar issues can also be seen in AB2020, and AL2020.

Assumption (ref).2 slightly relaxes Assumption 2 of Wang2019 and Assumption 1.iv of Chen2020. The Probit and Logit models are obviously covered by this assumption.

Assumptions (ref).3 strengthens Assumption (ref).3 by imposing mixing conditions. This idea is consistent with Assumption A.2 of SC2013 and Assumption (ref).3 of FPSY2019. We can then establish the asymptotic distribution of $\widehat{\theta}_i$ from this set of low level conditions.

Assumption (ref).4 is used for the derivation involving the inverse of the Hessian matrix. It is mild since both $\Omega_\gamma$ and $\Omega_u$ are diagonal matrices. Using this assumption, we can show that the diagonal elements in the Hessian matrix play a dominating role when studying the asymptotic properties of the estimators.

We are now ready to present the rates of convergence in the next lemma.

lemmaUnder Assumptions (ref)-(ref), as $(N,T)\to (\infty,\infty)$, \begin{enumerate} • $\max_{i\geq 1}\|\widehat{\theta}_i-\theta_{0i}\|=o_P(1)$, • $\max_{t\geq 1}\|\widehat{f}_t-f_{0t}\|=o_P(1)$, • $\frac{1}{N} \|\widehat{\Theta} -\Theta_0 \|^2 = O_P\left(\frac{1}{N\wedge T}\right)$, • $\frac{1}{T}\|\widehat{F}-F_0 \|^2 =O_P\left(\frac{1}{N\wedge T}\right)$, \end{enumerate} where $\widehat{\Theta} = (\widehat{B},\widehat{\Gamma})$ and $\Theta_0 $ is defined under (ref).

In Lemma (ref), the first two results strengthen Lemma (ref), and show the consistency of $\widehat{\theta}_i$ and $\widehat{f}_t$ uniformly in $i$ and $t$. Based on these two results, we are able to establish the rates of convergence in the last two results of Lemma (ref), which lead to the derivation of the asymptotic distributions below.

To present the asymptotic distributions, the following conditions are necessary.

assumptionAssume that $\Sigma_{u,i}$, $\Sigma_{\theta,i}$, and $\Sigma_{\gamma,t}$ are positive definite matrices for $\forall i$ and $\forall t$, where $\Sigma_{\theta,i}=\lim_{T\rightarrow\infty} \frac{1}{T}\sum_{t=1}^T\sum_{s=1}^T E [g_\varepsilon(z^0_{it}) g_\varepsilon(z^0_{is}) e_{it}e_{is} u^0_{it}u_{is}^{0\prime}]$.

Assumption (ref) guarantees the positive definiteness of the asymptotic covariances. Under the mixing conditions of Assumption (ref).3, we can show that $\Sigma_{\theta,i}$ is the asymptotic covariance matrix of $\frac{1}{\sqrt{T}}\frac{\partial \log L(B_0, F_0,\Gamma_0)}{\partial \theta_i}$. Thus, we are able to establish the asymptotic distributions in the following theorem.

theoremLet Assumptions (ref)-(ref) hold, and $(N,T)\to (\infty,\infty)$. \begin{enumerate} • If $T(\log T)^2/N^2\rightarrow0$, then \begin{eqnarray*} \sqrt{T}(\widehat{\theta}_i -\theta_{0i}) \to_D N(0,\Sigma_{u,i}^{-1}\Sigma_{\theta,i}\Sigma_{u,i}^{-1}), \end{eqnarray*} where $\widehat{\theta}_i =(\widehat{\beta}_i',\widehat{\gamma}_i')'$ and $\theta_{0i} =(\beta_{0i}',\gamma_{0i}')'$. • If $N(\log N)^2/T^2\rightarrow0$ and $\frac{1}{\sqrt{N}}\frac{\partial \log L(B_0, F_0,\Gamma_0)}{\partial f_t}\to_D N(0,\Sigma_{f,t})$ for $\forall t$, then \begin{eqnarray*} \sqrt{N}(\widehat{f}_t -f_{0t}) \to_D N(0,\Sigma_{\gamma,t}^{-1}\Sigma_{f,t}\Sigma_{\gamma,t}^{-1}). \end{eqnarray*} \end{enumerate}

The conditions $T(\log T)^2/N^2\to 0$ and $N(\log N)^2/T^2\to 0$ in the body of the theorem are similar to those in Theorem 1 of BN2013. One may regard the above theorem as the binary response counterpart of Theorem 1 of BN2013.

If the error terms $\varepsilon_{it}$ are i.i.d., we have the following consistent estimators for the unknown matrices involved in Theorem (ref):

eqnarray[eqnarray omitted — 554 chars of source]

where $\widehat{u}_{it}=(x_{it}',\widehat{f}_t')'$, $\mathsf{g}_{it} (w) =\frac{[y_{it}-G_\varepsilon(w )] g_\varepsilon(w) }{[1-G_\varepsilon(w)]G_\varepsilon(w)}$, and $\mathfrak{g} (w)= \frac{[g_\varepsilon(w)]^2 }{[1-G_\varepsilon(w)]G_\varepsilon(w )} $. While there are weak serial correlation or cross-sectional dependence involved in $\varepsilon_{it}$'s, the constructions of $\widehat{\Sigma}_{\theta,i} $ and $\widehat{\Sigma}_{u,i}$ need to be adjusted accordingly, which in fact is a quite complex problem itself even just in the literature of time series analysis, see, Chapter 2 of FanYao, for example. Therefore, we do not further purse these estimators when weak correlation is involved over $(i,t)$.

Finally, we consider a mean group type of estimator for $\beta_{0i}$'s, which is often studied when heterogeneous coefficients are involved. Traditionally, the models with heterogeneous coefficients always assume that

eqnarray[eqnarray omitted — 54 chars of source]

where $\eta_i$ is i.i.d. over $i$ and has mean 0. One of the most cited works is Pesaran2006. As $\eta_i$ is only indexed by $i$, the estimate of $\beta_0$ is always achieved at a slow rate $\frac{1}{\sqrt{N}}$.

Interestingly, from the perspective of hypothesis testing, another strand of the literature considers the so-called “small departure" (e.g., Assumption 3 of goncalves_2011 and Eq. (3.5) of zhang2012inference), which can be formalized using the following assumption.

assumptionThere is a vector of unknown parameters-of-interest $\beta_0$ such that $\beta_{0i} = \beta_0 + \frac{1}{N^{\alpha}} \, \eta_i$, where $0\le \alpha\leq \frac{1}{2}$, and $\eta_i$ is a vector of independent and identically distributed (i.i.d.) random errors with $E[\eta_i]=0$ and ${\rm Var}[\eta_i] =\Sigma_{2}>0$. Moreover, $\{\eta_i\}$ is independent of $\{(x_{it}, \varepsilon_{it},\gamma_{0i}, f_{0t}): i\geq 1, \, t\geq 1\}$.

Assumption (ref) imposes a localized version on $\beta_{0i}$, and it narrows down the usual “departure" of the form: $\beta_{0i} = \beta_0 + \eta_i$ as commonly assumed in the relevant literature. If we do consider a testing problem such as $H_0: \beta_{0i} ={\beta}_0$, the choice of $\beta_{0i} = {\beta}_0 + \frac{1}{\sqrt{N}} \, \eta_i$ is naturally considered as a sequence of local alternatives with an optimal rate of convergence of an order $N^{-1/2}$ in the conventional parametric setting (see, for example, dong_gao_2018 for more details about testing small departures). In what follows, we aim to bridge both strands of the literature and explore under what condition an optimal rate $\frac{1}{\sqrt{NT}}$ can be achieved.

That said, we now consider a mean group type of estimator for $\beta_{0i}$'s using Assumption (ref). Define ${\widehat{\beta}}= \frac{1}{N} \sum_{i=1}^N \widehat{\beta}_i$ and $\overline{\beta}_{0} = \frac{1}{N} \sum_{i=1}^N \beta_{i0}$. Observe that

eqnarray[eqnarray omitted — 347 chars of source]

Equation ((ref)) indicates that the fast rate of convergence is achievable under Assumption 5. Depending on the limiting behaviour of $\rho_1 \equiv \frac{\sqrt{T}}{N^{1-\alpha}}$ and $\rho_2 \equiv \sqrt{\frac{T}{N}}$, the form of the asymptotic distribution of $\sqrt{T} \, N^{\alpha} ({\widehat{\beta}} - {\beta}_0 )$ may depend on both terms of ((ref)). In the case of $\rho_2 \rightarrow 0$, the first term contributes to the asymptotic distribution. When $\rho_1\rightarrow 0$, the second term mainly contributes to the asymptotic distribution. In the case where $\frac{\rho_1}{\rho_2} \rightarrow c\in (0, \infty)$, both terms contribute to the asymptotic distribution. In the relevant literature under the standard setting: $\beta_{i0} = {\beta}_0 + \eta_i$, however, the second term always dominates the asymptotic distribution. We now establish the following results in Theorem (ref).

theoremLet Assumptions (ref)-(ref) hold, and $(N,T)\to (\infty,\infty)$. \begin{enumerate} • Consider the case of $0\le \alpha <\frac{1}{2}$. If $\frac{N^{\alpha + \frac{1}{2}}}{T} \rightarrow 0$, then \begin{eqnarray*} N^{\alpha + \frac{1}{2}} ({\widehat{\beta}} - {\beta}_0 ) \rightarrow_D N(0, \Sigma_{2}). \end{eqnarray*} • Consider the case of $\alpha = \frac{1}{2}$. If $\frac{T}{N} \rightarrow \infty$, then \begin{eqnarray*} N ({\widehat{\beta}} - {\beta}_0 - {\rm bias}(N) ) \rightarrow_D N(0, \Sigma_{2}), \end{eqnarray*} where ${\rm bias}(N) = O_P\left(\frac{1}{N}\right)$. • Consider the case of $\alpha = \frac{1}{2}$. If $\frac{T}{N} \rightarrow \rho \in (0, \infty)$, then as $(N, T)\rightarrow (\infty, \infty)$ \begin{equation} \sqrt{N\, T} ({\widehat{\beta}} - {\beta}_0 - {\rm bias}(N, T) ) \rightarrow_D N(0, \Sigma_{12}), \end{equation} where $\Sigma_{12} = \Sigma_{1} + \rho^2 \, \Sigma_{2}>0$, \begin{eqnarray*} \Sigma_1 = \lim_{N, T} \frac{1}{NT} \sum_{i,j=1}^N \sum_{t,s=1}^TE[ \mathsf{g}_{it} (z_{it}^0) \mathsf{g}_{it} (z_{js}^0) \Sigma_{u, i}^{(d_\beta)} u_{it}^0 u_{js}^{0\prime} \Sigma_{u, j}^{(d_\beta)\prime}], \end{eqnarray*} $\mathsf{g}_{it}(\cdot)$ is defined under (ref), and $\Sigma^{(d_\beta)}_{u, i}$ includes the first $d_\beta$ rows of $\Sigma^{-1}_{u,i}$. \end{enumerate}

Theorem (ref) shows that we can establish the asymptotic distributions for the case of $\frac{T}{N} \rightarrow \rho \in (0, \infty]$. For the case of $\frac{T}{N} \rightarrow 0$, however, we have not been able to establish an asymptotic distribution for ${\widehat{\beta}}$. This is because we cannot improve the higher-order term $O_P\left(\frac{1}{N\wedge T}\right)$ involved in the leading-order approximation:

eqnarray*[eqnarray* omitted — 269 chars of source]

as shown at the beginning of the proof of Theorem (ref).

The first result of Theorem (ref) shows that the rate of convergence can be as fast as $N^{-(\alpha + \frac{1}{2})}$ for the case of $0\leq \alpha <\frac{1}{2}$ and $\frac{N^{\alpha + \frac{1}{2}}}{T} \rightarrow 0$. When $\alpha=0$, it reduces to the standard rate as established in the relevant literature (see, for example, Pesaran2006). The third result of Theorem (ref) shows that the conventional parametric rate of $(N T)^{-1/2}$ is achievable, but there is a bias term involved due to the “trade-off" between the asymptotic bias and the asymptotic variance. In Appendix (ref), we further discuss how to deal with biases. When $\{\varepsilon_{it}\}$ is i.i.d. over $i$ and $t$, we may be able to estimate $\Sigma_1$ by

eqnarray*[eqnarray* omitted — 218 chars of source]

The discussions made under (ref) still apply here.

We next move on to provide an information criterion to select the number of factors, and present a numerical implementation procedure.

Selection of the Number of Factors

In connection with Lemma (ref), we define the following information criterion:

eqnarray[eqnarray omitted — 236 chars of source]

where $\widehat{z}_{it}^{\mathsf{d}}$'s are obtained using (ref) by setting the number of factors as $\mathsf{d}$, $\xi_{NT}\to \infty$, and $ \frac{ \xi_{NT}}{\sqrt{NT}}\to 0$. We estimate $d_f$ by minimizing (ref):

eqnarray[eqnarray omitted — 160 chars of source]

where $d_{\max} (\geq d_f)$ is a user specified large fixed constant.

In general, the criterion always has a form as follows.

eqnarray[eqnarray omitted — 104 chars of source]

Such a form has been widely adopted by a wide range of selection procedures, e.g., traditional AIC/BIC, LASSO (HHM2008), and factor analysis (BN2002). The logic is that when under-selection occurs, a bias or a contradictory result will arise. As a consequence, the measurement of estimation errors in (ref) will be significantly large. When over-selection occurs, the estimation error will be asymptomatically the same as the situation of the correct-selection. However, due to the efficiency loss caused by the over-selection, the penalty term comes to play, and will yield a rate of convergence much larger than that associated with the estimator errors. That is why $\xi_{NT}$ needs to satisfy certain conditions specified under (ref). In the literature, a variety of log forms have been proposed for $\xi_{NT}$. Whether $\xi_{NT}$ has an optimal form remains unclear, but $\log\sqrt{N+T}$ works well for different data generating processes in our simulations below. Similar discussions have been provided to $g(N,T)$ of BN2002.

Last but not least, for the choice of the measurement of estimation errors, two well adopted forms are likelihood type of presentation (Eq. (2.9) of ChuZhuWang) or mean squared errors (Eq. (6) of HHM2008). In view of the development of Lemma (ref), the two forms are almost equivalent in our case. We adopt the latter for simplicity, and summarize the asymptotic property in the next theorem.

theoremUnder Assumptions (ref) and (ref), suppose further that $\xi_{NT}\to \infty$ and $ \frac{ \xi_{NT}}{\sqrt{NT}}\to 0$. As $(N,T)\to (\infty,\infty)$, $\Pr(\widehat{\mathsf{d}}=d_f)\to 1$.

It is worth mentioning that Theorem (ref) only requires very limited conditions, i.e., Assumptions (ref) and (ref). For each given value of the number of factors, we can implement the estimation procedure of Section (ref), which then allows us to calculate the information criterion (ref) to select the number of factors.

Simulation

In this section, we conduct simulations to examine the theoretical findings of Section (ref), and specifically consider the following data generating process.

eqnarray*[eqnarray* omitted — 162 chars of source]

The factors and loadings are generated by $f_{0t,j}\sim U(-2.5, 2.5)$ and $\gamma_{0i,j} \sim U(0, 6)$, where $\gamma_{0i,j}$ and $f_{0t,j}$ stand for the $j^{th}$ elements of $\gamma_{0i}$ and $f_{0t}$ respectively and $j=1,\ldots, d_f$. In order to introduce correlation between the regressors and the factor structure, we let $x_{it,j} = N(0,1)+0.5( |\gamma_{0i,1}|+|f_{0t,1}|)$, where $x_{it,j}$ stands for the $j^{th}$ element of $x_{it}$, and $j=1,\ldots, d_\beta$. For the coefficients, let $\beta_{0i,j} =i/N$, where $i=1,\ldots, N$, $j=1,\ldots, d_\beta$, and $\beta_{0i,j}$ stands for the $j^{th}$ element of $\beta_{0i}$. We consider two distributions for the error term, and for each distribution we vary the time series correlation and cross-sectional dependence of the error term in order to examine the sensitivity of the method proposed.

Case 1 -- light tailed distribution

itemize• DGP 1: $\varepsilon_{it} \sim N(0,1)$, where $N(0,1)$ is the standard normal distribution. DGP 2: Let $\varepsilon_t = \rho_\varepsilon \cdot \varepsilon_{t-1} +\Sigma_\nu^{1/2}\cdot\nu_t$, where $\rho_{\varepsilon} =0.3$, $\Sigma_\nu = \{0.3^{|i-j|}\}_{N\times N}$, $\varepsilon_t= (\varepsilon_{1t},\ldots, \varepsilon_{Nt})'$, and $\nu_t= (\nu_{1t},\ldots, \nu_{Nt})'$ with each $\nu_{it}$ being an independent draw from $N(0,1)$. DGP 3: Let $\rho_\varepsilon$ of DGP 2 be 0.7, and keep the rest settings identical to those of DGP 2.

Case 2 -- heavy tailed distribution

itemize• We still consider three DGPs as in Case 1, but replace all normal distributions with logistic distributions.

For each generated dataset, we first select the number of factors using the information criterion defined in (ref), and then conduct the estimation. We repeat the above procedure $M$ times, and report the following values to evaluate the finite sample performance:

eqnarray[eqnarray omitted — 586 chars of source]

where $\widehat{\mathsf{d}}_j $, $\widehat{B}_j$ and $\widehat{F}_j$ stand for the estimated number of factors, the estimated value of $B_0=(\beta_{01},\ldots, \beta_{0N})'$, and the estimated value of $F_0=(F_{01},\ldots, F_{0T})'$ at the $j^{th}$ iteration respectively; and $\beta_{0i}^{(\ell)}$ and $\widehat{\beta}_{i,j}^{(\ell)} $ stand for the $\ell^{th}$ element of $\beta_{0i}$ and its estimate at the $j^{th}$ replication.

We comment on these measures. $P_c$, $P_u$, and $P_o$ measure the probabilities of correctly, under, and over selecting the number of factors respectively. As explained in Section (ref), $\widehat{F}_j$ yields a consistent estimation of $F_0$ up to a rotation matrix, so we measure the distance between $P_{\widehat{F}_j}$ and $ P_{F_0}$ in (ref). In stead of looking at $B_0$ as a whole, $\text{Std}_{\beta_{0}^{(\ell)}} $ examines the stability of the estimation procedure. Specifically, the quantity $ \sqrt{ \frac{1}{M}\sum_{j=1}^M (\widehat{\beta}_{i,j}^{(\ell)} -\beta_{0i}^{(\ell)})^2 } $ provides the standard deviation associated with the estimates of $\beta_{0i}^{(\ell)}$. As $i$ runs from 1 to $N$, we further take the average over $i$.

We let\footnote{Due to the restrictions of computational power, we no longer explore larger values of $M$. In practice, heavier tails require longer time to compute, which should be expected.} $d_\beta=2$, $d_f=2$, $N,T\in \{50,100, 150\}$, and $M=500$. Moreover, let $\xi_{NT}$ of (ref) be $\log\sqrt{N+T}$ throughout the numerical studies without loss of generality. The results are summarized in Table (ref) and Table (ref) below. First, we point out that the newly proposed methodology works well regardless the tail behaviour of the error term $\varepsilon_{it}$, because the values associated with DGPs 1-3 are roughly same across both tables. Second, in Table (ref), the values of $P_c$ go up to 1, as the sample sizes increase. The pattern is same for DGPs 1-3. It is noteworthy that when the sample sizes are relatively small, the information criterion tends to under select the number of factors, which has a clear impact on the values of $\text{RMSE}_{F_0}$ presented in Table (ref). Third, in Table (ref), the values of $\text{RMSE}_{B_0}$, $\text{RMSE}_{F_0}$, and $ \text{Std}_{\beta_{0}^{(1)}} $ converge to 0 in general as the sample sizes go up, which should be expected. The exception is the case with $T=50$, in which the values of $ \text{Std}_{\beta_{0}^{(1)}} $ increase slightly as $N$ increases. Again, the pattern is same for DGPs 1-3. Fourth, compared to the Logit model, the Probit model tends to yield smaller values of Std$_{\beta_0}^{(\ell)}$, which is due to the fact that the Probit model has thin tails.

table[table omitted — 2,100 chars of source]
table[table omitted — 3,227 chars of source]

A Case Study

In this section, we apply the model and methodology to portfolio analysis using daily returns data of S&P 500 stocks. In both economics and finance disciplines, a vast amount of efforts have been devoted to portfolio analysis and management. Among them, a fundamental one is to select the optimal weights associated with each stock when constructing a portfolio. Mathematically, it is realized by the next minimisation problem.

eqnarray[eqnarray omitted — 105 chars of source]

where $1_{N}$ is a $N\times 1$ vector of ones, $\Sigma_p $ is a positive definite matrix and is usually decided by the data of stock returns, and $w$ includes the weights assigned to each stock. The analytic solution to the above minimisation problem is

eqnarray[eqnarray omitted — 93 chars of source]

More often than not, one adopts all stocks to construct $\Sigma_p$ (e.g., CLL2019, ELW2019), which naturally falls into a category of high dimensional matrix estimation. Thereby, to boost the estimation accuracy, we see the increasing popularity of methods using rank reduction (PX2019), penalization (CLL2019), or both (FLM13) among others. By default, the aforementioned techniques account for all stocks in practice, though some of them may have relatively small weights compared to the others. As pointed out in CD2006 and Nyberg2011, the sign of stock market returns may be predictable even if the returns themselves are not predictable. Also, CD2006 mention that “As volatility moves, so too does the probability of a positive return: the higher the volatility, the lower the probability of a positive return". In connection with the fact that a primary goal of portfolio analysis is to minimize the volatility (ELW2019), a binary response panel data model with interactive fixed effects naturally marries the above studies by modelling the probabilities of positive returns, so we can drop those having low probabilities.

Data

The stock prices data are collected from \url{https://www.kaggle.com} over the time period between 2 January 2008 and 31 December 2018. After removing the companies which have missing stock returns during the whole time period, we end up with $319$ stocks ($N= 319$). We adopt the log-normalised CBOE volatility index (VIX) (\url{http://www.cboe.com}) as the regressor, which is widely viewed as a good indicator of market sentiment (e.g., CD2006, PX2019). Furthermore, we collect risk-free interest (RFI) data from the U.S. Department of the Treasury (\url{https://www.treasury.gov}) to construct the Sharpe ratio later on in order to evaluate the performance of the proposed method.

Empirical Analysis

Below, we conduct a rolling-window analysis, and focus on the out-of-sample forecast. For each window, we estimate $\Sigma_p$ with the sample information on the most recent 505 trading days (roughly two years) using the next model.

eqnarray*[eqnarray* omitted — 164 chars of source]

where 1 and 0 stand for positive and non-positive returns respectively. We always count the first available trading day of each rolling-window as 0, and use $t=0, \ldots, 503$ to implement estimation. $x_{it}$ includes the value of VIX at day $t$, and $y_{i,t+1}$ is the sign of stock $i$ at day $t+1$. For each estimation, we first select the number of factors, and then conduct the estimation. This will provide us the following quantities: $\widehat{\beta}_i$'s, $\widehat{\gamma}_i$'s and $\widehat{f}_t$ with $t=1,\ldots, 503$. Afterwards, we bring the value of $x_{it}$ at $t=504$ in the estimated model to forecast the probabilities of the stocks having positive returns at $t=505$. Since $f_{0t}$ with $t=504$ is still unknown, we replace it by $\widehat{f}_t$ with $t=503$ as an approximation. The reason is that although $f_{0t}$ may vary over $t$, we do not expect any sudden jump at any particular time point. Alternatively, to forecast the next period of the unknown factors, one may follow BBE to impose more structure such as a VAR process. By doing so, one cannot select the number of factors in each estimation. Also, the number of optimal lags should be considered before one can conduct a more comprehensive investigation as in BBE. In this study, we on longer pursue the empirical results along this line in order not to deviate from our main goal.

We keep the stocks with estimated probabilities of positive returns greater than or equal to 0.5 to construct $\Sigma_p$ of (ref), and then calculate the weight vector $w$ using (ref). If all estimated probabilities are less than 0.5, then record $w=0$ (i.e., no transaction made for the day). Finally, we calculate the weighted average of returns for $t=505$. We repeat the above forecasting process from the first available window till the end, and consider the Probit model for the error term. The results of using other distributions (e.g., $t$-distribution and Logit model) for the error terms are quite similar, so we focus on the results of Probit model below. We follow ELW2019 to consider the problem of estimating the global minimum variance portfolio, in the absence of short-sales constraints.

As a comparison, when predicting the sign of we also consider a model with fixed effects only (referred to as FE below) and the model of BL2017 (referred to as BL below). In addition, we also consider two traditional approaches which utilize the entire stocks. Specifically, we consider the equal-weighted portfolio (referred to as EW below), which is a standard benchmark and has been promoted by DGU2007, among others. Also, we use correlation matrix for $\Sigma_p$ (referred to as CM below), which is mentioned in ELW2019. We acknowledge that many other methods are available when constructing the portfolio, such as the penalization method in FLM13, the nonparametric approaches in CLL2019 and PX2019, the dynamic covariance matrix estimation in ELW2019, etc. It is extremely hard to exhaust all possible methods in one study, as it may lead to a comprehensive review paper. Thus, we no longer compare with these methods in this article.

In what follows, we report (1). the average of weighted returns across the entire available period (Mean), which is annualized by multiplying 252; (2). the standard deviation of weighted returns across the entire available period (Std), which annualized by multiplying $\sqrt{252}$; (3). the information ratio defined as the ratio of Mean to Std (IR); and (4) the Sharp ratio defined as the mean of returns minus the risk free interest normalised by Std (SR). The annualized Mean and Std are consistent with those defined in Section 6.2 of ELW2019. We refer interested readers to their paper for more relevant discussions. Moreover, since a strand of literature on portfolio analysis is interested in estimation from a factor model with a predetermined number of factors (PX2019), we report the results for the cases where the number of factors are fixed as $1,\ldots, 5$ respectively. We also report the results of the case where IC of (ref) is used to select the optimal number of factors for each estimation. The results are summarized in Table (ref).

table[table omitted — 1,000 chars of source]

In general, we are looking for a strategy, which generates a small value of Std, but large values of Mean, IR and SR. The CM approach yields the largest Mean, but also has the largest Std, which should be least preferred due to the high risk caused by the high volatility. The FE model yields the smallest Std, which in a sense should be expected due to simplicity of the model. It also produces the minimum values of Mean, IR, and SR respectively. The IFE models (with the “Optimal" number of factors and one factor respectively) outperform the rest models regarding the criteria IR and SR, which are commonly adopted to measure the performance of portfolios (e.g., PX2019, ELW2019). Overall, the newly proposed framework has a reasonably good performance.

Finally, we acknowledge the limit of the current empirical study. For example, one may adopt a VAR structure for the unobservable factors as in BBE, and investigate further the impulse responses and the optimal number of lags. Also, one may further conduct the penalization estimation for the stocks having high probability of getting positive returns, which bridges the literature of binary response models (e.g., Chen2020, Wang2019) and the literature of high-dimensional covariance matrix estimation (e.g., FLM13). In order not to deviate from our main goal, we do not pursue these results in the current study.

Conclusion

In this paper, we investigate binary response models for heterogeneous panel data with interactive fixed effects by allowing both the cross-sectional dimension and the time dimension to diverge. From the modelling perspective, our setting is similar to BL2017, but we do not require a specific structure on the regressors, which allows us to avoid putting any restriction between the number of regressors and the number of unobservable factors. Our investigation establishes a link between a maximum likelihood estimation and a least squares approach. As a consequence, the identification restrictions provided in Bai2009 and Moon are readily to be applied to binary response models with very minor modifications. We further establish asymptotic distributions for the unobservable factors and their loadings, which can be considered as the binary response counterpart of those established in BN2013. In addition, we provide a simple information criterion to detect the number of factors. Last but not least, we conduct intensive numerical studies to examine the finite sample performance of the newly proposed model and methodology, and demonstrate the practical relevance.

From a practical perspective, the framework can be applied to predict the probability of corporate failure (CCL2014), conduct credit rating analysis (JJW2015), etc. Following CD2006 and Nyberg2011, in the empirical study, we focus on the sign prediction of stock returns, and then use the results of sign forecast to to conduct portfolio analysis. By implementing rolling-window out of sample forecasts, we demonstrate the practical relevance of the paper.

In the future work, it might be interesting to consider a network model such as those considered in Yanetal2019 and Dzemski. We conjuncture that a result like Lemme (ref) might be achievable to simplify the asymptotic development. Moreover, a high dimensional model with sparse coefficients such as that in ChuZhuWang is also worth to be investigated.

Acknowledgements

Gao acknowledges financial support from the Australian Research Council Discovery Grants Program under Grant Numbers: DP170104421 and DP200102769. Peng acknowledges the Australian Research Council Discovery Grants Program for its financial support under Grant Number DP210100476.

{