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.
90,722 characters · 16 sections · 55 citation commands
Linear Multidimensional Regression with Interactive Fixed-Effects
Models of multidimensional data -- data with more than two dimensions -- are becoming increasingly important in econometric analysis as large data sets with a multidimensional structure become available. For example, consider studying demand elasticities with consumption data that varies by product, $i$, store, $j$, repeated over time, $t$.\footnote{A non-exhaustive list of related examples can be found in the introduction of matyas2017econometrics. } In this context, analysts may be concerned with shifts in taste preferences unobserved by the econometrician related to unobserved characteristics in each dimension - like a cultural or sporting event that impacts prices and demand heterogeneously over product, store, and time. Whilst these are not explicitly observed by the econometrician, fixed-effects are a useful tool to infer unobserved heterogeneity and control for such variation.
Additive fixed-effects in multidimensional data can at most accommodate variation in unobserved heterogeneity over a subset of dimensions with any of the fixed-effects terms. For example, in the three-dimensional model, additive effects only control for variation over $ij$, $it$ and $jt$, but not jointly over all $ijt$. Hence, if heterogeneity is across all dimensions interactively, additive fixed-effects cannot sufficiently control for unobserved heterogeneity. This paper develops tools to control for unobserved heterogeneity in the form of interactive fixed-effects that controls for variation that interacts over all dimensions of the data.
Consider $\beta$ estimation in the interactive fixed-effects model with three dimensions:
where all terms in $\sum_{\ell = 1}^L \lambda_{i\ell} \delta_{j\ell} \gamma_{t\ell}$ are unobserved.\footnote{Bounded $L$ is only needed for terms interacting over all dimensions, and for the parametric rate of consistency. Less sparse models with $L$ unbounded are possible, but lead to slower convergence. } Considering three dimensions is without loss of generality for the methods used herein. The object $L$ is often difficult to calculate, see haastad1989tensor. However, for purposes of this analysis it is only necessary to know what is called the multilinear rank, defined in Section (ref). The multilinear rank is the set of matrix-ranks of the different matrices the array can be transformed into along each indice. For example, the multilinear rank of the interactive fixed-effects in three dimensions is a vector of three potentially different ranks - the first being the rank when the first dimensions are the rows, and the second and third dimensions are jointly the columns, and so on for the other multilinear rank entries. In the presence of covariates, estimation of interactive fixed-effect rank is still an open area of research in the standard panel setting. Simulations suggest estimation with more than the required components performs well.\footnote{In standard panel data overestimating rank leads to valid inference, see MoonWeidner2015. } Additive fixed-effects are omitted for brevity but are subsumed by the interactive fixed-effect term or can be removed with a simple within transformation.
Let $X_{ijt}$ be arbitrarily correlated with interactive fixed-effects, $\sum_{\ell = 1}^L \lambda_{i\ell} \delta_{j\ell} \gamma_{t\ell}$, but independent of the noise term, $\varepsilon_{ijt}$. The challenge to estimating $\beta$ is isolating variation in $X_{ijt}$ that is not correlated with the interactive fixed-effects term. This paper develops multidimensional weighted-within transformations that project out this unobserved heterogeneity and shows settings where standard factor methods work well. The weighted within transformation is a novel contribution. Group fixed-effects from bonhomme2022discretizing, and in freeman2023linear are also possible.
This paper makes two main contributions to the literature. The first is to show that the three or higher dimensional model can be couched in a standard two-dimensional panel data model, and to derive sufficient conditions for consistency using factor model methods from Bai2009,MoonWeidner2015, albeit at slow rates of convergence. The second contribution is to introduce weighted fixed-effects methods, which, when combined with a double debias procedure, achieves the parametric rate of convergence and asymptotic normality for $\beta$ estimates. This is yet to be shown in the three dimensional case using existing methods. Simulations corroborate these theoretical findings and an empirical demand estimation application demonstrates the estimator in practice.
The novel estimator proposed in this paper can be described as an extension to the usual within transformation. Consider additive fixed-effects of the form $a_{ij} + b_{it} + c_{jt}$. These can be projected out using the transformation,
applied equivalently to $X_{ijt}$, where the bar variables denote the average taken over the “dotted” index for the entire sample. That is, $\Bar{{Y}}_{\cdot jt} := \frac{1}{N_1}\sum_{i = 1}^{N_1} Y_{ijt}$, $\Bar{{Y}}_{\cdot\cdot t} := \allowbreak \frac{1}{N_1N_2} \sum_{i = 1}^{N_1}\sum_{j = 1}^{N_2}\allowbreak Y_{ijt}$, etc. In the presence of interactive fixed-effects, the simple within transformation demeans each fixed-effect term to leave,
which is clearly not controlled for if the $\{\lambda,\delta,\gamma\}$ terms admit heterogeneity over $i,j,t$.
Consider instead an extension of the within transformation that uses weighted means instead of uniform means. With an abuse of notation, consider,
where the bar variables combined with the star indices denote weighted means for that observation. For example, $\Bar{Y}_{i^*jt} = \sum_{i^\prime}w_{i,i^\prime}Y_{i^\prime jt}$ for weights across $i^\prime$ for each $i$. In turn, the remainder from the interactive fixed-effects term is,
If appropriate weights are used, the weighted within transformation can project out the interactive fixed-effects. The same can be true if the weighted means are replaced with cluster means for appropriate cluster assignments, similar to a group fixed-effects estimator. Hence, with a relatively small change to how the within transformation is performed, much more general fixed-effects can be considered in the linear model.
The model for interactive fixed-effects has precedent in the standard two-dimensional panel data setting, for example in Bai2009,Pesaran2006,
The interactive term $\sum_{\ell = 1}^L\lambda_{i\ell} f_{t\ell}$ also sufficiently captures variation in additive individual and time effects without the need to specify these separately. For multidimensional applications the problem in (ref) can be transformed to a two dimensional problem and estimated as (ref) directly using the transformed data. Indeed, Section (ref) displays useful preliminary convergence rates for this approach. However, convergence rates can suffer severely from the over-parametisation implied by (ref). Without strong assumptions on the sparsity of the fixed-effects, only a slow rate of convergence can be guaranteed for this approach. However, when this approach is used to construct preliminary estimators of the fixed-effects, the parametric rate of consistency is possible when these are combined with the novel weighted within transformations. kuersteiner2020dynamic consider a similar dimension specific projection in the time dimension in the presence of error correlations, in particular see Supplement D.5.2.
On top of this slow convergence rate for the two dimensional transformation, finite sample bias may also arise when $L$ is large and only a subset of the unobserved heterogeneity parameters are low-dimensional. For unbiased estimates of $\beta$, transforming the multidimensional array to a matrix then estimating (ref) requires either: (a) all the fixed-effects are low-dimensional; or (b) that a known subset of the fixed effects are low-dimensional, which is a strong assumption. Alternatively, whilst the weighted transformations also require that a subset of the fixed-effect parameters are low-dimensional, the analyst does not need to know which ones are. In this sense, the novel weighted estimator is robust. A concrete example is considered in a simulation exercise, and there is evidence in the empirical application of differences in the rank of fixed-effects in each dimension.
The beer demand elasticity application uses Dominick's supermarket data for Chicago, 1991-1995, where price and quantity vary over product, store, and fortnight. For comparison, barley price is used as an instrument to estimate elasticities. IV estimates show demand is strongly negatively elastic (-3.39), however, estimates are very imprecise. This could be because the instrument only varies over one dimension, time, or, e.g., because prices of other inputs are also highly variable so the instrument does not explain much price variation. Estimates from the weighted within transformation also show demand is downward sloping and elastic (-3.12), but with much better precision. Factor model estimates are highly sensitive to how the data is transformed into two-dimensions. The estimates from the weighted-within transformation are similar to the own-price elasticities estimated in hausman1994competitive.
The technical component of this paper relates to the numerical analysis literature on low-rank approximations of multidimensional arrays. As pointed out in de2008tensor, the low-rank approximation problem in the tensor setting is not well-posed, hence presents technical difficulties. See kolda2009tensor for a summary of the multidimensional array decomposition problem, and vannieuwenhoven2012new, rabanser2017introduction. Therefore, it is necessary to innovate on this tensor low-rank problem to find appropriate analytical results. To this end, this paper utilises well-posed components of two-dimensional sub-problems for use in nuisance parameter applications. The advantage in this paper is that fixed-effects only need to be differenced out at a sufficiently fast asymptotic rate, and the low-rank tensor problem does not need to be solved explicitly.\footnote{elden2011perturbation, along with related papers, suggest a reformulation of the low multilinear rank problem that may have promising applications. }
The paper is organised as follows: Section (ref) motivates the estimator with a beer demand estimation empirical application; Section (ref) introduces the general model; Section (ref) details the estimators and convergence results; Section (ref) discusses the convergence results and some practical considerations; Section (ref) displays the simulation results; and Section (ref) concludes.
The proposed methods are applied to estimate the demand elasticity for beer. Price and quantity for beer sales is taken from the Dominick's supermarket dataset for the years 1991-1995 and is related to supermarkets across the Chicago area. Price and quantity vary over three dimensions in this example -- product ($i$), store ($j$) and fortnights ($t$). Fixed-effects that interact across all three dimensions can control for temporal taste shocks to beer consumption that differ over both product and store. Take for instance a large sporting event (temporary $t$ shock) that changes preferences differently across locations ($j$) and across certain subsets of sponsored beer ($i$). Analysts may want some way to control for such heterogeneity that interacts over all dimensions, even if they do not explicitly observe such sponsorship variables. For example, in the stadiums for the many NBA finals playoffs the Chicago Bulls played in the early 1990's, Miller Lite beer advertisements could be seen alongside advertisements for a substitute product, Canadian Club whisky. This suggests these events attracted large marketing campaign spends for these and other beer substitute brands that most likely also included price offers at local supermarkets. Whilst the impact of these advertisements and price offers on the demand for or price of beer is not clear and, further, that it is reasonably safe to assume the econometrician does not observe the plethora of marketing campaigns around these events, the analyst would most likely still want to control for aggregate shocks like these.
In their critique of the fixed-effects approach to demand estimation, berry2021foundations note the approach in its simplistic additive form cannot control effects that interact over products $\times$ markets, see Section 2.5.1 in berry2021foundations. This limitation is significantly relaxed with interactive terms over all three dimensions. Indeed, Table (ref) suggests even products $\times$ markets interactions are limiting - considering products $\times$ markets $\times$ time substantially shifts point estimates, cf. “Additive Fixed-effects” versus “Weighted-within” estimates.
Models for demand estimation ideally account for endogenous variation in prices and quantity. The classic instrumental variable approach is to find a supply shifter that shifts the supply curve, allowing the econometrician to trace out the slope of the demand curve. A popular instrument in the estimation of beer demand is the commodity price for barley, one of the product's main ingredients, see e.g. saleh2014simple,tremblay1995advertising, richards2021dynamic. Since the price of barley is arguably not driven by the demand for it by any one supplier of beer, it can be a useful variable to instrument for price shifts. In the following, it is taken as given that the price of barley is exogenous with respect to the fixed-effects and noise term, $\varepsilon$.
For valid inference, the instrument should be highly correlated with the price of beer. In this dataset, correlation between the price of barley, which varies over only month, every second $t$, and price of beer depends on how beer price is first aggregated. If beer price is first integrated over $i$, $j$, and to the monthly level, such that it only varies over every second $t$, then it is highly correlated with the price of barley, at 0.79. However, if beer price is not aggregated at all it is only correlated at 0.001. This heuristic suggests there are important product and store level price drivers for beer that are not accounted for by fluctuations in the price of barley, or indeed by price fluctuations only over time. Standard errors for the IV estimator below are much larger than other estimators, which may be explained partly by this loss in effective sample size.\footnote{The effective sample size for the IV estimator drops from $N_1N_2T$ to just $T$. }
Table (ref) refers to estimates from the following model,
No additional controls are included here since they are low-dimensional and subsumed by additive fixed-effects. Estimates for pooled OLS are positive, a contradiction that demand curves are downward sloping. IV estimates are negative, but standard errors are large. As noted above, the effective sample size for the second stage of the instrumental variable approach is orders of magnitude smaller because the instrument only varies over time, specifically every second fortnight. Estimates under the additive model for fixed-effects are also negative, with much smaller standard errors than the IV estimate.
The matrix method is implemented in all three dimensions, i.e., transform the three-dimensional array into three different matrices - with products as rows, with stores as rows, and time as rows. Five factors are estimated for each. The results suggest the factor model for fixed-effects is very sensitive to how the array is organised into a two-dimension problem, with vastly different conclusions drawn from each specification. Simulations in Section (ref) suggest an explanation - sparsity of the fixed-effects in the interaction term can be heterogeneous, and imply a different factor model rank when the tensor of data is flattened, or matricised, in different dimensions.
This paper's novel estimator, the weighted-within, estimates elasticities similar to the IV point estimates, but with much smaller standard errors.\footnote{ To choose bandwidths: For $h \leq 0.2$, coefficient and variance estimates become increasingly unstable, thus this was the lower bound of consideration. Estimates did not change in an economically meaningful way with $h$, so the bandwidth was chosen to minimise standard errors. } Weighted within estimates are similar to the own-price elasticity from Table 1 in hausman1994competitive, which average around $-2.5$.
Figure (ref) displays estimates and confidence bands for the factor model with products as rows, and the weighted-within estimator, allowing the number of interactive terms to increase. The figure displays two things. First, there is evidence that up to approximately 7-8 interactive terms may be appropriate. Second, there is a substantial shift in estimate with the weighted-within estimator across the number of interactive terms. This implies a persistent debias with the weighted-within transformation. HAC confidence intervals for the weighted-within estimator are also substantially tighter than the factor model.
Let $\beta^0$ denote the true parameter value for the slope coefficients. The model in full dimensional generality is, \footnote{ For example, in index notation this model can be written as,
with $\mathcal{A}_{i_1, i_2, \dots, i_d} = \sum_{\ell = 1}^L \varphi^{(1)}_{i_1\ell} \dots \varphi^{(d)}_{i_d\ell}$. }
where $\boldsymbol{Y}, \boldsymbol{X}_k, \boldsymbol{\varepsilon} \in\mathbb{R}^{N_1\times N_2\times\dots \times N_d}$. $\boldsymbol{\mathcal{A}} = \sum_{\ell = 1}^L \varphi^{(1)}_\ell \circ \dots\circ \varphi^{(d)}_\ell$ where $\varphi^{(n)}_\ell\in\mathbb{R}^{N_n}$ for each $n = 1, \dots, d$ and “$\circ$” is the outer product. $L$ is bounded and fixed. $\boldsymbol{\varepsilon}$ is a noise term independent of all $\boldsymbol{X}_k$ and all unobserved fixed-effects terms. Take $i_n\in\{1,\dots,N_n\}$ for all $n\in\{1,\dots,d\}$ as the dimension specific index, where $N_n$ is the sample size of dimension $n$. The regressors $\boldsymbol{X}_k$ may be arbitrarily correlated with $\boldsymbol{\mathcal{A}}$. Throughout this paper all dimensions are considered to grow asymptotically, that is $N_n \rightarrow \infty$ for all $n$, however, for the asymptotic theory only a subset of two dimensions need to grow asymptotically.
Model (ref) can be seen as a natural extension of the Bai2009 model to three (or more) dimensions with $\boldsymbol{\mathcal{A}}$ interpreted as a “higher-dimensional” factor stucture. Similar to this strand of the literature, all terms in $\boldsymbol{\mathcal{A}}$ are considered fixed nuisance parameters. There are potentially many extensions to the factor model setting in Bai2009 to the higher dimension case. This paper starts with what seems the most natural extension.
$\mathcal{A}$ incorporates additive fixed effects that vary in any strict subset of dimensions. For example, in the three dimensional setting it can control for the additive terms, $a_{ij} + b_{it} + c_{jt}$. These can be controlled for using $L = \min\{N_1, N_2\} + \min\{N_1, N_3\} + \min\{N_2, N_3\}$, with the first $\min\{N_1, N_2\}$ terms $\sum_{\ell = 1}^{\min\{N_1, N_2\}} \varphi_{i\ell}^{(1)}\varphi_{j\ell}^{(2)} = a_{ij}$ by making $\varphi_{t\ell}^{(3)}$ constant for $\ell = 1,\dots, \min\{N_1, N_2\}$, and so on for the $ b_{it}$ and $c_{jt}$.\footnote{Strictly speaking, this is not allowed if $L$ is bounded and fixed. Bounded and fixed $L$ is really only required for any fixed-effect term that varies over all $i,j,t$. } These could also be controlled for directly using the standard within-transformation before considering the model in (ref).
For a $d$-order tensor, $\mathbf{A}$, a factor-$n$ flattening, denoted as $\mathbf{A}_{(n)}$, is the rearrangement of the tensor into a matrix with dimension $n$ varying along the rows and the remaining dimensions simultaneously varying over the columns. That is, $\mathbf{A}_{(n)}\in\mathbb{R}^{N_n\times N_{n+1}N_{n+2}\dots N_1\dots N_{n-1}}$. The Frobenius norm, $\|\cdot\|_F$, of a matrix or tensor is the entry-wise norm, $\|\mathbf{A}\|_F^2 = \sum_{i_1 = 1}^{N_1} \dots \sum_{i_d = 1}^{N_d} A_{i_1\dots i_d}^2$. The spectral norm, denoted $\| \cdot\|$, is the largest singular value of a matrix. For a $d$-order tensor, $\mathbf{A}$, the multilinear rank, denoted $\mathbf{r}$, is a vector of matrix ranks after factor-$n$ flattening in each dimension, with each component of this vector $r_n = rank\big(\mathbf{A}_{(n)}\big)$. Tensor rank, different to multilinear rank, is defined as the least number of outer products of vectors to replicate the tensor. That is, for tensor $\mathbf{A}$ and vectors $u_\ell^{(n)}\in\mathbb{R}^{N_n}$, tensor rank is the smallest $L$ such that $\mathbf{A} = \sum_{\ell = 1}^L u_\ell^{(1)}\circ \dots\circ u_\ell^{(d)}$, where $\circ$ is the outer product of a vector. The notation $a\lesssim b$ means the asymptotic order of $a$ is bounded by the asymptotic order of $b$. The $n$-mode product of a tensor $\mathbf{A}$ and matrix $B$ is denoted $\mathbf{A}\times_n B$ and has elements $$(\mathbf{A}\times_n B)_{i_1,\dots,j,\dots,i_d} = \sum_{i_n = 1}^{N_n} A_{i_1,\dots, i_n, \dots, i_d} B_{j,i_n}, $$ which is equivalent to saying the flattening $(\mathbf{A}\times_n B)_{(n)} = B\mathbf{A}_{(n)}$.
The singular value decomposition of a matrix, ${A}\in \mathbb{R}^{N_1\times N_2}$ is
where $U$ is the matrix of left singular vectors, $u_r$, $V$ is the matrix of right singular vectors, $v_r$, and $\Sigma$ is a diagonal matrix of singular values, $\sigma_r$, with values running in descending order down the diagonal. For a rank-$r$ matrix, the first $r$ entries on the diagonal of $\Sigma$ are strictly positive and the remaining entries are zero.
Take the approximation problem,
It is well known from the Eckart-Young-Mirsky theorem that the solution to this approximation problem is the first $k$ terms of the singular value decomposition, i.e. $\sum_{r = 1}^{k} \sigma_r u_r v_r^\prime$. The Eckart-Young-Mirsky theorem effectively picks out the row and column subspaces that best explain variation in the matrix $A$ as the leading columns of the matrix $U$, respectively of $V$. The sum of squared error at the minimiser is thus $ \sum_{r = k + 1}^{\min\{N_1,N_2\}} \sigma_r^2$. This is commonly called a low-rank approximation and forms the cornerstone for estimation of unobserved heterogeneity in the factor model and interactive fixed-effects models in BaiNg2002,Bai2009, MoonWeidner2015 amongst others.
The Eckart-Young-Mirsky theorem, however, does not extend to the three or higher dimensional setting, see de2008tensor for details. To avoid this complication, the multidimensional problem can be translated to the two-dimensional setting to utilise the Eckart-Young-Mirsky theorem, or the fixed-effects parameters need to be shrunk separately, as is done with the weighted-within transformation.
This section starts with a brief explanation of the main estimation approach to be used for inference. This main estimator draws on the double debias approach summarised in chernozhukov2022locally. Sections (ref)-(ref) detail two preliminary estimation steps required for the main estimator, and Section (ref) details how these fit together in the main estimator to produce an asymptotically normal estimator centred at the true value. The two preliminary estimators may be of independent interest. The second preliminary estimator in Section (ref), and main estimator in Section (ref) are novel to this paper. The first preliminary estimator in Section (ref) is also novel, but is a very simple adaptation of two-dimensional estimators.
The main estimator used for inference is presented first as an infeasible estimator, with the preliminary steps making it feasible. Consider the conditional expectation, $\mathbb{E}\left(\boldsymbol{Y}|\boldsymbol{\mathcal{A}}\right) = \Gamma_{Y}$, $\mathbb{E}\left(\boldsymbol{X}|\boldsymbol{\mathcal{A}}\right) = \Gamma_{X}$, such that,
with $\mathbb{E}\left(\boldsymbol\eta|\boldsymbol{\mathcal{A}}\right) = 0$, $\mathbb{E}\left(\boldsymbol{\varepsilon}|\boldsymbol{X},\boldsymbol{\mathcal{A}}\right) = 0$. The display for $\boldsymbol{X}$ is general, in that if $\boldsymbol{X}$ is unrelated to $\boldsymbol{\mathcal{A}}$, $\boldsymbol{\Gamma}_{X} = 0$. Hence, the display in (ref) is a representation, not an imposed model. If the correlation between $\boldsymbol{X}$ and $\boldsymbol{\mathcal{A}}$ is weak in the sense that the empirical mean of $\boldsymbol{\Gamma}_{X_k}^2$ converges to zero at arbitrary rate, then the following orthogonality conditions are not required and the convergence result from the preliminary estimator is sufficient. In this sense, the estimation approach below is robust to the existence of $\boldsymbol{\Gamma}_{X_k}$, and how prevalent it is in the generating process for each $X_k$ up to weak regularity conditions.
The moment conditions used to estimate $\beta$ are,
These conditions are Neyman Orthogonal with respect to $\Gamma_Y,\Gamma_X$ in that,
This property allows error to be multiplicative over estimation error in $\Gamma_Y,\Gamma_X$ in the feasible version of the estimator, since only second-order terms appear in the error. The weighted-within transformation applied to $Y$ and $X$ in the main estimator is automatically endowed with this Neyman orthogonal property, since it is simply a linear smoother across both $Y$ and $X$. That is, for linear smoother $S$, $SY = SX\beta^0 + S\mathcal{A} + S\varepsilon =\hat \Gamma_X \beta^0 + \hat{\mathcal{A}} + S\varepsilon = \hat\Gamma_Y$.
The (infeasible) inference corrected estimator, $\widehat\beta_{IC}^{(infeasible)}$, is
where ${\rm vec}_K(\boldsymbol{X} - \boldsymbol{\Gamma}_{X})\in \mathbb{R}^{\prod_n N_n \times K}$ is shorthand for the matrix of vectorised covariates, with each column the vectorised $(\boldsymbol{X}_k - \boldsymbol{\Gamma}_{X_k})$. The ${\rm vec}$ notation for $(\boldsymbol{Y} - \boldsymbol{\Gamma}_{Y})$ is the standard vectorisation. This estimator is robust to fixed-effects appearing in the equation for $Y$ and/or $X$, and is less sensitive to misspecification of the fixed-effects in each. Terms $\Gamma_Y$ and $\Gamma_X$ do not need to be explicitly estimated for the feasible estimator, only differenced out from $Y$, respectively $X$.
The second preliminary estimator in Section (ref) below, in particular the projections involved there, provide these objects. The first preliminary estimator in Section (ref) is required for the second preliminary estimator. As is shown in Section (ref), this estimator is robust to fixed-effects appearing in very general form in $X$. Indeed, slow consistency rates for ${\boldsymbol{\Gamma}}_{X}$ estimates are sufficient under weak regularity in $\varepsilon$.
This section provides a description of some matrix methods that can be applied directly to the multidimensional model and stipulates the assumptions required for consistency. kapetanios2021estimation employ a similar approach for three-dimensional arrays in conjunction with the Pesaran2006 common correlated effects estimator. babii2022tensor employ a similar matricisation procedure as that detailed below, but are interested in inference on the fixed-effect parameters.
Consider recasting the multidimensional array problem into a two dimensional panel problem by flattening $Y$ and $X$ in the $n$-th dimension,
where $Y_{(n)}, X_{(n)} ,\varepsilon_{(n)} \in \mathbb{R}^{N_n\times \prod_{n^\prime \neq n}N_{n^\prime}}$, $\varphi^{(n)}$ is an $N_n\times r_{n}$ matrix and $\Gamma_n$ is an $\prod_{n^\prime \neq n}N_{n^\prime}\times r_{n}$ matrix that accounts for variation in the remaining $\varphi^{(n^\prime)}$ for all $n^\prime\neq n$. The term $r_{n}$ is indexed by the dimension $n$ because it may vary non-trivially according to the flattened dimension. This is then the model described in (ref), that is, the standard linear model with factor structure unobserved heterogeneity as studied in Bai2009.
The two-dimensional fixed-effects estimator for a given flattening, $n$, optimises
Then $\widehat{\beta}^{2D}_{(n)} = \operatorname*{argmin}_{\beta} R(\beta, \widehat{r}_n, n)$ is the slope estimate for the two-dimensional setup. The analyst must choose both the dimension to flatten in, $n$, and the rank of the estimated interactive fixed-effects term, $\widehat{r}_n$. It is well known that the minimum in (ref) is achieved using the leading $\widehat{r}_n$ terms from the singular value decomposition of the error term, $Y_{(n)} - X_{(n)}^\prime \beta$. This gives $\widehat{\varphi}^{(n)}$ as the first ${\widehat{r}_n}$ columns of $\widehat{U}\widehat{\Sigma}$ and $\widehat{\Gamma}_n$ as the first ${\widehat{r}_n}$ columns of $\widehat{V}$ where $\widehat{U}$, $\widehat{\Sigma}$ and $\widehat{V}$ are the terms from (ref) of the singular value decomposition of $Y_{(n)} - X_{(n)}^\prime \beta$. Because this error term is a function of $\beta$, an iteration is required between estimating $\beta$ and finding the singular value decomposition of the error term. This is a well studied iteration; for details see Bai2009 or MoonWeidner2015.
In the following assumptions let $\widehat{r}_n$ be the estimated number of factors for the $(n)$-flattening of the regression line when applying the least square methods in (ref). Also, let $\mathcal{L}\subset \{1,\dots,d\}$ be a non-empty subset of the dimensions. In the following, the multilinear rank of $\boldsymbol{\mathcal{A}}$ is restricted such that it is low-rank along at least one of the flattenings.
Assumptions (ref), (ref) and (ref) are standard regularity assumptions already well established in the literature, e.g. see MoonWeidner2015. Assumption (ref).(i) ensures that the covariates have bounded norms, for example having bounded second moments. Assumption (ref).(ii) allows for some weak correlation across dimensions, see MoonWeidner2015, or is otherwise implied if the noise terms are independently distributed with bounded fourth moments, see latala2005some. Assumption (ref) is implied if $X_{i_1,i_2,\dots,i_d;k}\varepsilon_{i_1,i_2,\dots,i_d}$ are zero mean, bounded second moment and only admits weak correlation across dimensions for each $k = 1,\dots, K$. Assumption (ref) states that, after factor projection, the set of covariates still collectively admit full-rank variation.
Assumption (ref) is new and asserts that there exists at least one flattening of the interactive term, $\boldsymbol{\mathcal{A}}$, that is low-dimensional or simply low-rank. This requires that at least one of the unobserved terms $\varphi^{(n)}$ is low dimensional. Note that not all dimensions must satisfy Assumption (ref) for the below result. If the correct dimension is chosen then variation from the interactive term can be sufficiently projected out using the factor model approach. This makes up the statement of the following Proposition.
Proposition (ref) follows directly from MoonWeidner2015 since the flattening procedure reduces the problem to the standard linear interactive fixed-effects model. This result only applies to estimates in the dimension(s) that satisfy the low-rank assumption in Assumption (ref), i.e., the analyst has chosen the correct dimension to flatten. The constraint $\widehat{r}_{n^*}\geq r_{n^*}$ can also be changed to $\widehat{r}_{n^*}\geq c$; however, this is more conservative than required for the statement of the result. This constraint does not require knowledge of $r_{n^*}$, just that the number of estimated factors is greater than or equal the true number. The estimation procedure from Proposition (ref) can also be augmented to flatten over multiple indices, but makes Assumption (ref) harder to satisfy. For example, take the tensor $\boldsymbol{\mathcal{A}}$ flattened over the first two indices as $\mathcal{A}_{(1,2)}\in \mathbb{R}^{N_1N_2\times\prod_{n\notin\{1,2\}}N_n}$. If the parameters $\varphi^{(n)}$ for $n = 3,\dots,d$ are high-dimensional, Assumption (ref) is only satisfied when both $\varphi^{(1)}$ and $\varphi^{(2)}$ and their product space is low-dimensional. However, flattening along multiple dimensions can improve the convergence rate in Proposition (ref) to $O_p\left(\min\{N_{1}N_{2}, \prod_{n\notin\{1,2\}}N_n\}\right)^{-1/2}$.
Presented here is a simplified version of the estimator, where weights are formed from normed difference over a whole vector. This presents a curse of dimensionality that is solved with an iterative version of the estimator. Indeed, appendix (ref) presents an iterative version of the estimator that is theoretically more tractable, albeit more complicated. The iterated version is a backfitting version of the estimator presented here. Simplifying assumptions are made in this section to avoid complexities related to backfitters.
Let $\widehat{\varphi}^{(n)}_{i_n}\in\widehat{\Phi}_n$ generically denote a known or estimated proxy for fixed-effect ${\varphi}^{(n)}_{i_n}$. Let $\mathcal{W}$ be an ordered set of weight matrices, where the $n$\textsuperscript{th} item ${W}_{n}\in \mathbb{R}^{N_n \times N_n}$ has elements,
where $k$ is a kernel function, and $h_n$ is a bandwidth parameter. The weighted-within transformation in (ref) can be generalised with the following series of $n$-mode products,
likewise also on each $\mathbf{X}_k$, where $M_n = \mathbb{I}_{N_n} - W_{n}$ and $\times_n$ is the $n$-mode product defined in Section (ref). Define the weighted-within estimator as,
Proxies, $\widehat{\varphi}^{(n)}_{i_n}$, estimated from matrix method in Section (ref) are considered for inference in Section (ref), and is a more challenging setting than when proxies are observed. Further discussion of proxy estimation is relegated to Section (ref).
Assumption (ref).(i)-(iv) are standard kernel function restrictions. Assumption (ref).(v) regularises the kernel function to decay sufficiently quickly.
Assumption (ref).(i) is analogous to Assumption (ref), appropriated to the weighted-within projection. This can be verified if $\mathbf{X}$ admits a component with sufficiently independent variation across $i_1,\dots,i_d$, e.g. an i.i.d. noise component. Assumption (ref).(ii) requires weak exogeneity in the covariates after the weighted-within transformation, which can be viewed as similar to Assumption (ref). Assumption (ref).(ii) is achieved with weakly dependent distributed errors and a sample split, similar to freeman2023linear, see Section (ref) including Lemma (ref). It is, however, a strong simplifying assumption in view of the curse of dimensionality when kernel functions are evaluated over normed vectors, not scalars. Assumption (ref).(ii) could be stated as $O_p(\prod_nN_n^{-3/4}h_n^{-L/2})$ to account for the $L$-dimensional vector. However, the backfitting approach in Appendix (ref) alleviates this succinctly, so this complication is ignored here.
Proposition (ref) shows convergence for the kernel estimator is bounded by convergence of the proxy estimates. Smaller bandwidths $h_n$ lead to higher finite sample variance, and also make Assumption (ref).(i)-(ii) more difficult to satisfy, hence it is set no smaller than $O(C_n^{-1})$. Conditions in Section (ref), and discussion in Section (ref) suggest $O_p(C_{n^*}^{-1})$ can be $O_p\big(1/\sqrt{N_{n^*}}\big)$ under regularity conditions imposed in Bai2009. Hence, the parametric rate is attainable if $\mathcal{M}^\prime = \{1,\dots,d\}$ and $L$ is fixed and bounded. Implicit in the conditions for $C_{n^*}^{-1} = 1/\sqrt{N_{n^*}}$ is that the correct multilinear rank of the interactive fixed-effects is known, at least for the dimensions $\mathcal{M}^\prime$. This can likely be relaxed to the case where the upper bound on the multilinear rank is known, which would follow from the results in MoonWeidner2015 used in Proposition (ref), but left for further research.
This section establishes asymptotic normality for the main estimator, named inference corrected estimator.\footnote{The procedures here may not in general translate to the matrix low-rank approximation estimator, since convergence rates can be too slow for that estimator. } Proposition (ref) establishes an upper bound on the kernel weighted fixed-effect estimator convergence rate that can be refined to exactly the parametric rate for the fixed-effect asymptotic bias component. To ensure the bias from the fixed-effect term converges sufficiently quickly a further correction is used, proposed here.
Consider again conditional expectations, $\mathbb{E}\left(\boldsymbol{Y}|\boldsymbol{\mathcal{A}}\right) = \Gamma_{Y}$, $\mathbb{E}\left(\boldsymbol{X}|\boldsymbol{\mathcal{A}}\right) = \Gamma_{X}$, such that,
with $\mathbb{E}\left(\boldsymbol\eta|\boldsymbol{\mathcal{A}}\right) = 0$, $\mathbb{E}\left(\boldsymbol{\varepsilon}|\boldsymbol{X},\boldsymbol{\mathcal{A}}\right) = 0$. As before, the display for $\boldsymbol{X}$ is general, in that if $\boldsymbol{X}$ is unrelated to $\boldsymbol{\mathcal{A}}$, $\boldsymbol{\Gamma}_{X} = 0$. However, for identification $\boldsymbol\eta$ must admit positive sum of squares.
Repeated here, the (infeasible) inference corrected estimator, $\widehat\beta_{IC}^{(infeasible)}$, is
where ${\rm vec}_K(\boldsymbol{X} - \boldsymbol{\Gamma}_{X})\in \mathbb{R}^{\prod_n N_n \times K}$ is shorthand for the matrix of vectorised covariates, with each column the vectorised $(\boldsymbol{X}_k - \boldsymbol{\Gamma}_{X_k})$. The ${\rm vec}$ notation for $(\boldsymbol{Y} - \boldsymbol{\Gamma}_{Y})$ is the standard vectorisation. This is so far infeasible because ${\boldsymbol{\Gamma}}_{X}$, and ${\boldsymbol{\Gamma}}_{Y}$ need to be estimated.
Consider estimators for ${\boldsymbol{\Gamma}}_{X}$, and ${\boldsymbol{\Gamma}}_{Y}$ as $\widehat{\boldsymbol{\Gamma}}_{X}$, and $\widehat{\boldsymbol{\Gamma}}_{Y}$, respectively, and $\widehat{\Omega}_X = N^{-1}{\rm vec}_K(\boldsymbol{X} - \widehat{\boldsymbol{\Gamma}}_{X})^\prime {\rm vec}_K(\boldsymbol{X} - \widehat{\boldsymbol{\Gamma}}_{X})$. Then, with $\widehat{\boldsymbol{\eta}}:= \boldsymbol{X} - \widehat{\boldsymbol{\Gamma}}_{X}$,
where $\widehat{\boldsymbol{\mathcal{A}}}$ is the preliminary estimator of $\boldsymbol{\mathcal{A}}$ used to form $\widehat{\boldsymbol{\Gamma}}_{Y}$. If $\widehat{\Omega}_X^{-1}N^{-1}{\rm vec}(\boldsymbol{\mathcal{A}} - \widehat{\boldsymbol{\mathcal{A}}}) = o_p\left({\prod_{n = 1}{N_n^{-1/2}}}\right)$, then standard asymptotic normality arguments follow. Hence, convergence rates on $({\boldsymbol{\Gamma}}_{X} - \widehat{\boldsymbol{\Gamma}}_{X})$, and $(\boldsymbol{\mathcal{A}} - \widehat{\boldsymbol{\mathcal{A}}})$ are studied.
From Section (ref), with the kernel weighted estimator, $\Big(\big(\prod_{n = 1}{N_n^{-1}}\big){\rm vec}(\boldsymbol{\mathcal{A}} - \widehat{\boldsymbol{\mathcal{A}}})^\prime {\rm vec}(\boldsymbol{\mathcal{A}} - \widehat{\boldsymbol{\mathcal{A}}})\Big)^{1/2} = O_p(\prod_n h_n)$, where $h_n \gtrsim N_n^{-1/2}$ is the bandwidth used for dimension $n$. Define,
Under Assumption (ref), $\widehat{\Omega}_X^{-1} = O_p(1)\mathbbm{1}_{K\times K}$, and iid $\boldsymbol\varepsilon, \boldsymbol\eta$ it can be shown,
and under weak dependence structures in $\boldsymbol\varepsilon, \boldsymbol\eta$,
By Cauchy-Schwarz $\xi_{X\mathcal{A}} \leq \xi_\mathcal{A}\xi_{{X}}$ which is $O_p(N^{-1/2}\xi_{{X}})$ by Proposition (ref). Hence, under iid noise terms, $O_p(\xi_{X\mathcal{A}}) = o_p(N^{-1/2})$ and $\{\xi_{\mathcal{A}}, \xi_{{X}}\} = o_p(1)$, asymptotic bias is sufficiently small. Under weak dependence, $O_p(\xi_{X\mathcal{A}}) = o_p(N^{-1/2})$, and $\{\xi_{\mathcal{A}}, \xi_{{X}}\} = o_p(N^{-1/4})$, is required. Proposition (ref) establishes the upper bound on convergence for $\xi_{\mathcal{A}}$ can be $O_p(N^{-1/2})$, such that for any $\xi_X = o_p(N^{-1/4})$, $o_p(N^{-1/2})$ convergence rate for asymptotic bias is achieved.
Dependence between either $\boldsymbol\eta$, respectively $\boldsymbol{\varepsilon}$, and $\widehat{\boldsymbol{\Gamma}}_{X}$ or $\widehat{\boldsymbol{\mathcal{A}}}$ can occur because the estimator $\widehat{\boldsymbol{\Gamma}}_{X}$, and/or $\widehat{\boldsymbol{\mathcal{A}}}$, may be functions of $\boldsymbol\eta$, respectively $\boldsymbol{\varepsilon}$. To break this dependence, a simple sample splitting procedure can be implemented, e.g., as outlined in freeman2023linear, which is transferable to the independently distributed error setting.
Below regularity conditions from Bai2009,BaiNg2002 ensure estimates of the fixed-effects used as proxies in the kernel weights converge at the sufficient rate.
Below are sampling assumptions to achieve asymptotic normality.
Assumption (ref) implies there exists a consistent estimate of $\mathbb{E}\left(\boldsymbol{X}|\boldsymbol{\mathcal{A}}\right)$ at rate $N^{-1/4}$. For iid idiosyncratic terms $\boldsymbol\varepsilon, \boldsymbol\eta$ this can be relaxed to consistency at arbitrary rate. Hence, in the context of Proposition (ref), this is not an overly strong assumption on $\widehat{\boldsymbol{\Gamma}}_{X_k}$. Indeed, by Proposition (ref) and Corollary (ref) the convergence rates in Assumption (ref) are $\xi_{\mathcal{A}} = O_p(1/\sqrt{N})$. Assumption (ref) details a reprieve from requiring $\xi_{\mathcal{A}} = O_p(1/\sqrt{N})$ if indeed $\xi_X$ converges at a faster rate.
The next assumption restricts $\boldsymbol\eta$ and $\boldsymbol{\varepsilon}$ to be independent of $\widehat{\boldsymbol{\Gamma}}_{X}$ and $\widehat{\boldsymbol{\mathcal{A}}}$. As mentioned already, a sample splitting device can achieve this.\footnote{Appendix (ref) details a sample splitting device that can be used for $\widehat{\boldsymbol{\Gamma}}_{X}$ and $\widehat{\boldsymbol{\mathcal{A}}}$ in the presence of weakly dependent idiosyncratic terms. } Let $\protect\mathpalette{\protect\independenT}{\perp}$ denote stochastic independence.
Independence between $\boldsymbol{\eta}_{k,i_1\dots i_d}$ and ${\boldsymbol{\mathcal{A}}}_{i_1\dots i_d}$ in Assumption (ref).(i) is made for a cleaner derivation, but can be weakly dependent under further regularity conditions on $\boldsymbol{\eta}_{k,i_1\dots i_d}$. Assumption (ref).(ii) implies in (ref), $\widehat{\Omega}_X^{-1}{\rm vec}_K(\widehat{\boldsymbol{\eta}})^\prime{\rm vec}(\boldsymbol{\varepsilon}) = \widehat{\Omega}_X^{-1}{\rm vec}_K({\boldsymbol{\eta}})^\prime{\rm vec}(\boldsymbol{\varepsilon}) + o_p(N^{-1/2})$.
With Proposition (ref) assumptions, Assumptions (ref)-(ref), and with $C_n^{-1} \asymp h_n = O_p(1/\sqrt{N_n})$ for all $n = 1,\dots,d$ and $\mathcal{M} = \{1,\dots,d\}$, maintaining that $N:= \prod_{n = 1}{N_n}$,
In the presence of both heteroskedasticity and correlation in all dimensions, the most flexible specification to be considered here, a central limit theorem is required for $(N)^{-1/2}{\rm vec}_K({\boldsymbol{\eta}})^\prime{\rm vec}(\boldsymbol{\varepsilon})$. In this case the variance is,
Stronger independent and identically distributed assumptions are considered first, which are sequentially relaxed for stronger results later in the section.
Consistent estimation of the variance term is possible under the same set of assumptions. Estimate $\widehat{\boldsymbol{\varepsilon}} = \boldsymbol{Y} - \boldsymbol{X}\cdot\widehat\beta_{IC} - \widehat{\boldsymbol{\mathcal{A}}}$. Then, $\frac{1}{N}\sum_{i_1\dots i_d}\widehat{\varepsilon}_{i_1\dots i_d}^2 = \sigma_\varepsilon^2 + o_p(1)$ follows from the consistency of $\widehat{\boldsymbol{\Gamma}}_{Y}$, $\widehat{\boldsymbol{\Gamma}}_{X}$, and $\widehat\beta_{IC}$. Likewise $\mathbb{E}\left[ \boldsymbol{\eta}_{i_1\dots i_d}^{\,} \boldsymbol{\eta}_{i_1\dots i_d}^{\prime} \right]$ can be estimated consistently from the sample analog $\frac{1}{N}\sum_{i_1\dots i_d} \left(\boldsymbol{X}_{i_1\dots i_d} - \widehat{\boldsymbol{\Gamma}}_{X, _{i_1\dots i_d}}\right) \left(\boldsymbol{X}_{i_1\dots i_d} - \widehat{\boldsymbol{\Gamma}}_{X, _{i_1\dots i_d}}\right)^\prime $.
Assumptions on the error can be used to satisfy Lyapunov conditions. Instead, a central limit theorem is assumed, as is common in the interactive fixed-effects literature.\footnote{See for example Assumption E in Section 5 in Bai2009. }
Under correlation, independence in Assumption (ref) needs a stronger sample split than from freeman2023linear -- a novel sample split device is devised in Appendix (ref).
Assumption (ref) can be achieved if, e.g. $\{\xi_X,\xi_\mathcal{A}\} = o_p(N^{-1/4})$, without further restrictions on $\xi_{\sigma,N}$. If $\varepsilon_{i_1\dots i_d}$ are i.i.d., or i.n.i.d., $\xi_{\sigma,N}=0$, satisfying the assumption.
From Bai2009, $\widetilde{\Sigma}$ can be estimated with the newey1987hac kernel approach for time dimension, in conjunction with a partial sample estimator for cross-sections.
This section serves to discuss the results in Section (ref), motivate further some of the chosen methods, and provide some methods to estimate proxies for kernel weights.
Proposition (ref) takes for granted the dimension to flatten across admits a low rank interactive fixed-effect term for the least square method in Bai2009. Established diagnostics in BaiNg2002, ahn2013eigenvalue and hallin2007determining can be used to determine the number of factors in the pure factor model, which could be repeated along the different flattenings in the multidimensional setting. However, establishing the interactive fixed-effect rank in the presence of covariates is an open area of research.
Consider the factor model applied to a flattening that may not be low-rank. First, assume $\varphi^{(1)}$ varies in a high-dimensional parameter space, e.g. with $N_1 < N_2N_3$, $\varphi^{(1)} \in \mathbb{R}^{N_1\times N_1}$ and $\Gamma\in\mathbb{R}^{N_2N_3\times N_1}$ with each column mutually orthogonal for both these matrices. Then $\varphi^{(1)}\Gamma^\prime$ is full-rank and the factor model will not control for this. On the contrary, consider $\varphi^{(1)} \in \mathbb{R}^{N_1\times N_1}$ where all columns are linearly dependent. Then the matrix $\varphi^{(1)} \Gamma^\prime$ is rank-1 regardless of $L$ and of how $\varphi^{(2)}$ and $\varphi^{(3)}$ vary, thus can be projected with a factor model estimated with 1 factor. This situation is exemplified in simulations in Section (ref). The matrix method requires knowledge of this low rank dimension - the kernel weighted-within does not.
Knowing the multilinear rank, MoonWeidner2015 show a factor model with at least $r_n$ factors results in consistent $\beta$ estimates.\footnote{This actually requires knowing an upper bound on multilinear rank, not the actual multilinear rank. } However, with generic heteroskedasticity and weak correlations, the factor model can only be bounded by a slow rate of convergence, which is exemplified in simulations below. In the three dimensional setting with proportional dimensions, this bound is as slow as $N^{-1/6}$, where $N$ is total sample size. This is too slow to use for standard inference procedures.\footnote{Debias do exist, but current proposals would achieve at best $N^{-1/3}$ convergence. There may be split sample debias innovations that can achieve sufficiently small asymptotic bias, but this is not explored. }
Discussed here are some functionals of multidimensional arrays that are useful for estimating proxies to form kernel weights. For this purpose, it is important to find proxies that isolate variation in each dimension since weights are formed separately for each index.
Consider for any of the dimensions $n$ the corresponding matrix of left singular vectors from above, $U^{(n)}$, estimated with noise $\varepsilon_{ijt}$. That is, each $U^{(n)}$ are calculated from the matrix $\boldsymbol{\mathcal{V}}_{(n)} = \boldsymbol{\mathcal{A}}_{(n)} + \boldsymbol{\varepsilon}_{(n)}$. Under regularity conditions on the noise term $\boldsymbol{\varepsilon}$, the left singular vectors from this decomposition consistently estimate $U^{(n)}$, up to rotations. In the three dimensional case, define the $r_n$-vector $\widehat{U}^{(n)}_{i}$ as the $i$-{th} row of the left singular matrix of $\boldsymbol{\mathcal{A}} + \boldsymbol{\varepsilon}$ flattened in the $n$\textsuperscript{th} dimension. BaiNg2002 show the following “up-to-rotation” consistency result:
Hence, if the true error, $\boldsymbol{\mathcal{V}} = \boldsymbol{\mathcal{A}} + \boldsymbol{\varepsilon}$, is observed, this establishes consistency of $\widehat{U}^{(n)}_{i_n} $ as proxies. Matrices, $H^k_n$, in Lemma (ref) can be ignored as these do not impact convergence rates, see end of Section (ref). However, the true error, $\boldsymbol{\mathcal{V}} = \boldsymbol{\mathcal{A}} + \boldsymbol{\varepsilon}$, is not observed. Instead, the observed error $\widehat{\mathcal{V}}_{ijt} = Y_{ijt} - X_{ijt}\widehat{\beta} = X_{ijt}(\beta^0 - \widehat{\beta}) + \mathcal{A}_{ijt} + \varepsilon_{ijt}$, depends on the estimate $\widehat{\beta}$, hence should be accounted for. Proposition A.1 in Bai2009 provides,
if $\widehat{U}^{(n)}_{i_n}$ come from least squares in Bai2009. Ignoring $H^{k}_n$, $\frac{1}{N_n}\sum_{i_n = 1}^{N_n}\left\lVert\widehat{U}^{(n)}_{i_n} - \varphi^{(n)}_{i_n}\right\rVert^2 = O_p\left( {N_n}^{-1}\right)$ if each dimension sample size grows at a roughly comparable rate and the convergence rate for $\|\beta^0 - \widehat{\beta}\|$ in Proposition (ref) is used.\footnote{ To be precise, the comparable rate restriction for any dimension $n = 1,\dots,d$ is $N_n \lesssim \prod_{n^\prime \neq n} N_{n^\prime}$. The comparable rate commonly used in the two-dimensional panel literature is $N\approx T$. Hence, in the multidimensional problem, the relative rate requirement is milder than the two-dimensional case. } Then, $C_n^{-2}$ from Proposition (ref) is $1/N_n$ and the parametric rate for the kernel weighted estimator can be established.
As a special case, consider the weighted-within transformation with a linear kernel, $k: \mathcal{U}\times \mathcal{U} \to \mathbb{R}$ using estimated singular vectors $\widehat{U}_{i_n}^{(n)}$: $k(u_i, u_j) = u_i'u_j$. Then in the transformation,
each $W_n$ is replaced by $\widehat{U}^{(n)}\widehat{U}^{(n)\,'}$. By construction, these singular vectors are orthonormal, hence, this is equivalent to $W_n = \widehat{U}^{(n)}\left(\widehat{U}^{(n)\,'}\widehat{U}^{(n)} \right)^{-1} \widehat{U}^{(n)\,'}$. Then,
is simply a multilinear projection of $\{\widehat{U}^{(n)}\}_{n=1}^d$. The error can be very simply bounded in this case by arguments presented in Appendix (ref). This paper explores the more general case under generic kernel functions since it likely has applications extending past the parametric model for interactive effects. However, the linear kernel case is included for completeness.
Two simulation exercises are considered. The first analyses performance when sample size is allowed to increase. The second analyses fixed sample properties of the estimators when Assumption (ref), which governs fixed-effect multilinear rank, is violated in some dimensions.
Data is generated as,
with, $\mathcal{A}_{ijt} = \sum_{\ell = 1}^{2} \lambda_{i\ell} \gamma_{j\ell} f_{t\ell} $. Also,
Error $\varepsilon_{ijt}$ admits heteroskedasticity, and correlation in each dimension of the data. The order of $i$ and $j$, are randomised to simulate an unknown cross correlation structure.
The results depicted in Figure (ref) show two specifications for $\rho$. The left panel depicts $\rho = 1$, and the right panel depicts $\rho = 1/3$. Estimators considered use 2 estimated factors such that the multilinear rank is correctly predicted. The left panel shows that whilst the matrix methods, labelled Factor, look consistent, they converge at a rate much slower than the empirical bounds. The right panel shows more egregious convergence when the bias problem is made worse with $\rho = 1/3$. Indeed, for the sample sizes considered, the empirical bounds do not ever cover the true parameter value in either panel.
This contrasts with the results for the kernel weighted fixed-effect projection. Finite sample bias is negligible compared to variance in both the left and right panel. This estimator has roughly the same variance as the matrix methods, and for the sample sizes considered, the empirical bounds are always valid.
Table (ref) displays coverage and compares the weighted estimator, also with the inference correction from Section (ref), and the factor model applied to all dimensions. The W Inf coverage, which is coverage for the inference corrected estimator, is always better than the uncorrected estimator, and achieves nominal coverage for larger sample sizes when the multilinear rank is exactly or over estimated. There is evidence of slight over-coverage when the rank is overestimated, but this is expected when undersmoothing. Coverage for the uncorrected estimator improves with overestimation of the rank. The factor model undercovers with progressively worse coverage as the sample size grows.
Table (ref) shows simulation results for the following DGP,
with $\mathcal{A}_{ijt} = \sum_{\ell = 1}^{N_1} \lambda_{i\ell} \gamma_{j\ell} f_{t\ell} $ and all other parameters the same as (ref). $\mathcal{A}_{ijt}$ is normalised to have unit variance. $\boldsymbol{\mathcal{A}}$ is specified such that it is rank 1 when flattened in the first dimension and rank $N_1$ when flattened in either dimension two or three. That is, the multilinear rank is $\mathbf{r} = (1,N_1,N_1)$. To achieve this, the matrix $\lambda$ is designed to be rank-1 and the matrices $\gamma$ and $f$ are designed to be rank-$N_1$.
In Table (ref), the estimators OLS and Fixed-effects are simply the pooled OLS estimator and the pooled OLS estimator after additive fixed-effects are projected out, respectively. As expected both of these have bias. The factor model is used after first flattening along each dimension as Factor(dim = $n$), where $n$ is the dimension used for flattening. In each case, 2 factors are projected. The results show the bias is close to zero when the correct dimension is flattened over (the first dimension in this case) and very poor bias when the incorrect dimension is used (the second and third dimensions). Lastly, the weighted differencing estimator is estimated with Gaussian kernel function with various bandwidths; which are standardised to be equivalent to standard deviations of the proxy measures. Weighted estimators show a clear bias-variance trade-off, but all with bias of the same order as the correct factor model.
This analysis is repeated for the four dimensional case in Table (ref), where the first and second dimensions admit low-dimensional unobserved interactive fixed-effects parameters. The simulations suggest similar results as the three dimensional case, where the factor models perform well when flattened in the low-dimensional dimensions (first and second) and poorly in the high-dimensional dimensions (third and fourth).
This paper develops methods to generalise the interactive fixed-effect to multidimensional datasets with more than two dimensions. Theoretical results show that standard matrix methods can be applied to this setting but require additional knowledge of the data generating process and generally have slow convergence rates. Nonetheless, these provide useful preliminary estimates. The multiplicative interactive fixed-effect error from the kernel weighted method show an improvement on the convergence rate of slope coefficient estimates to the parametric rate, and suggest a more robust approach to projecting fixed-effects. Simulations show finite sample properties when the sample size is allowed to grow and when it is fixed. These simulations exemplify how sensitive the two-dimensional methods are to model specification issues as simple as how to organise the dataset. They also show the robustness of the kernel weighted fixed-effects estimators without having to make these same specifications. A method for inference with the weighted fixed-effects estimator is also introduced.
The model is applied to a demand model for beer consumption. The application demonstrates that simply applying the two-dimensional factor model approach is sensitive to how dimensions are arranged to suit these estimators. The weighted-within transformation estimates elasticities close to an instrumental variable point estimate, but with substantially better precision.