The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
84,567 characters
Efficient Estimation by Fully Modified GLS with an Application to the Environmental Kuznets Curve
\maketitle
\begin{abstract}
\noindent
This paper develops the asymptotic theory of a Fully Modified Generalized Least Squares estimator for multivariate cointegrating polynomial regressions. Such regressions allow for deterministic trends, stochastic trends and integer powers of stochastic trends to enter the cointegrating relations. Our fully modified estimator incorporates: (1) the direct estimation of the inverse autocovariance matrix of the multidimensional errors, and (2) second order bias corrections. The resulting estimator has the intuitive interpretation of applying a weighted least squares objective function to filtered data series. Moreover, the required second order bias corrections are convenient byproducts of our approach and lead to standard asymptotic inference. We also study several multivariate KPSS-type of tests for the null of cointegration. A comprehensive simulation study shows good performance of the FM-GLS estimator and the related tests. As a practical illustration, we reinvestigate the Environmental Kuznets Curve (EKC) hypothesis for six early industrialized countries as in \cite{wagnergrabarczykhong2019}.
\bigskip
\noindent
\textbf{JEL Classification}: C12, C13, C32, Q20
\bigskip
\noindent
\textbf{Keywords}: Cointegrating Polynomial Regression, Cointegration Testing, Environmental Kuznets Curve, Fully Modified Estimation, Generalized Least Squares
\end{abstract}
\newpage
\section{Introduction}
In recent years, there has been an increasing interest in the theoretical properties and theoretical justifications of nonlinear cointegrating relations. For theoretical properties we refer to the textbook treatise by \cite{wang2015} , the recent review article by \cite{tjostheim2020}, and the extensive references found in either of them. Theoretical justifications are in some cases refinements of existing economic theory, e.g. nonlinear cointegration among bond yields with different times to maturity due to yield-dependent risk premia as discussed in \cite{breitung2001}, or nonlinear purchasing power parity due to transaction/transportation costs and trade barriers (e.g. \cite{hongphillips2010}). In other cases, economic theory postulates a nonlinear cointegrating relation from the outset. A popular example of the latter is the Environmental Kuznets curve described in \cite{grossmankrueger1995}.\footnote{There is no direct reference to Kuznets in the original paper by \cite{grossmankrueger1995}. But their nonlinear relations between environmental indicators and per capita GDP do remind strongly of the inverted U-shaped between income inequality and economic growth proposed by Kuznets (1901-1985). The term `environmental Kuznets curve' was used later.}
There are three branches of literature on the estimation of such nonlinear cointegrating relations. First, the papers by \cite{parkphillips1999} and \cite{parkphillips2001} are concerned with nonlinear cointegration analysis of a parametric form. Second, there is a literature on nonparametric kernel estimation of nonlinear cointegrating equations, see for example \cite{wangphillips2009} or \cite{lihillipsgao2017}. The third approach is reminiscent of a nonparametric sieve estimation with power polynomial basis. That is, one estimates a cointegrating relation containing integer powers of integrated regressors. \cite{wagnerhong2016} named this a cointegrating polynomial regression (CPR). The multivariate seemingly unrelated regressions extension is available in \cite{wagnergrabarczykhong2019}. Our model specification builds on this Seemingly Unrelated Cointegrating Polynomial Regression (SUCPR) setup.
We make two theoretical contributions to the literature on cointegrating polynomial regressions. First, we propose the Fully Modified Generalized Least Squares (FM-GLS) estimator. This estimator requires two main steps: (1) It employs the inverse covariance matrix of the $2nT$-dimensional innovation vector, that is, the covariance matrix of the vector which stacks the $n$ disturbances in the cointegrating equations and the $n$ disturbances driving the $I(1)$ regressors over the time span $T$. The estimation of this inverse covariance matrix is based on the Modified Cholesky Decomposition (MCD) originating from Pourahmadi (1999). The approach is computationally simple because the required quantities are obtained from the coefficients and prediction error variances of best linear least squares predictors. In our setting this translates into estimating multiple VAR models up to some maximum lag order $q$. Sufficient conditions for consistency are provided. (2) We exploit the previous results to correct the second-order biases, resulting in improved efficiency and standard chi-square inference. Also note that the approach differs from the linear cointegration results in \cite{markogakisul2005} and \cite{moonperron2005} since our bias corrections do not rely on leads and lags augmentation. Second, a multi-equation cointegration specification asks for a multivariate cointegration test. Building upon the work by \cite{choisaikkonen2010}, we propose three such tests. The first test uses pre-filtered residuals to account for serial correlation, whereas the other two are direct multivariate generalizations of the KPSS-type of test in \cite{wagnerhong2016}. The estimator and cointegration tests are subsequently studied by Monte Carlo simulation. In our simulations, the FM-GLS estimator has a higher estimation accuracy and its implied Wald test has better size control and higher size-adjusted power. We find by simulation that prefiltering improves the size control of the cointegration tests but has an adverse effect on power. In the empirical application there is a surprisingly large spread in the widths of the confidence intervals. It turns out that FM-SUR, and to a lesser degree FM-SOLS, underestimates the parameter uncertainty compared to FM-GLS.
The plan of this paper is as follows. Section \ref{sec:model} introduces the model and the modified Cholesky block decomposition. This decomposition is the main ingredient for the fully modified GLS estimator. The related asymptotic theory and stationarity tests are discussed in Section \ref{sec:asymptotictheory} whereas a finite sample simulation study is presented in Section \ref{sec:simulations}. The empirical application can be found in Section \ref{sec:empappl} where we look at the environmental Kuznets curve. Section \ref{sec:conclusion} concludes. All proofs are collected in the Appendices.\footnote{The Appendices contains the proofs of all the results that are related to the generalized least squares estimator. Supplementary material is available on the websites of the authors.}
Some words on notation. $C$ denotes a generic positive constant. The integer part of the number $a\in \SR^{+}$ is denoted by $[a]$. For a vector $\vx\in\SR^n$, its dimension is abbreviated by $\dim(\vx)$ and its $p$-norm by $\|\vx\|_p=(\sum_{i=1}^{n}|x_i|^p)^{1/p}$. When applied to a matrix, $\|\mA\|_p$ signifies the induced norm defined by $\|\mA\|_p=\sup_{\vx\neq \vzeros} \|\mA\vx\|_p/\|\vx\|_p$. The subscripts are omitted whenever $p=2$, e.g. $\|\vx\|=\left(\sum_{i=1}^{n}|x_i|^2\right)^{1/2}$ and $\|\mA\|=\left(\lambda_{max}\left(\mA'\mA\right)\right)^{1/2}$ where $\lambda_{max}(\cdot)$ is the largest eigenvalue. Similarly, $\lambda_{min}(\cdot)$ denotes the smallest eigenvalue. The Frobenius norm is denoted as $\|\cdot\|_{\calF}$. The $(n\times n)$ identity matrix is written as $\mI_n$. The $i^{th}$ row or $i^{th}$ column of an arbitrary matrix $\mA$ are selected using $\col_i(\mA)$ and $\operatorname{row}_i(\mA)$, respectively. The Kronecker product is denoted ``$\otimes$''. We use the symbol ``$\wto$'' to signify weak convergence and the symbol ``$\stackrel{d}{=}$'' for equality in distribution. The stochastic order and strict stochastic order relations are indicated by $O_p(\cdot)$ and $o_p(\cdot)$.
\section{The Model}\label{sec:model}
As in \cite{wagnergrabarczykhong2019}, we study a system of seemingly unrelated cointegrating polynomial regressions (SUCPR), that is
\begin{equation}
\begin{aligned}
\vy_t =\mZ_t' \vbeta+\vu_t,\qquad\qquad \text{for }t=1,2,\ldots, T,
\label{eq:basemodel}
\end{aligned}
\end{equation}
where the dependent variable $\vy_t:=[y_{1t},y_{2t},\ldots,y_{nt}]'$ and innovations $\vu_t:=[u_{1t},u_{2t},\ldots,u_{nt}]'$ are $(n\times 1)$ random vectors. For the cross-sectional unit $i$, we use as explanatory variables: (1) deterministic components such as an intercept and polynomial time trends up to order $d_i$, and (2) integer powers of the $I(1)$ regressors $x_{it}$ up to degree $s_i$. Defining $\vd_{it}=[1,t,\ldots,t^{d_i}]'$, $\vs_{it}=[x_{it},\ldots,x_{it}^{s_{i}}]$, and $\vz_{it}=[\vd_{it}',\vs_{it}']'$, we subsequently collect all explanatory variables in the block diagonal matrix $\mZ_t=\diag[\vz_{1t},\ldots,\vz_{nt}]$. We are interested in the $d$-dimensional parameter vector $\vbeta$ where $d=\sum_{i=1}^n (d_i+s_i+1)$. Overall, each cross-sectional unit in \eqref{eq:basemodel} specifies a single cointegrating relation containing polynomials in deterministic and stochastic trends. For each $i$, the highest orders of these polynomials, i.e. $d_i$ and $s_i$, are assumed to be fixed and known. We do not allow for cointegration in the cross-sectional dimension.
The innovation series $\{\vu_t\}$ is allowed to exhibit dependencies over time and across series. We assume that these dependencies can be modeled by a stationary VAR($\infty$) process, that is
\begin{equation}
\bm{\mathcal{A}}(L) \vu_t= \left(\mI_n- \sum_{j=1}^\infty \mA_j L^j \right)\vu_t=\veta_t,
\end{equation}
(see Assumption \ref{assumpt1:linearproc} for further details). Efficient estimation of the parameter vector $\vbeta$ now requires the use of generalized least squares (GLS). Our \cite{zellner1962}-type GLS estimator relies on the inverse of the $(nT\times nT)$ matrix $\mSigma_{\vu}=\E(\vu\vu')$ where $\vu=[\vu_1',\vu_2',\ldots,\vu_T']'$. In this paper, we directly estimate $\mSigma_{\vu}^{-1}$ using a multivariate extension of the modified Cholesky decomposition by \cite{pourahmadi1999}. This extension was named the \emph{Modified Cholesky Block Decomposition} (MCBD\label{abr:MCBD}) by \cite{kimzimmerman2012} and \cite{kohligarciapourahmadi2016}. The latter papers used the MCBD to parametrize the covariance matrix of multivariate longitudinal data. As in \cite{beutnerlinsmeekes2019}, we use the MCBD for the time series application mentioned above, i.e. the computation of $\mSigma_{\vu}^{-1}$. The decomposition is closely related to linear minimum MSE predictors.
We define
\begin{equation}
\begin{aligned}
\mA(\ell) &=
\begin{bmatrix}
\mA_1(\ell) & \cdots & \mA_\ell(\ell)
\end{bmatrix}
= \operatornamewithlimits{arg\;min}_{(\mTheta_1,\ldots,\mTheta_\ell)\in \SR^{n\times n\ell}} \E\left\|\vu_t-\mTheta_1 \vu_{t-1}-\cdots-\mTheta_\ell \vu_{t-\ell} \right\|^2, \\
\mS(\ell) &= \E\left[\vu_t-\mA_{1}(\ell)\vu_{t-1}-\cdots-\mA_{\ell}(\ell)\vu_{t-\ell} \right]\left[\vu_t-\mA_{1}(\ell)\vu_{t-1}-\cdots-\mA_{\ell}(\ell)\vu_{t-\ell} \right]',
\end{aligned}
\label{eq:populationMCD}
\end{equation}
and $\mS(0)=\E(\vu_t^{}\vu_t')$. The inverse of the covariance matrix $\mSigma_{\vu}$ is then given by
\begin{equation}
\mSigma_{\vu}^{-1}=\bm{\mathcal{M}}_{\vu}' \bm{\mathcal{S}}_{\vu}^{-1}\bm{\mathcal{M}}_{\vu}^{},
\label{eq:modCholdecomp}
\end{equation}
where $ \bm{\mathcal{S}}_{\vu} = \diag \Big(\mS(0),\mS(1),\ldots,\mS(T-1) \Big)$,
\begin{equation}
\begin{aligned}
\bm{\mathcal{M}}_{\vu} &=\left[\bm{m}_{\vu}^{ij} \right]_{1\leq i,j\leq T},\text{ with }\bm{m}_{\vu}^{ij}=
\begin{cases}
\mZeros_{n\times n},&\text{if}\quad i<j,\\
\mI_n,&\text{if}\quad i=j,\\
-\mA_{i-j}(i-1),&\text{if}\quad 2\leq i\leq T,\;1\leq j\leq i-1,
\end{cases}
\end{aligned}
\end{equation}
and the $\mA_{j}(i)$ follow from the partitioning $\mA(\ell)=\big[\mA_{1}(\ell),\ldots,\mA_{\ell}(\ell)\big]$.
Weak stationarity of $\{\vu_t\}$ implies that the block elements of $\bm{\mathcal{M}}_{\vu}$ being far below the main diagonal are small. This suggests a banding approach in which small elements are replaced by zeros. More specifically, we construct a \emph{Banded Inverse Autocovariance Matrix} (BIAM) \label{ref:BIAM} as
\begin{equation}
\mSigma_{\vu}^{-1}(q)=\bm{\mathcal{M}}_{\vu}'(q)\bm{\mathcal{S}}_{\vu}^{-1}(q)\bm{\mathcal{M}}_{\vu}^{}(q),
\label{eq:bimam}
\end{equation}
where $1\leq q\ll T$ is called the banding parameter, $\bm{\mathcal{S}}_{\vu}(q)=\diag\Big(\mS(0),\mS(1),\ldots,\mS(q),\ldots,\mS(q)\Big)$ and $\bm{\mathcal{M}}_{\vu}(q)=\left[\bm{m}_{\vu}^{ij}(q) \right]_{1\leq i,j\leq T}$ with
\begin{equation}
\vm_{\vu}^{ij}(q)=\begin{cases}
\mZeros_{n\times n},&\text{if}\quad i<j\;\text{or}\; \{q+1<i\leq T,\;1\leq j\leq i-q-1\}\\
\mI_n,&\text{if}\quad i=j\\
-\mA_{i-j}(i-1),&\text{if}\quad 2\leq i\leq q,\;1\leq j\leq i-1\\
-\mA_{i-j}(q),&\text{if}\quad q+1\leq i\leq T,\;i-q\leq j\leq i-1.
\end{cases}
\label{eq:Mmatrixconstrucion}
\end{equation}
\begin{example} \label{example:VAR1}
Consider a stationary $n$-dimensional VAR($3$) process specified as $\vu_t= \sum_{j=1}^3 \mA_j \vu_{t-j}+\veta_t$ with $\veta_t\stackrel{i.i.d.}{\sim}(\vzeros,\mSigma_{\eta\eta})$. For $T=4$, the MCBD $\mSigma_{\vu}^{-1}=\bm{\mathcal{M}}_{\vu}' \bm{\mathcal{S}}_{\vu}^{-1}\bm{\mathcal{M}}_{\vu}^{}$ is based on
\begin{equation}
\bm{\mathcal{M}}_{\vu}^{}
=
\left[
\begin{smallmatrix}
\mI_n & \mZeros & \mZeros & \mZeros \\
-\mA_1(1) & \mI_n & \mZeros & \mZeros \\
-\mA_2(2) & -\mA_1(2) & \mI_n & \mZeros \\
-\mA_3 &-\mA_2 & -\mA_1 & \mI_n
\end{smallmatrix}
\right]
,
\qquad\qquad
\bm{\mathcal{S}}_{\vu}=
\left[
\begin{smallmatrix}
\mS(0) \\
& \mS(1) \\
& & \mS(2) \\
& & & \mSigma_{\eta\eta}
\end{smallmatrix}
\right].
\end{equation}
Alternatively, with banding parameter $q=2$, the related banded inverse autocovariance matrix is $\mSigma_{\vu}^{-1}(2)=\bm{\mathcal{M}}_{\vu}'(2)\bm{\mathcal{S}}_{\vu}^{-1}(2)\bm{\mathcal{M}}_{\vu}^{}(2)$ with
\begin{equation}
\bm{\mathcal{M}}_{\vu}^{}(2)
=
\left[
\begin{smallmatrix}
\mI_n & \mZeros & \mZeros & \mZeros \\
-\mA_1(1) & \mI_n & \mZeros & \mZeros \\
-\mA_2(2) & -\mA_1(2) & \mI_n & \mZeros \\
\mZeros &-\mA_2(2) & -\mA_1(2) & \mI_n
\end{smallmatrix}
\right]
,
\qquad\qquad
\bm{\mathcal{S}}_{\vu}(2)=
\left[
\begin{smallmatrix}
\mS(0) \\
& \mS(1) \\
& & \mS(2) \\
& & & \mS(2)
\end{smallmatrix}
\right].
\end{equation}
\end{example}
The model of \eqref{eq:basemodel} can be stacked over time to yield the representation $\vy=\mZ \vbeta+\vu$ with $\vy=[\vy_1',\vy_2',\ldots,\vy_T']'$, $\mZ=[\mZ_1,\mZ_2,\ldots,\mZ_T]'$ and $\vu$ as before. For the moment, \emph{we will assume $\mSigma_{\vu}^{-1}(q)$ to be known} and focus on the following estimator:
\begin{equation}
\widehat{\vbeta}_{GLS}:=\left(\mZ' \mSigma_{\vu}^{-1}(q) \mZ \right)^{-1} \mZ' \mSigma_{\vu}^{-1}(q) \vy.
\end{equation}
A discussion on the properties of this infeasible estimator is informative because: (1) the incurred estimation error of an appropriately constructed estimator$\widehat{\mSigma_{\vu}^{-1}}(q)$ will be asymptotically negligible, and (2) we can suppress the effect of banding by letting $q$ increase with sample size.
Two remarks related to $\widehat{\vbeta}_{GLS}$ are instructive. First, the GLS estimator differs from the usual least squares estimator $\widehat{\vbeta}_{OLS}:=\left(\mZ' \mZ \right)^{-1} \mZ' \vy$ by a weighing with the inverse covariance matrix $\mSigma_{\vu}^{-1}(q)$. It is well documented in standard econometric textbooks (e.g. chapter 7 of \cite{davidsonmackinnon2004}) that this weighing may lead to substantial efficiency gains. Second, it is illustrative to substitute the Modified Cholesky Decomposition of $\mSigma_{\vu}^{-1}(q)$ into the definition of this infeasible GLS estimator. The result is $\widehat{\vbeta}_{GLS}=(\mZ_{filt}' \bm{\mathcal{S}}_{\vu}^{-1}(q) \mZ_{filt}^{} )^{-1}\mZ_{filt}' \bm{\mathcal{S}}_{\vu}^{-1}(q) \vy_{filt}$ where $\mZ_{filt}^{}=\bm{\mathcal{M}}_{\vu}(q)\mZ$, and $\vy_{filt}=\bm{\mathcal{M}}_{\vu}(q) \vy$. The premultiplications by $\bm{\mathcal{M}}_{\vu}(q)$ have the effect of filtering and take care of serial correlation. $\bm{\mathcal{S}}_{\vu}^{-1}(q)$ applies scaling and rotation to account for the correlations between the series. The following univariate autoregressive setting exemplifies this intuition.
\begin{example} \label{example:praiswinston}
A regression model $y_t=\beta t +u_t$ has AR($1$) innovations $u_t=\rho u_{t-1}+\eta_t$ where $\eta_t\stackrel{i.i.d.}{\sim}(0,\sigma^2)$ and $|\rho|<1$. Taking $n=1$, the expressions of Example \ref{example:VAR1} are easily adapted to yield:
\begin{equation}
\bm{\mathcal{S}}_{\vu}^{}=
\diag\left(\frac{\sigma^2}{1-\rho^2},\sigma^2,\ldots,\sigma^2 \right),
\qquad
\bm{\mathcal{M}}_{\vu}^{}\vy=
\begin{bmatrix}
1 \\
-\rho &1 \\
\vdots & \ddots &\ddots \\
0 & \cdots & -\rho & 1
\end{bmatrix}
\begin{bmatrix}
y_1 \\
y_2 \\
\vdots \\
y_T
\end{bmatrix}
=
\begin{bmatrix}
y_1 \\
y_2-\rho y_1 \\
\vdots \\
y_T-\rho y_{T-1}
\end{bmatrix},
\end{equation}
and a similar transformation for the linear trend. The implied GLS estimator coincides with the estimator from \cite{praiswinston1954}.
\end{example}
\section{Asymptotic Theory} \label{sec:asymptotictheory}
In this section, we present the asymptotic results. More specifically, we derive: (1) the limiting distribution of the GLS estimator, (2) the fully modified GLS (FM-GLS) estimator that corrects for second order bias terms, (3) a Wald test statistic, and (4) several multivariate KPSS-type of tests for the null of cointegration. We will also compare this FM-GLS estimator with the two fully modified estimators defined in Proposition 1 of \cite{wagnergrabarczykhong2019}. The following assumption will facilitate the development of the asymptotic theory.
\begin{assumption}[Innovation Processes]\label{assumpt1:linearproc}
The innovations processes in the model satisfy the following assumptions:
\begin{enumerate}[(a)]
\item The process $\vzeta_t^{}=[\veta_t',\vepsi_t']'$ is an independent and identically distributed (i.i.d.) sequence with $\E(\vzeta_t^{} \vzeta_t')=\left[\begin{smallmatrix} \mSigma_{\eta\eta} & \mSigma_{\eta \epsilon} \\ \mSigma_{\epsilon \eta} & \mSigma_{\epsilon \epsilon} \end{smallmatrix} \right] \succ 0$ and $\E(\|\vzeta_t\|^{2r})\leq C_r<\infty$ for some constant $C_r>0$ and some $r>2$.
\item $\det\big(\bm{\mathcal{A}}(z)\big)\neq 0$ for all $|z|\leq 1$ and $\sum_{j=0}^\infty j \| \mA_j \|_\calF<\infty$.
\item $\diff \vx_t=\vv_t$ admits the VAR($\infty$) process $\bm{\mathcal{D}}(L)\vv_t=\vepsi_t$, where $\bm{\mathcal{D}}(L)=\mI_n-\sum_{j=1}^\infty \mD_j L^j$. Moreover, $\det\big(\bm{\mathcal{D}}(z)\big)\neq 0$ for all $|z|\leq 1$ and $\sum_{j=0}^\infty j \, \|\mD_j \|_{\calF}<\infty$.
\end{enumerate}
\end{assumption}
The stationary VAR($\infty$) specifications for $\{\vu_t\}$ and $\{\vv_t\}$ are natural given the linear minimum MSE predictor formulae that underly the definitions of the MCBD and BIAM. Moreover, the conditions in Assumption \ref{assumpt1:linearproc} ensure that the lag polynomials $\bm{\mathcal{A}}(L)$ and $\bm{\mathcal{D}}(L)$ are invertible (see for example Theorem 7.4.2 of \cite{hannandeistler2012}), thereby showing that our Assumption \ref{assumpt1:linearproc} is similar to the linear processes assumptions that are regularly adopted in the literature on nonlinear cointegration, cf. \cite{choisaikkonen2010}, \cite{wagnerhong2016}, and \cite{wagnergrabarczykhong2019}. The assumption $\det\big(\bm{\mathcal{D}}(1)\big)\neq 0$ rules out cointegration among the components of $\{\vx_t\}$.
Under Assumption \ref{assumpt1:linearproc}(a), an invariance principle holds for $\vzeta_t$, i.e. $\frac{1}{T^{1/2}} \sum_{t=1}^{[rT]} \vzeta_t\wto \bm{B}_{\vzeta}(r)\equiv \left[\begin{smallmatrix} \bm{B}_\eta(r)\\ \bm{B}_\epsilon(r) \end{smallmatrix} \right]$ where $\bm{B}_{\vzeta}$ denotes an $2n$-dimensional Brownian motion with covariance matrix $\left[\begin{smallmatrix} \mSigma_{\eta\eta} & \mSigma_{\eta \epsilon} \\ \mSigma_{\epsilon \eta} & \mSigma_{\epsilon \epsilon} \end{smallmatrix} \right] $. Moreover, Assumptions \ref{assumpt1:linearproc}(b)-(c) justify the use of the Beveridge-Nelson decomposition (\cite{phillipssolo1992}). A functional central limit theorem for linear processes is thus also applicable to $\vxi_t^{}=[\vu_t',\vv_t']'$, that is
\begin{equation}
\frac{1}{T^{1/2}} \sum_{t=1}^{[rT]} \vxi_t \wto \bm{B}_{\vxi}(r)
\equiv
\begin{bmatrix}
\bm{B}_u(r) \\
\bm{B}_v(r)
\end{bmatrix}
\equiv
\begin{bmatrix}
\bm{\mathcal{A}}(1) & \mZeros \\
\mZeros & \bm{\mathcal{D}}(1)
\end{bmatrix}^{-1}
\begin{bmatrix}
\bm{B}_\eta(r)\\
\bm{B}_\epsilon(r)
\end{bmatrix},
\label{eq:brownianequivalences}
\end{equation}
where the Brownian motion $\bm{B}_{\vxi}(r)$ of dimension $2n$ has covariance matrix
\begin{equation}
\mOmega=
\begin{bmatrix}
\mOmega_{uu} & \mOmega_{uv} \\
\mOmega_{vu} & \mOmega_{vv}
\end{bmatrix}
=
\begin{bmatrix}
\bm{\mathcal{A}}(1) & \mZeros \\
\mZeros & \bm{\mathcal{D}}(1)
\end{bmatrix}^{-1}
\begin{bmatrix} \mSigma_{\eta\eta} & \mSigma_{\eta \epsilon} \\ \mSigma_{\epsilon \eta} & \mSigma_{\epsilon \epsilon} \end{bmatrix}
\begin{bmatrix}
\bm{\mathcal{A}}(1)' & \mZeros \\
\mZeros & \bm{\mathcal{D}}(1)'
\end{bmatrix}^{-1}.
\label{eq:covariancetransform}
\end{equation}
Apart from this long-run covariance matrix $\mOmega=\sum_{h=-\infty}^\infty \E\big(\vxi_t^{}\vxi_{t+h}'\big)$, we also introduce the one-sided long-run covariance matrix $\mDelta= \left[\begin{smallmatrix} \mDelta_{uu} & \mDelta_{uv} \\ \mDelta_{vu} & \mDelta_{vv} \end{smallmatrix} \right]= \sum_{h=0}^\infty \E\big(\vxi_t^{}\vxi_{t+h}'\big)$. The Brownian motion defined by $\bm{B}_{u.v}=\bm{B}_u-\mOmega_{uv}^{} \mOmega_{vv}^{-1} \bm{B}_v^{}$ is by construction orthogonal to $\bm{B}_v$. Its $(n\times n)$ covariance matrix equals $\mOmega_{u.v}^{}=\mOmega_{uu}^{}-\mOmega_{uv}^{} \mOmega_{vv}^{-1} \mOmega_{vu}^{}$.
\subsection{Infeasible GLS}
We start our analysis assuming that the $(nT\times nT)$ covariance matrix $\mSigma_{\vu}(q)$ is a known quantity for each $q$. The modified Cholesky block decomposition of page \pageref{abr:MCBD} can now be used to derive the limiting distribution of this infeasible GLS estimator. A insightful exposition of our results requires further notation.
\begin{enumerate}[(a)]
\item Introduce scaling matrices: $\mG_{\vd_i,T}:=T^{-1/2}\diag[1,T^{-1},\ldots,T^{-d_i}]$ for the time trends, and $\mG_{\vs_i,T}:=T^{-1/2}\diag[T^{-1/2},T^{-1},\dots,T^{-s_i/2}]$ for the stochastic trends. Moreover, we define $\mG_T:=\diag\left[\mG_{1,T},\dots,\mG_{n,T}\right]$, where $\mG_{i,T}:=\diag\left[\mG_{\vd_i,T},\mG_{\vs_i,T}\right]$.
\item Let $\vd_i(r):=\left[1,r,\dots,r^{d_i}\right]'$, $\bm{B}_{s_i}(r):=\left[B_{v_i}^{}(r),B_{v_i}^2(r),\dots,B_{v_i}^{s_i}(r)\right]'$ and $\vj_i(r):=\left[\vd_i(r)',\bm{B}_{s_i}(r)'\right]'$. Define $d\times n$ block-diagonal random matrix $\mJ(r):=\diag\left[\vj_1(r),\dots,\vj_n(r)\right]$.
\item $\vb_i=\Big[\vzeros_{d_i+1}',1,2\int_{0}^{1}B_{v_i}^{}(r)dr,\dots,s_i\int_{0}^{1}B_{v_i}^{s_i-1}(r)dr\Big]'$.
\end{enumerate}
Finally, we use $\bm{B}_{v_j}$ as shorthand notation for the $j^{th}$ component of $\bm{B}_v$.
\begin{theorem}[Limiting Distribution of the infeasible GLS Estimator]\label{thm:infeasGLS}
If Assumption \ref{assumpt1:linearproc} holds, and if $q=q(T)$ satisfies $\frac{1}{q}+\frac{q}{T}\to 0$ as $T\to\infty$, then
\begin{equation}
\begin{aligned}
\mG_T^{-1}\left(\widehat{\vbeta}_{GLS}-\vbeta\right) &\wto\left(\int_0^1 \mJ(r)\mOmega_{uu}^{-1}\mJ(r)'dr\right)^{-1}\\
&\qquad\times\left(\int_0^1\mJ(r)\mOmega_{uu}^{-1}d\bm{B}_{u.v}^{}(r)+\int_0^1 \mJ(r)\mOmega_{uu}^{-1}\mOmega_{uv}^{}\mOmega_{vv}^{-1}d\bm{B}_{v}^{}(r)+{\bm{\mathcal{B}}}\right),
\end{aligned}
\label{eq:infeasibleGLSlimiting}
\end{equation}
where ${\bm{\mathcal{B}}}=\left[{\bm{\mathcal{B}}}_{1}',\dots,{\bm{\mathcal{B}}}_{n}'\right]'$ and ${\bm{\mathcal{B}}}_{i}=\operatorname{row}_i\Big(\mSigma_{\epsilon\eta}\Big)\col_i\Big(\mSigma_{\eta\eta}^{-1}\Big)~\vb_i$.
\end{theorem}
The limiting result in \eqref{eq:infeasibleGLSlimiting} coincides with the limiting distribution of the MSUR estimator, $\widetilde{\vbeta}_{MSUR}:= \left(\mZ'(\mI_T\otimes \widehat{\mOmega}_{uu}^{-1}) \mZ \right)^{-1} \Big( \mZ'(\mI_T\otimes \widehat{\mOmega}_{uu}^{-1}) \vy\Big)$, as reported in \cite{wagnergrabarczykhong2019}, see their Proof of Proposition 1. The equivalence of these limiting distributions is caused by the facts that: (1) applying a linear filter to an integrated series only affect its long-run variance (e.g. \cite{phillipspark1988}), and (2) the previous statement remaining true when applying a linear filter to higher integer powers of integrated series. The terms $\int_0^1 \mJ(r)\mOmega_{uu}^{-1}\mOmega_{uv}^{}\mOmega_{vv}^{-1}d\bm{B}_{v}^{}(r)$ and ${\bm{\mathcal{B}}}$ in \eqref{eq:infeasibleGLSlimiting} reflect the presence of second order bias terms caused by serial correlation and endogeneity. In Section \ref{sec:fullymodifiedinf}, we introduce the fully modified (FM) correction that adjust these bias terms and leads to standard inference. We first introduce a feasible version of the GLS estimator.
\subsection{Consistent Estimation of $\mSigma_{\vu}^{-1}(q)$ and Feasible GLS}\label{subsec:covarianceconsistency}
Up to this point we have discussed the infeasible estimator $\widehat{\vbeta}_{GLS}:=\left(\mZ' \mSigma_{\vu}^{-1}(q) \mZ \right)^{-1} \mZ' \mSigma_{\vu}^{-1}(q) \vy$. A feasible GLS approach requires a consistent estimator of the $(nT\times nT)$ matrix $\mSigma_{\vu}^{-1}(q)$. Several authors, e.g. \cite{wupourahmadi2009} and \cite{mcmurypolitis2010}, have constructed consistent estimators of large covariance matrices using banding or tapering to reduce the number of unknown parameters. Direct usage of their results poses two difficulties because: (1) numerical inversion of large matrices is computationally expensive for large $nT$, and (2) matrix inversion might even be impossible because the estimated covariance matrix cannot be guaranteed to be positive definite. In the light of the such considerations, we will estimate $\mSigma_{\vu}^{-1}(q)$ directly and ensure it to be positive definite. The approach is the sample counterpart of the BIAM described on page \pageref{ref:BIAM}. That is, we replace true innovations by first stage OLS residuals $\widehat{\vu}_t=\vy_t-\mZ_t\widehat{\vbeta}_{OLS}$, and subsequently minimise a sample moment in estimated residuals rather than the population mean squared forecasting error. This method was previously used by \cite{chengingyu2015} and \cite{ICG2016} for univariate time series. For a multivariate time series, we define
\begin{equation}
\begin{aligned}
\widehat{\mA}(\ell)&=
\begin{bmatrix}
\widehat{\mA}_1(\ell) & \cdots & \widehat{\mA}_\ell(\ell)
\end{bmatrix}
=
\operatornamewithlimits{arg\;min}_{(\mTheta_1,\ldots,\mTheta_\ell)\in \SR^{n\times n\ell}}\sum_{t=\ell+1}^{T}\left\|\widehat{\vu}_t-\mTheta_1 \widehat{\vu}_{t-1}-\cdots-\mTheta_\ell \widehat{\vu}_{t-\ell} \right\|^2\\
\widehat{\mS}(\ell)&=\frac{1}{T-\ell}\sum_{t=\ell+1}^{T}\left[\widehat{\vu}_t-\widehat{\mA}_1(\ell) \widehat{\vu}_{t-1}-\cdots-\widehat{\mA}_\ell(\ell)\widehat{\vu}_{t-\ell}\right]\left[\widehat{\vu}_t-\widehat{\mA}_1(\ell) \widehat{\vu}_{t-1}-\cdots-\widehat{\mA}_\ell(\ell)\widehat{\vu}_{t-\ell}\right]',
\end{aligned}
\label{eq:sampleMCD}
\end{equation}
$1\leq \ell\leq q$, and $\widehat{\mS}(0)=\frac{1}{T}\sum_{t=1}^{T}\widehat{\vu}_t^{}\widehat{\vu}_t'$. Similarly to \eqref{eq:bimam}-\eqref{eq:Mmatrixconstrucion}, we subsequently construct the matrices $\widehat{\bm{\mathcal{M}}}_{\vu}(q)=\left[\widehat{\vm}_{\vu}^{ij}(q) \right]_{1\leq i,j\leq T}$ and $\widehat{\bm{\mathcal{S}}}_{\vu}(q)=\diag\left(\widehat{\mS}(0),\widehat{\mS}(1),\ldots,\widehat{\mS}(q),\ldots,\widehat{\mS}(q)\right)$, and obtain the implied multivariate BIAM estimator as
\begin{equation}
\widehat{\mSigma_{\vu}^{-1}}(q)=\widehat{\bm{\mathcal{M}}}_{\vu}'(q)\widehat{\bm{\mathcal{S}}}_{\vu}^{-1}(q)\widehat{\bm{\mathcal{M}}}_{\vu}^{}(q).
\label{eq:BIAM}
\end{equation}
\begin{assumption}\label{assump3:residuals}
For $\widehat{\vu}=[\widehat{\vu}_1',\ldots,\widehat{\vu}_T']'$ and $\vu=[\vu_1',\ldots,\vu_T']'$, assume $\|\widehat{\vu}-\vu\|^{2}=O_p(1)$.
\end{assumption}
\begin{assumption}\label{assump4:bandingparameter}
Assume $q=q_T$ satisfies $\frac{1}{q_T}+\frac{q_T^3}{T}\to 0$ as $T\to\infty$.
\end{assumption}
Assumption \ref{assump3:residuals} requires the residuals to be sufficiently close to the true innovations. It is a rather mild assumption and it is satisfied if residuals are computed by least squares. Assumption \ref{assump4:bandingparameter} places constraints on the banding parameter $q_T$. First, Assumption \ref{assump4:bandingparameter} requires the banding parameter to diverge with sample size. This ensures that no nonzero elements are (asymptotically) set to zero. Moreover, the assumption $q_T^3/T\to 0$ establishes an upper bound for the growth rate of $q_T$. The definition of $\widehat{\mA}(\ell)$, see \eqref{eq:sampleMCD}, shows that we are fitting a vector autoregression (VAR) of increasing lag order to the residuals. Identical rate requirements are reported by \cite{lewisreinsel1985} when they derive consistency and asymptotic normality results when finite VAR models are fitted to infinite order VAR processes. The following theorem shows the consistent estimation of $\mSigma_{\vu}^{-1}$ and implies that the infeasible and feasible GLS estimator have the same limiting distribution.
\begin{theorem}[Consistent Estimation of $\mSigma_{\vu}^{-1}$]\label{thm:consistentDECOMP}
If Assumptions \ref{assumpt1:linearproc}-\ref{assump4:bandingparameter} hold, then
\begin{equation}
\begin{aligned}
\left\|\widehat{\mSigma_{\vu}^{-1}}(q_T)-\mSigma_{\vu}^{-1}\right\|
&\leq\Big\|\widehat{\mSigma_{\vu}^{-1}}(q_T)-\mSigma_{\vu}^{-1}(q_T)\Big\|+\Big\|\mSigma_{\vu}^{-1}(q_T)-\mSigma_{\vu}^{-1}\Big\| \\
&=O_p\left(\,\sqrt{q_T^3/T} \, \right)+O\left(\frac{1}{\sqrt{q_T}}\sum_{s=q_T+1}^{\infty}s\left\|\mA_{s}\right\|_{\calF}\right) \pto 0 \text{ as $T\rightarrow \infty$}.
\end{aligned}
\end{equation}
\end{theorem}
\subsection{Fully Modified Inference} \label{sec:fullymodifiedinf}
The asymptotic results of Theorem \ref{thm:infeasGLS} is not immediately useful for statistical inference. There are two difficulties. First, the second order bias dislocates the limiting distribution which can translate into substantial finite sample bias. This leads to a loss in efficiency. Second, possible dependencies between the Brownian motions $\bm{B}_u$ and $\bm{B}_v$ cause the limiting distribution to depend on nuisance parameters. Critical values would therefore be nuisance parameter dependent as well.
These two issues have received extensive attention in the linear cointegration literature. A (non-exhaustive) list of solution methods is: joint modeling as in \cite{phillips1991}, \citeauthor{saikkonen1992}'s\ (\citeyear{saikkonen1992}) dynamic least squares, and the integrated modified OLS and fixed-b approaches by \cite{vogelsangwagner2014}. We rely on the fully modified (FM) approach advocated by \cite{phillipshansen1990} and \cite{phillips1995}. The idea is a twofold modification of the estimator: (1) second order bias terms are subtracted, and (2) a transformation of the dependent variable is introduced to obtain a zero-mean Gaussian mixture limiting distribution. Recently, \cite{wagnergrabarczykhong2019} have proposed two estimators within the framework of seemingly unrelated cointegrating polynomial regressions. These estimators, FM-SOLS and FM-SUR, rely on kernel estimators of the one- and two-sided long-run covariance matrix (see Theorem \ref{thm:fullymodifiedOLSGLS}). As such, we introduce the following assumption.
\begin{assumption}[Consistent Estimation of Long-run Covariance Matrices]\label{assumpt:consistentlongrun}
$\widehat{\mOmega}$ and $\widehat{\mDelta}$ are consistent kernel estimators of the long-run covariance matrix $\mOmega$ and the one-sided long-run covariance matrix $\mDelta$, respectively.
\end{assumption}
\cite{andrews1991} and \cite{neweywest1994} use kernel estimators for long-run covariance estimation. Their method involves the calculation of weighted sums of the autocovariance matrices of the residuals. These weights are determined by a kernel function and bandwidth parameter. Our Assumption \ref{assumpt:consistentlongrun} is easily satisfied by imposing suitable conditions on the kernel function and bandwidth parameter. We refer to \cite{phillips1995} and \cite{jansson2002} for an enumeration of such conditions.
Alternatively, we can obtain consistent one- and two-sided long-run covariance estimators within the BIAM framework of Section \ref{subsec:covarianceconsistency}.\footnote{An overview of the procedure is given here. Section \ref{detailsFMinference} in the Supplement provides further details.} This approach resembles \cite{berk1974}. The GLS estimator and its FM counterpart are thus constructed within a single framework. The estimators are as follows. For all $t=1,2,\ldots,T$, we first stack $\widehat{\vu}_t$ and $\Delta \vx_t = \vv_t$ in the $2n$-dimensional vector $\widehat{\vxi}_t=[\widehat{\vu}_t',\Delta \vx_t']'$. Since the BIAM estimator is fitting VAR processes up to order $q_T$, we will use the estimated VAR($q_T$) approximations to define the long-run covariance estimators. For $\mOmega$, the estimator is $\widehat{\mOmega}_{q_T}=\left(\mI_{2n}-\sum_{j=1}^{q_T} \widehat{\mF}_j^{}(q_T) \right)^{-1} \widehat{\mSigma}_{q_T}\left(\mI_{2n}-\sum_{j=1}^{q_T} \widehat{\mF}_j'(q_T) \right)^{-1}\label{eq:twosidedLRV}$, where $\widehat{\mSigma}_{q_T}=\widehat{\mS}(q_T)$ and $\widehat{\mF}_j(q_T)$ denote respectively the estimated prediction error variance and the coefficient matrix of the $j^{th}$ lag when a VAR($q_T$) is fitted to $\{\widehat{\vxi}_t \}_{t=1}^T$. The population one-sided long-run covariance matrix is $\mDelta= \sum_{h=0}^\infty \E\big(\vxi_t^{}\vxi_{t+h}'\big)$. It is thus intuitive to approximate this quantity by a finite sum of estimated covariance matrices of $\{\widehat{\vxi}_t \}_{t=1}^T$. These covariance matrices are nothing but subblocks of the matrix $\widehat{\mSigma}_{\vxi}^{}(q_T)=\widehat{\bm{\mathcal{M}}}_{\vxi}^{-1}(q_T)\widehat{\bm{\mathcal{S}}}_{\vxi}^{}(q_T)\widehat{\bm{\mathcal{M}}}_{\vxi}^{-1\prime}(q_T)$.\footnote{We use $\mSigma_{\vxi}$ to denote the $(2nT\times 2nT)$ matrix $\E(\vxi \vxi')$ where $\vxi=[\vxi_1',\vxi_2',\ldots,\vxi_T']'$. The matrices $\widehat{\bm{\mathcal{M}}}_{\vxi}(q)$ and $\widehat{\bm{\mathcal{S}}}_{\vxi}(q)$ are defined similarly to respectively $\widehat{\bm{\mathcal{M}}}_{\vu}(q)$ and $\widehat{\bm{\mathcal{S}}}_{\vu}(q)$ (see page \pageref{eq:BIAM}). The matrix $\widehat{\bm{\mathcal{M}}}_{\vxi}(q_T)$ is lower triangular with identity matrices on the main diagonal. Therefore, its matrix inverse exists and is fast to compute.} We therefore use
\begin{equation}
\widehat{\mDelta}_{q_T,r_T}= \mQ_{r_T}' \widehat{\mSigma}_{\vxi}^{}(q_T) \mQ_1^{}
\end{equation}
where $\mQ_r=\left[\mZeros_{2n \times 2n}, \cdots, \mZeros_{2n \times 2n}, \mI_{2n}, \cdots, \mI_{2n}\right]'$ is an $\big(2n T \times 2n \big)$ block matrix of zeros of which the last $r$ blocks have been replaced by identity matrices. To ensure consistency, we place the following rate restriction on the number of included autocovariance matrices.
\begin{assumption} \label{assumpt:rT}
As $T\rightarrow\infty$, $r_T\rightarrow\infty$, $\frac{r_T^{}q_T^3}{T}\rightarrow 0$, and $r_T=O(q_T)$.
\end{assumption}
Definitions and limiting results for FM estimators are presented in Theorem \ref{thm:fullymodifiedOLSGLS}. The FM-SOLS, FM-SUR and FM-GLS estimator all depend on estimators for $\mDelta$ and $\mOmega$. It is only the consistency of these estimators that is relevant for the asymptotic analysis, not whether the kernel or BIAM approach is employed. As such, we will not complicate notation by introducing additional notation to indicate whether the kernel or BIAM approach is used. In subsequent theorems, simulation results and the empirical application we will use kernel estimators for FM-SOLS and FM-SUR, and the BIAM approach for FM-GLS. This seems to be the logical choice for these estimators.
\begin{theorem}\label{thm:fullymodifiedOLSGLS}
For $i=1,\ldots,n$, define $\widehat{\vb}_i=\left[\vzeros_{d_i+1}',T, 2 \sum_{t=1}^T x_{it},\ldots, s_i \sum_{t=1}^T x_{it}^{s_i-1} \right]'$, for $i=1,\ldots,n$. Also, define the $(n\times n)$ matrix $\widehat{\mDelta}_{vu}^+$ as the (implied) consistent estimator of $\mDelta_{vu}^+=\mDelta_{vu}^{}-\mDelta_{vv}^{}\mOmega_{vv}^{-1} \mOmega_{vu}^{}$.
\begin{enumerate}[(a)]
\item Define the FM-SOLS estimator as
\begin{equation}
\widehat{\vbeta}_{SOLS}^+ = \left( \mZ' \mZ \right)^{-1} \left( \mZ' \vy^+ -\widehat{\mA} \, \right),
\end{equation}
where $\vy^+:=[\vy_{1}^{+\prime},\vy_{2}^{+\prime},\ldots,\vy_{T}^{+\prime}]'$ with $\vy_t^+ = \vy_t^{}-\widehat{\mOmega}_{uv}^{} \widehat{\mOmega}_{vv}^{-1} \diff \vx_t^{}$, and $\widehat{\mA}:= [\widehat{\mA}_1',\ldots,\widehat{\mA}_n']'$ with $\widehat{\mA}_i= \widehat{\mDelta}_{v_i u_i}^+\widehat{\vb}_i$ and $\widehat{\mDelta}_{v_i u_i}^+$ being the $i^{th}$ element on the main diagonal of $\widehat{\mDelta}_{v u}^+$. If Assumptions \ref{assumpt1:linearproc} and \ref{assumpt:consistentlongrun} hold, then
\begin{equation}
\mG_T^{-1} \left( \widehat{\vbeta}_{SOLS}^+ - \vbeta \right) \wto \left( \int_0^1 \mJ(r) \mJ(r)' dr \right)^{-1} \int_0^1 \mJ(r) d\bm{B}_{u.v}(r).
\end{equation}
\item Define the FM-SUR estimator as
\begin{equation}
\widehat{\vbeta}_{SUR}^+ = \left(\mZ' (\mI_T\otimes \widehat{\mOmega}_{u.v}^{-1}) \mZ \right)^{-1} \left( \mZ' (\mI_T\otimes \widehat{\mOmega}_{u.v}^{-1}) \vy^+ - \widetilde{\mA}^* \right),
\end{equation}
where $\widetilde{\mA}^*:= [\widetilde{\mA}_1^*,\ldots,\widetilde{\mA}_n^*]$ with $\widetilde{\mA}_i^*= \operatorname{row}_i\left(\widehat{\mDelta}_{vu}^+\right) \col_i\left( \widehat{\mOmega}_{u.v}^{-1} \right) \widehat{\vb}_i$. If Assumptions \ref{assumpt1:linearproc} and \ref{assumpt:consistentlongrun} hold, then
\begin{equation}
\mG_T^{-1} \left( \widehat{\vbeta}_{SUR}^+ - \vbeta \right) \wto \left( \int_0^1 \mJ(r) \mOmega_{u.v}^{-1} \mJ(r)' dr \right)^{-1} \int_0^1 \mJ(r) \mOmega_{u.v}^{-1} d\bm{B}_{u.v}^{}(r).
\end{equation}
\item Define the FM-GLS estimator as
\begin{equation}
\widehat{\vbeta}_{FGLS}^+=\Big(\mZ'\widehat{\mSigma_{\vu}^{-1}}(q)\mZ\Big)^{-1}\left[\mZ'\widehat{\mSigma_{\vu}^{-1}}(q)\vy-\mZ'\Big(\mI_T\otimes \widehat{\mOmega}_{uu}^{-1}\widehat{\mOmega}_{uv}^{}\widehat{\mOmega}_{vv}^{-1}\Big)\vv-\widehat{{\bm{\mathcal{B}}}}^{+}\right],
\label{eq:FMGLSestimator}
\end{equation}
where $\vv := [\diff \vx_1',\ldots,\diff \vx_T']' = [\vv_1',\ldots,\vv_T']'$, and $\widehat{{\bm{\mathcal{B}}}}^{+}=\big[\widehat{{\bm{\mathcal{B}}}}_1^{+\prime},\dots,\widehat{{\bm{\mathcal{B}}}}_n^{+\prime}\big]'$ with
$$\widehat{{\bm{\mathcal{B}}}}_{i}^{+}=\left[\operatorname{row}_i\Big(\widehat{\mSigma}_{\epsilon\eta}\Big)\col_i\Big(\widehat{\mSigma}_{\eta\eta}^{-1}\Big)-\operatorname{row}_i\Big(\widehat{\mDelta}_{vv}\Big)\col_i\Big(\widehat{\mOmega}_{vv}^{-1}\widehat{\mOmega}_{vu}^{}\widehat{\mOmega}_{uu}^{-1}\Big)\right]~\widehat{\vb}_i.
$$
If Assumptions \ref{assumpt1:linearproc}-\ref{assump4:bandingparameter} and \ref{assumpt:rT} hold, then
\begin{equation}
\mG_T^{-1} \left( \widehat{\vbeta}_{FGLS}^+ - \vbeta \right) \wto \left(\int_0^1 \mJ(r)\mOmega_{uu}^{-1}\mJ(r)'dr\right)^{-1} \int_0^1 \mJ(r)\mOmega_{uu}^{-1}d\bm{B}_{u.v}^{}(r).
\end{equation}
\end{enumerate}
\end{theorem}
The FM-GLS estimator is new to the seemingly unrelated CPR literature, whereas the FM-SOLS and FM-SUR estimators have recently appeared in \cite{wagnergrabarczykhong2019}. Theorem \ref{thm:fullymodifiedOLSGLS} indicates that all three estimators have a zero-mean Gaussian mixture limiting distribution implying that standard inference is applicable for each. However, we also see from Theorem \ref{thm:fullymodifiedOLSGLS} that the limiting distributions are generally different because different types of weighing are used in the construction of the estimators.\footnote{There are special cases in which some (pairs of) estimators become asymptotically equivalent. For example, if $n=1$, then all estimators are asymptotically equivalent because the weighting matrices $\mOmega_{u.v}^{-1}$ and $\mOmega_{uu}^{-1}$ are now scalars. Also, under exogeneity, we have $\mOmega_{uu} = \mOmega_{u.v}$ and the FM-SUR an FM-GLS estimators share the same limiting distribution.}
For completeness, we also detail how the FM-GLS estimator can be used to test linear hypotheses. A formal presentation of such a result is more involved because of the different convergence rates of the individual parameter estimators. That is, the parameters with the lowest convergence rate will dominate the asymptotic distribution and one should take care to avoid a degenerate limiting distribution. We will rule out such complications by considering hypothesis tests on individual parameters.\footnote{For general linear hypothesis, we refer the reader to \cite{simsstockwatson1990} where a reordering based on convergence rates is used to establish the limiting distribution of the Wald $F$ statistic for general linear hypothesis. The same approach is applicable in our setting but we will not explore this in greater detail.} Therefore, let $\mR$ denote a $(k \times s)$ selection matrix in which every row contains a single 1 and zeros otherwise. The null hypothesis $\mR \vbeta= \vr$ can be tested using the standard chi-squared limiting distribution of the Wald statistic (Theorem \ref{thm:test}). These tests are practically relevant. For example, exclusion restrictions of the type $\mR \vbeta=\vzeros$ allow us to test whether the cointegrating relation is linear.
\begin{theorem}\label{thm:test}
Consider the null hypothesis $H_0:\mR\vbeta=\vr$, which imposes $k$ linearly independent restrictions. Under the assumptions of Theorem \ref{thm:fullymodifiedOLSGLS}(c), the Wald test statistic
\begin{equation}
\calW=\Big(\mR\widehat{\vbeta}_{FGLS}^+-\vr\Big)'\widehat{\mPhi}^{-1}\Big(\mR\widehat{\vbeta}_{FGLS}^+-\vr\Big)\wto \chi_k^2,
\end{equation}
where $\widehat{\mPhi}=\mR\left[\mZ' \left( \mI_T \otimes \widehat{\mOmega}_{uu}^{-1} \right) \mZ \right]^{-1}\left[\mZ'\Big(\mI_T\otimes \widehat{\mOmega}_{uu}^{-1}\widehat{\mOmega}_{u.v}^{}\widehat{\mOmega}_{uu}^{-1}\Big)\mZ\right]\left[\mZ' \left( \mI_T \otimes \widehat{\mOmega}_{uu}^{-1} \right) \mZ \right]^{-1}\mR'$.
\end{theorem}
\subsection{Testing the Null of Cointegration} \label{sec:cointtest}
Stationarity tests are used to avoid spurious regressions and to verify the correct specification of the cointegrating relation. To test for stationarity of the seemingly unrelated cointegrating polynomial regressions (SUCPR) errors, we combine the test statistic from \cite{nyblomharvey2000} with the sub-sampling approach found in \cite{choisaikkonen2010} and \cite{wagnerhong2016}. We consider three test statistics. To treat all test statistics in a unified framework, we define
\begin{equation}
\bm \varphi_{j,b}(\{\vx\}) = \left[ \vx_j' , \sum_{s=j}^{j+1} \vx_s',\ldots, \sum_{s=j}^{j+b-1} \vx_s' \right]',
\end{equation}
that is, a vector of length $nb$ stacking the cumulative sums of $\{\vx_j,\ldots, \vx_{j+b-1} \}$. If the true innovations $\vu_1,\ldots,\vu_T$ were observed, then we could use the full-sample KPSS-type of test statistic $\frac{1}{T^2} \bm \varphi_{1,T}(\{\vu\})'(\mI_T\otimes \widehat{\mOmega}_{uu}^{-1})\bm \varphi_{1,T}(\{\vu\})=\tr\left[ \widehat{\mOmega}_{uu}^{-1} \frac{1}{T^2} \sum_{t=1}^T \left( \sum_{s=1}^t \vu_s \right)\left( \sum_{s=1}^t \vu_s \right) ' \right]$ to test for stationarity of the innovations. Under the null of stationarity, this test statistic would converge weakly to $\int_0^1 \| \bm W(r) \|^2 dr$ with $\bm W(r)$ denoting an $n$-dimensional standard Brownian motion. This limiting distribution is free of nuisance parameters and the cumulative distribution function is available as a series expansion (see the Supplement).
The innovations $\vu_1,\ldots,\vu_T$ are only available when cointegrating relations are pre-specified. If these coefficients are estimated, then this additional parameter uncertainty will contaminate the limiting distribution with nuisance parameters.\footnote{There are exceptions. \cite{shin1994} reports a nuisance parameter free limiting distribution for a single-equation linear cointegrating relation. This remains true if only a single integrated variable enters the cointegrating regression with a higher power, see Proposition 5 in \cite{wagnerhong2016}.} The idea behind the subsampling approach is to construct a test statistic incorporating $b = b_T$ residuals while computing parameter estimators from all $T$ observations. If $b_T$ increases slowly with sample size, then the parameter estimation error will be negligible relative to the randomness in the errors and the asymptotic distribution remains $\int_0^1 \| \bm W(r) \|^2 dr$.
The three KPSS-type of test are based on the following residuals: $\hat{\vu}_{t,SOLS}^+= \vy_t^+- \mZ_t \widehat{\vbeta}_{SOLS}^+$, $\hat{\vu}_{t,SUR}^+= \vy_t^+ - \mZ_t \widehat{\vbeta}_{SUR}^+$, and $\hat{\vu}_{t,FGLS}=\vy_t- \mZ_t \widehat{\vbeta}_{FGLS}^+$. The test statistic are:
\begin{equation}
K_{j,b_T}^i = \frac{1}{b_T^2} \bm \varphi_{j,b_T}(\{\hat{\vu}_i^+\})' \left(\mI_{b_T} \otimes \widehat{\mOmega}_{u.v}^{-1} \right)\bm \varphi_{j,b_T}(\{\hat{\vu}_i^+\}),\qquad\qquad \text{for }i\in\{SOLS, SUR\},
\label{eq:KPSSresidualtype}
\end{equation}
and
\begin{equation}
K_{j,b_T}^{BIAM} = \frac{1}{b_T^2} \bm \varphi_{j,b_T}(\{\hat{\vu}_{FGLS}\})' \widehat{\mSigma_{\vu}^{-1}}(q_T,b_T) \bm \varphi_{j,b_T}(\{\hat{\vu}_{FGLS}\}),
\label{eq:KPSSbiam}
\end{equation}
where $\widehat{\mSigma_{\vu}^{-1}}(q_T,b_T)$ is the $(n b_T\times n b_T)$ submatrix of $\widehat{\mSigma_{\vu}^{-1}}(q_T)$ obtained by selecting the rows and columns related to all time indices in the set $\{ n(T-b_T)+1, n(T-b_T)+2 ,\ldots,nT\}$. The test statistic in \eqref{eq:KPSSbiam} fits naturally into the FM-GLS estimation framework.
\begin{theorem}\label{thm:kpss_subtest}
Let the assumptions from Theorem \ref{thm:fullymodifiedOLSGLS} hold.
\begin{enumerate}[(a)]
\item If $\frac{1}{b_T}+\frac{b_T}{T}\to 0$ as $T\to \infty$, then
$$
K_{j,b_T}^i \wto \int_0^1 \| \bm{W}(r) \|^2 dr, \qquad 1 \leq j \leq T-b_T+1, \qquad \text{for }i\in\{SOLS, SUR \}.
$$
\item If $\frac{q_T}{b_T}+\frac{b_T}{T}\to 0$ as $T\to \infty$, then $ K_{j,b_T}^{BIAM}\wto \int_0^1 \| \bm{W}(r) \|^2 dr$ for any $1 \leq j \leq T-b_T+1$.
\end{enumerate}
\end{theorem}
A sample of size $T$ allows for up to $M=\lfloor T/b_T\rfloor$ series of nonoverlapping blocks of residuals of length $b_T$. Similarly to \cite{choisaikkonen2010}, we apply the Bonferroni procedure to use all these series and thereby increase power. The approach is applicable to any of the three test statistics in Theorem \ref{thm:kpss_subtest}. As such, we keep the notation general and use a generic $K_j$ to denote a test statistic based on the $j^{th}$ subseries, $j=1,2,\ldots,M$. In the Bonferrroni procedure we compute $K_{max}=\left\{K_1,K_2,\ldots,K_M\right\}$ and do not reject the null hypothesis whenever $K_{max}\leq c_{\alpha/M}$ with $c_{\alpha/M}$ defined by $\mathbb{P}\left(\int_0^1 \| \bm{W}(r) \|^2 dr\geq c_{\alpha/M}\right)=\alpha/M$. The Bonferroni inequality implies $\lim_{T\to \infty} \mathbb{P}\left(K_{max} \leq c_{\alpha/M} \right) \geq 1- \lim_{T\to\infty} \sum_{j=1}^M\mathbb{P}\left(K_j>c_{\alpha/M} \right)=1- \alpha$ and we see that the probability of a type-I error does not exceed the significance level $\alpha$.
\begin{remark}
We suggest to follow \cite{choisaikkonen2010} in terms of the implementation of the subsampling approach. That is, the block size $b_T$ is selected using the minimum volatility rule by \cite{romanowolf2001}. For this particular block size we subsequently select subsamples by taking non-overlapping blocks from alternatively the start and the end of the sample.
\end{remark}
\begin{remark}
The limiting results in Theorem \ref{thm:kpss_subtest} guarantee a correct asymptotic size. Our simulations show (1) that these tests have power against various alternative hypotheses and (2) that power increases with sample size. A theoretical investigation of the power properties is outside of the scope of this paper.
\end{remark}
\section{Simulations} \label{sec:simulations}
We now study the finite sample performance of the estimators and stationarity tests. First, we compare the FM-GLS estimator with the FM-SOLS and FM-SUR estimators from \cite{wagnergrabarczykhong2019}. All long-run covariance matrices are computed using a Bartlett kernel and the automatic bandwidth selection approach due to \cite{andrews1991}. For FM-GLS, the banding parameter $q_T$ is selected using the subsampling and risk-minimization approach explained in section 5 from \cite{bickellevina2008}.\footnote{More details concerning the implementation can be found in the Supplement.} Infeasible counterparts of the estimator are constructed assuming the knowledge of the true serial correlation and/or cross-sectional dependence pattern. These estimators are denoted by infSOLS, infSUR, and infGLS. Second, we look at the cointegration tests. We consider three test statistics: $K^{SOLS}$ and $K^{SUR}$ use the residuals as in \eqref{eq:KPSSresidualtype}, whereas $K^{BIAM}$ employs the pre-filtered residuals from \eqref{eq:KPSSbiam}. All tests are implemented with minimum volatility block size selection and Bonferroni correction.
We consider $T\in\{100,200,500\}$ and $n\in\{3,5\}$. All tests are performed at a nominal significant level of $5\%$. For stationary processes, a presample of 200 observations is used to remove the influence of the starting values. All results are based on $2.5\times 10^4$ Monte Carlo replicates.
\subsection{Monte Carlo Designs}
We generate data according to a quadratic seemingly unrelated CPR. That is, we adopt the DGP in \eqref{eq:basemodel} with $\vz_{it}=\big[1,t,x_{it},x_{it}^2\big]'$. The integrated variables satisfy $\vx_0=\vzeros$ and $\diff \vx_t = \vv_t$. We explore two error processes.
\bigskip
\noindent\textbf{Setting A (Errors as in \cite{wagnergrabarczykhong2019})}: As a benchmark, we revisit the simulation setting in \cite{wagnergrabarczykhong2019} and generate innovations according to
\begin{equation}
\vu_t=\rho_1\vu_{t-1}+\vepsi_{t}+\rho_2\ve_t,\qquad \vv_t=\ve_t+0.5\ve_{t-1},
\end{equation}
where $\vepsi_t\stackrel{i.i.d.}{\sim}\rN\big(\vzeros,\mSigma(\rho_3)\big)$, $\ve_t\stackrel{i.i.d.}{\sim}\rN\big(\vzeros,\mSigma(\rho_4)\big)$ and
\begin{equation}\label{eq:toeplitz_structure}
\mSigma(\rho)=\begin{bmatrix}
1 & \rho & \cdots & \rho\\
\rho & 1 & \cdots & \rho\\
\vdots & \vdots & \ddots & \vdots\\
\rho & \rho & \cdots & 1\\
\end{bmatrix}
\end{equation}
is a symmetric Toeplitz matrix. The parameter $\rho_1$ controls the level of serial correlation and $\rho_2$ measures the degree of endogeneity. The parameters $\rho_3$ and $\rho_4$ indicate the extent of correlation across equations induced through $\vepsi_t$ and $\ve_t$, respectively. For simplicity, we assume identical values $\rho_1=\rho_2=\rho_3=\rho_4=\rho\in\{0,0.3,0.6,0.8\}$. The true coefficient vector is $\vbeta=\big[\vbeta_1',\dots,\vbeta_n'\big]'$, where $\vbeta_i=[1,1,5,\beta_{i,4}]'$ with $\beta_{i,4}=-0.3$, $i=1,\dots,n$.
\bigskip
\noindent\textbf{Setting B (VARMA Errors)}: To further investigate the importance of serial correlation, we consider a second specification of the innovation process:
\begin{equation}
\vu_t=\mLambda_1\vu_{t-1}+\veta_t+\mLambda_2\veta_{t-1},\qquad \vv_t=\mLambda_3\vv_{t-1}+\vepsi_{t},
\end{equation}
where $\veta_t$ and $\vepsi_{t}$ are generated as $\left[\begin{smallmatrix}
\veta_t\\
\vepsi_t
\end{smallmatrix}\right]\stackrel{i.i.d.}{\sim}
\rN\big(\vzeros,\mSigma(\theta)\big)$ and $\mSigma(\theta)\in \SR^{2n\times 2n}$ as in \eqref{eq:toeplitz_structure} but with parameter $\theta$. The matrices $\mLambda_i$ ($i=1,2,3$) are generated independently and similarly to \cite{chang2004}. That is, we take the following three steps:
\begin{enumerate}[(a)]
\item Generate an $n\times n$ random matrix $\mU_i$ from $\text{U}[0,1]$ and construct the orthogonal matrix $\mH_i=\mU_i^{}\Big(\mU_i'\mU_i^{} \Big)^{-1/2}$.
\item Generate $n$ eigenvalues $\lambda_{i1},\dots,\lambda_{in}\stackrel{i.i.d.}{\sim}\text{U}\big[\underline{\lambda},\bar{\lambda}\big]$.
\item Let $\mL_i=\diag\left(\lambda_{i1},\dots,\lambda_{in}\right)$ and compute $\mLambda_i=\mH_i^{}\mL_i^{}\mH_i'$.
\end{enumerate}
The parameter $\theta\in\{0.3,0.5\}$ governs regressor-error correlation and cross-equation correlation. The amount of serial correlation is specified through $\underline{\lambda}$ and $\bar{\lambda}$. The three scenarios $\big(\underline{\lambda},\bar{\lambda}\big)\in \big\{\left(0.1,0.5\right),\left(0.5,0.8\right),\left(0.8,0.95\right)\big\}$ steadily increase the autocorrelation in the generated data.
\bigskip
\noindent\textbf{Setting C (Cointegration Tests)}: We continue to construct innovations according to Setting B. Moreover, we fix $\left[\begin{smallmatrix}
\veta_t\\
\vepsi_t
\end{smallmatrix}\right]\stackrel{i.i.d.}{\sim}
\rN\big(\vzeros,\mSigma(\theta)\big)$ with $\theta=0.3$, and we construct the matrices $\mLambda_2$ and $\mLambda_3$ using $\big(\underline{\lambda},\bar{\lambda}\big)=(0.1,0.5)$. The eigenvalues of $\mLambda_1$ are varied to explore both size and power properties. We always estimate a \emph{quadratic} seemingly unrelated CPR.
\begin{description}
\item[Size DGP.] We generate the eigenvalues of $\mLambda_1$ as before. That is, take $\lambda_{11},\dots,\lambda_{1n}\stackrel{i.i.d.}{\sim}\text{U}\big[\underline{\lambda},\bar{\lambda}\big]$, where $\big(\underline{\lambda},\bar{\lambda}\big)\in \left\{\left(0.1,0.5\right),\left(0.5,0.8\right),\left(0.8,0.95\right)\right\}$.
\item[Power DGP1.] We set $\lambda_{1j}=1$ for $1\leq j\leq J_1$ and generate $\lambda_{1j}\stackrel{i.i.d.}{\sim}\text{U}\big[0.1,0.5\big]$ for $J_1+1\leq j\leq n$. The integer $J_1\in\{1,2,n\}$ represents the number of unit roots in $\{\vu_t\}$.
\item[Power DGP2.] The eigenvalues of $\mLambda_1$ are sampled as $\lambda_{11},\dots,\lambda_{1n}\stackrel{i.i.d.}{\sim}\text{U}\big[0.1,0.5\big]$, and the first $J_2\in\{1,2,n\}$ series follow a cubic SUCPR specification:
\begin{equation*}
y_{it}=\begin{cases}
1+t+5x_{it}-0.3x_{it}^2+0.01x_{it}^3+u_{it},&\quad 1\leq i\leq J_2,\\
1+t+5x_{it}-0.3x_{it}^2+u_{it},& \quad J_2+1\leq i\leq n.
\end{cases}
\end{equation*}
\item[Power DGP3.] We again take $\lambda_{11},\dots,\lambda_{1n}\stackrel{i.i.d.}{\sim}\text{U}\big[0.1,0.5\big]$ and construct
\begin{equation*}
y_{it}=\begin{cases}
\sum_{s=1}^{t} u_{is},&\quad 1\leq i\leq J_3,\\
1+t+5x_{it}-0.3x_{it}^2+u_{it},& \quad J_3+1\leq i\leq n,
\end{cases}
\end{equation*}
where $J_3\in\{1,2,n\}$ represents for the number of equations that specify a spurious relation.
\end{description}
Overall, the Power DGPs 1-3 consider: missing $I(1)$ regressors, omitted higher order powers of the $I(1)$ regressor $x_{it}$, and spurious regressions, respectively.
\subsection{Discussion of the Simulation Results}
Tables \ref{table:efficiencyQSUCPR} and \ref{tab:efficiency_varma} report the empirical mean squared error (MSE) for both feasible and infeasible estimators. As results are qualitatively similar across equations, we only report on the estimators for $\beta_{1,4}$ (the coefficient in front of $x_{1t}^2$). The column with FGLS contains the numerical value of the MSE and the MSEs of all other estimators are expressed relative to this benchmark. Values above 1 indicate a better performance of FM-GLS. We make the following observations:
\begin{enumerate}[(a)]
\item The FM-GLS estimator generally has the lowest MSE among all feasible estimators. These efficiency gains are small at low levels of endogeneity and serial correlation, but become sizeable at higher levels. Moreover, the Monte Carlo outcomes for the infeasible estimators indicate that these gains remain when the estimators are informed about the true endogeneity and serial correlation properties. It is thus the GLS weighting of the data that improves estimation accuracy.
\item There is one particular instance in Table \ref{tab:efficiency_varma} in which the performance of the FM-GLS estimator has a high MSE, namely the case of high persistency $\big(\underline{\lambda},\bar{\lambda}\big)=(0.8, 0.95)$, high endogeneity $\theta=0.5$, and small sample size $T=100$. This is caused by an inaccurate BIAM estimator resulting from the combination of a small sample size, high endogeneity, and high persistency. The problem disappears when $T$ increases.
\end{enumerate}
The subsequent set of simulations evolves around hypothesis testing, see Table \ref{table:sizeQSUCPR} and Figures \ref{fig:size_sc}-\ref{fig:power_jointtest_n5}. The errors are simulated using Setting A and we use the following Wald-type test statistics: the Wald-SOLS and Wald-SUR tests as developed in Proposition 2 in \cite{wagnergrabarczykhong2019}, and the Wald-FGLS test from Theorem \ref{thm:test}. We consider: (\textit{i}) the single equation test $H_0:\beta_{1,4}=-0.3$ against the two-sided alternative $H_1: \beta_{1,4}\neq -0.3$, and (\textit{ii}) the joint test $H_0: \beta_{1,4}=\beta_{2,4}=\ldots= \beta_{n,4}=-0.3$ against the alternative which rejects when at least one coefficient is unequal to $-0.3$. Some general remarks regarding size and size-corrected power are as follows.
\begin{enumerate}[(a)]
\setcounter{enumi}{2}
\item The Wald tests are typically oversized but the three tests react differently to changes in $\rho$. Increases in $\rho$ result in an increasing size for the SOLS and SUR version of the Wald test, whereas increases in $\rho$ lead to size decreases for the GLS type of Wald test. Overall, the GLS test provides better size control.
\item In Figures \ref{fig:size_sc} and \ref{fig:size_endo}, we vary the serial correlation parameter $\rho_1$ and the endogeneity parameter $\rho_2$ separately. Overall, variation in $\rho_1$ has a larger influence on size with the SUR test being most sensitive and the GLS test being least sensitive.
\item For all three Wald-type of tests, the size of the tests improves with sample size $T$.
\item The ordering in terms of size-corrected power is the same throughout Figures \ref{fig:power_singletest_n3}-\ref{fig:power_jointtest_n5}. That is, size-corrected power is lowest for Wald-SOLS, increases for Wald-SUR, and is highest for the Wald-FGLS test.
\end{enumerate}
The simulation results for the KPSS-type of cointegration tests can be found in Table \ref{tab:ct_tests}. The general conclusions are as follows.
\begin{enumerate}[(a)]
\setcounter{enumi}{6}
\item The empirical sizes of the $K^{SOLS}$ and $K^{SUR}$ tests are similar. We see: very conservative results for low serial correlation, decent size for medium serial correlation, and strongly oversized tests at high levels of serial correlation. These findings are completely in line with the simulation results that are reported in table 3 of \cite{choisaikkonen2010}. The same behaviour is observed for the $K^{BIAM}$ test but the deviations from the 5\% level are less extreme.
\item The power of the $K^{FOLS}$, $K^{SUR}$, and $K^{BIAM}$ tests behaves as expected: (1) power always increases with sample size, and (2) power increases when more unit roots, more misspecified equations, or more spurious relationships are incorporated in the DGP. The $K^{BIAM}$ test has the lowest power among the three tests. This is caused by the fact that the filter can nearly difference the data and hence make it appear more stationary.
\end{enumerate}
\section{Empirical Application}\label{sec:empappl}
The Environmental Kuznets Curve (EKC) conjectures an inverted U-shaped relation between environmental degradation and income per capita. That is, there is an initial decline in environmental quality with increasing economic activity, but beyond a certain turning point (caused by e.g. industrial transformation and increasing environmental awareness), economic growth goes hand in hand with environmental improvement. A more detailed description and historical overview of the EKC can be found in \cite{stern2004} and \cite{stern2017}, respectively. The implications of further economic growth on pollution, e.g. the emission of greenhouse-gases, are also key in understanding the future of global warming (\cite{nordhaus2013}).
We builds upon and compare with \cite{wagnergrabarczykhong2019}. That is, we look at carbon dioxide $(\text{CO}_2)$ emissions and GDP as proxies for environmental pollution and economic development (both per capita and in logarithms), respectively. The data is collected from the Maddison Project Database (MPD) and the homepage of the Carbon Dioxide Information Analysis Center (CDIAC).\footnote{The Maddison Project Database, \cite{madison2018}, contains the data on population size and real GDP. The data on $\text{CO}_2$ originates from \cite{cdiac2017}. We follow the official guidelines and multiply by $3.667$ and $10^6$ to convert the reported fossil-fuel emissions into total carbon dioxide emissions.} As in \cite{wagnergrabarczykhong2019}, we consider Austria (AT), Belgium (BE), Finland (FI), the Netherlands (NL), Switzerland (CH) and the United Kingdom (UK). Our yearly data spans the period from 1870 to 2014. We refer to the latter paper for a discussion of the stationarity properties of all series as well as the motivation for this particular set of countries. Overall, the dataset consist of $n=6$ countries with $T=145$ time series observations each. Such a panel with small $n$ and large $T$ is ideally suited for our FMGLS approach since the multivariate banded inverse autocovariance matrix remains computable.
We estimate the quadratic model specification:
\begin{equation}
e_{it}^{}=\beta_{i,1}^{}+\beta_{i,2}^{} t +\beta_{i,3}^{}g_{it}^{}+\beta_{i,4}^{} g_{it}^2+u_{it}^{},\quad i=1,2,\ldots,6,\quad t=1,2,\ldots,145,
\label{eq:ekc_quadraticmodel}
\end{equation}
where $e_{it}$ and $g_{it}$ are $\text{CO}_2$ emissions and GDP, respectively. As the first step in our analysis we employ the multivariate stationarity tests of Section \ref{sec:cointtest} to check this model specification (Table \ref{tab:multivariateKPSS}). All three tests reject the null of cointegration at a 5\% level signalling inappropriateness of the quadratic formulation. Figure \ref{fig:KPSSresiduals} shows the residuals on which these tests are based. What stands out in these graphs is the erratic behaviour of the series around the two world wars. Based on this fact, and to be able to compare to \cite{wagnergrabarczykhong2019}, we will continue the analysis using model \eqref{eq:ekc_quadraticmodel} and the given collection of countries. Before doing so, it will be worthwhile to discuss the time series properties of these residuals.
We consider the series $\{\hat \vu_{t,FGLS}\}$ in the remainder of this section but the other residuals series will provide qualitatively similar outcomes. When fitting the VAR($p$) models with $1\leq p \leq 8$ to these residuals, the BIC information criterion selects a lag order of $p=1$. The absolute eigenvalues of the estimated coefficient matrix are $(0.55,0.55,0.51,0.31,0.31,0.11)$, and the estimate for the error correlation matrix is
$$
\begin{blockarray}{c c c c c c c}
& AT & BE & FI & NL & CH & UK \\
\begin{block}{c [c c c c c c]}
AT & 1 & 0.16 & 0.10 & 0.16 & 0.18 & 0.22\\
BE & \bullet & 1 & 0.51 & 0.09 & 0.23 & 0.27\\
FI & \bullet & \bullet & 1 & 0.13 & 0.26 & 0.22\\
NL & \bullet & \bullet & \bullet & 1 & 0.22 & 0.10 \\
CH & \bullet & \bullet & \bullet & \bullet & 1 & 0.18\\
UK & \bullet & \bullet & \bullet & \bullet & \bullet & 1 \\
\end{block}
\end{blockarray}.
$$
There is thus serial and cross-sectional correlation to be exploited by the FM-GLS estimator.
The FM-SOLS, FM-SUR and FM-GLS estimation results of Model \eqref{eq:ekc_quadraticmodel} are reported in Table \ref{tab:estimationresults}. An inspection of the coefficient estimates and their confidence intervals reveals that: (1) $\beta_{i,3}$ is positive for each country, (2) $\beta_{i,4}$ is negative for each country, and (3) all coefficients are significant at the 5\% level. All these three facts are in line with the EKC hypothesis.\footnote{This is non-surprising because \cite{wagnergrabarczykhong2019} have selected the current set of countries because they display the EKC behaviour. Also, our estimation results are slightly different from those in \cite{wagnergrabarczykhong2019} due to the additional data for 2014, possible data updates, and/or differences in the bandwidth selection of the long-run covariance matrices.} Accordingly, there exists a turning point after which further per capita economic growth reduces per capita carbon dioxide emissions. The numerical values for the turning points are heterogeneous between countries.
The widths of the confidence intervals for $\beta_{i,3}$ and $\beta_{i,4}$ display a similar pattern. From shortest to longest, the ordering is always FM-SUR, FM-SOLS, and FM-GLS and we also see how widths vary substantially between methods. To uncover the origin of these findings we conduct one final simulation study with a parameter specification that closely mimics the properties of the dataset.\footnote{The details of this simulation DGP are provided in Section \ref{empiricalillustration} of the Supplement. A visualisation of the data and the model fit are also provided there.} The average empirical coverage probabilities of asymptotic 95\% confidence intervals are 78.0\%, 66.5\% and 89.0\% for FM-SOLS, FM-SUR, and FM-GLS, respectively. In other words, the calculated confidence intervals are generally too short. By reverse engineering it turns out that the confidence intervals should be scaled by factors of 1.67, 2.15 and 1.24 to bring them back to the desired nominal level. Overall, the applied researcher should be careful when using the confidence intervals as indications for parameter uncertainty.
\section{Conclusion} \label{sec:conclusion}
We proposed a framework to conduct inference on cointegrating polynomial regressions. Parameters are obtained using a Fully Modified GLS estimator and we studied a cointegration test that is based on filtered residuals. Monte Carlo simulations revealed the advantages and disadvantages of these methods. The empirical researcher should realize that all estimation approaches have a tendency to underestimate parameter uncertainty and thus provide confidence intervals that are too small. The FM-GLS estimator suffers the least from this problem. Several interesting questions are left for future research. From a theoretical viewpoint, it is interesting to study the behaviour of the modified Cholesky decomposition (and BIAM) when the series under consideration is nonstationary. This would give insights into the behaviour of: (1) the FM-GLS estimator while estimating spurious regressions, and (2) the power properties of the cointegration tests. From a practical viewpoint, there seems a need to obtain more acurate standard errors of the parameter estimators.
\section*{Acknowledgements}
This paper has been presented at the 2018 CFE meeting in Pisa, the NESG 2019 conference in Amsterdam, and the $6^{th}$ RCEA Time Series Econometrics Workshop in Larnaca. We would like to thank conference participants, especially Peter Pedroni and Peter Phillips, for useful comments and suggestions. We extend our thanks to Eric Beutner, Dick van Dijk, Richard Paap, Franz Palm, and Stephan Smeekes for their valuable feedback on earlier versions of this manuscript. All remaining errors are our own.
\clearpage
\bibliographystyle{chicagoa}
\newpage
\begin{thebibliography}{}
\bibitem[\protect\citeauthoryear{Abadir and Magnus}{Abadir and
Magnus}{2005}]{abadirmagnus2005}
Abadir, K.~M. and J.~R. Magnus (2005).
\newblock {\em Matrix Algebra}.
\newblock Cambridge University Press.
\bibitem[\protect\citeauthoryear{Anderson and Darling}{Anderson and
Darling}{1952}]{andersondarling1952}
Anderson, T.~W. and D.~A. Darling (1952).
\newblock Asymptotic theory of certain ``goodness of fit'' criteria based on
stochastic processes.
\newblock {\em The Annals of Mathematical Statistics\/}~{\em 23}, 193--212.
\bibitem[\protect\citeauthoryear{Andrews}{Andrews}{1991}]{andrews1991}
Andrews, D. W.~K. (1991).
\newblock Heteroskedasticity and autocorrelation consistent covariance matrix
estimation.
\newblock {\em Econometrica\/}~{\em 59}, 817--858.
\bibitem[\protect\citeauthoryear{Berk}{Berk}{1974}]{berk1974}
Berk, K.~N. (1974).
\newblock Consistent autoregressive spectral estimates.
\newblock {\em The Annals of Statistics\/}~{\em 2}, 489--502.
\bibitem[\protect\citeauthoryear{Beutner, Lin, and Smeekes}{Beutner
et~al.}{2019}]{beutnerlinsmeekes2019}
Beutner, E., Y.~Lin, and S.~Smeekes (2019).
\newblock {GLS} estimation and confidence sets for the date of a single break
in models with trends.
\newblock Working Paper.
\bibitem[\protect\citeauthoryear{Bickel and Levina}{Bickel and
Levina}{2008}]{bickellevina2008}
Bickel, P.~J. and E.~Levina (2008).
\newblock Regularized estimation of large covariance matrices.
\newblock {\em The Annals of Statistics\/}~{\em 36}, 199--227.
\bibitem[\protect\citeauthoryear{Boden, Marland, and Andres}{Boden
et~al.}{2017}]{cdiac2017}
Boden, T., G.~Marland, and R.~Andres (2017).
\newblock Global, regional, and national fossil-fuel $\text{CO}_2$ emissions.
\newblock Carbon Dioxide Information Analysis Center, Oak Ridge National
Laboratory, U.S. Department of Energy, Oak Ridge, Tenn., U.S.A.
\url{http://cdiac.ess-dive.lbl.gov/trends/emis/tre_coun.html}.
\bibitem[\protect\citeauthoryear{Bolt, Inklaar, de~Jong, and van Zanden}{Bolt
et~al.}{2018}]{madison2018}
Bolt, J., R.~Inklaar, H.~de~Jong, and J.~L. van Zanden (2018).
\newblock Rebasing "maddison": New income comparisons and the shape of long-run
economic development.
\newblock
\url{https://www.rug.nl/ggdc/historicaldevelopment/maddison/releases/maddison-project-database-2018}.
\bibitem[\protect\citeauthoryear{Breitung}{Breitung}{2001}]{breitung2001}
Breitung, J. (2001).
\newblock Rank tests for nonlinear cointegration.
\newblock {\em Journal of Business \& Economic Statistics\/}~{\em 19},
331--340.
\bibitem[\protect\citeauthoryear{Chang, Park, and Phillips}{Chang
et~al.}{2004}]{chang2004}
Chang, Y., J.~Y. Park, and P.~C.~B. Phillips (2004).
\newblock Bootstrap unit root tests in panels with cross-sectional dependency.
\newblock {\em Journal of Econometrics\/}~{\em 120}, 263--293.
\bibitem[\protect\citeauthoryear{Cheng and Pourahmadi}{Cheng and
Pourahmadi}{1993}]{chengpourahmadi1993}
Cheng, R. and M.~Pourahmadi (1993).
\newblock Baxter's inequality and convergence of finite predictors of
multivariate stochastic processes.
\newblock {\em Probability Theory and Related Fields\/}~{\em 95}, 115--124.
\bibitem[\protect\citeauthoryear{Cheng, Ing, and Yu}{Cheng
et~al.}{2015}]{chengingyu2015}
Cheng, T.-C.~F., C.-K. Ing, and S.-H. Yu (2015).
\newblock Toward optimal model averaging in regression models with time series
errors.
\newblock {\em Journal of Econometrics\/}~{\em 189}, 321--334.
\bibitem[\protect\citeauthoryear{Choi and Saikkonen}{Choi and
Saikkonen}{2010}]{choisaikkonen2010}
Choi, I. and P.~Saikkonen (2010).
\newblock Tests for nonlinear cointegration.
\newblock {\em Econometric Theory\/}~{\em 26}, 682--709.
\bibitem[\protect\citeauthoryear{Davidson}{Davidson}{1994}]{davidson1994}
Davidson, J. (1994).
\newblock {\em Stochastic Limit Theory}.
\newblock Oxford University Press.
\bibitem[\protect\citeauthoryear{Davidson and MacKinnon}{Davidson and
MacKinnon}{2004}]{davidsonmackinnon2004}
Davidson, R. and J.~G. MacKinnon (2004).
\newblock {\em Econometric Theory and Methods}.
\newblock Oxford University Press.
\bibitem[\protect\citeauthoryear{de~Jong}{de~Jong}{2002}]{dejong2002}
de~Jong, R.~M. (2002).
\newblock Nonlinear estimators with integrated regressors but without
exogeneity.
\newblock mimeo Michigan State University.
\bibitem[\protect\citeauthoryear{Findley and Wei}{Findley and
Wei}{1993}]{findleywei1993}
Findley, D.~F. and C.-Z. Wei (1993).
\newblock Moment bounds for deriving time series {CLT}'s and model selection
procedures.
\newblock {\em Statistica Sinica\/}~{\em 3}, 453--480.
\bibitem[\protect\citeauthoryear{Grossman and Krueger}{Grossman and
Krueger}{1995}]{grossmankrueger1995}
Grossman, G.~M. and A.~B. Krueger (1995).
\newblock Economic growth and the environment.
\newblock {\em The Quarterly Journal of Economics\/}~{\em 110}, 353--377.
\bibitem[\protect\citeauthoryear{Hamilton}{Hamilton}{1994}]{hamilton1994}
Hamilton, J.~D. (1994).
\newblock {\em Time Series Analysis}.
\newblock Princeton University Press.
\bibitem[\protect\citeauthoryear{Hannan and Deistler}{Hannan and
Deistler}{2012}]{hannandeistler2012}
Hannan, E. and M.~Deistler (2012).
\newblock {\em The Statistical Theory of Linear Systems}.
\newblock Society for Industrial and Applied Mathematics.
\bibitem[\protect\citeauthoryear{Hong and Phillips}{Hong and
Phillips}{2010}]{hongphillips2010}
Hong, S.~H. and P.~C.~B. Phillips (2010).
\newblock Testing linearity in cointegrating relations with an application to
purchasing power parity.
\newblock {\em Journal of Business \& Economic Statistics\/}~{\em 28}, 96--114.
\bibitem[\protect\citeauthoryear{Ing, Chiou, and Guo}{Ing
et~al.}{2016a}]{ICG2016}
Ing, C.-K., H.-T. Chiou, and M.~Guo (2016a).
\newblock Estimation of inverse autocovariance matrices for long memory
processes.
\newblock {\em Bernoulli\/}~{\em 22}, 1301--1330.
\bibitem[\protect\citeauthoryear{Ing, Chiou, and Guo}{Ing
et~al.}{2016b}]{ingchiouguo2016}
Ing, C.-K., H.-T. Chiou, and M.~Guo (2016b).
\newblock Estimation of inverse autocovariance matrices for long memory
processes.
\newblock {\em Bernoulli\/}~{\em 22}, 1301--1330.
\bibitem[\protect\citeauthoryear{Jansson}{Jansson}{2002}]{jansson2002}
Jansson, M. (2002).
\newblock Consistent covariance matrix estimation for linear processes.
\newblock {\em Econometric Theory\/}~{\em 18}, 1449--1459.
\bibitem[\protect\citeauthoryear{Kim and Zimmerman}{Kim and
Zimmerman}{2012}]{kimzimmerman2012}
Kim, C. and D.~L. Zimmerman (2012).
\newblock Unconstrained models for the covariance structure of multivariate
longitudinal data.
\newblock {\em Journal of Multivariate Analysis\/}~{\em 107}, 104--118.
\bibitem[\protect\citeauthoryear{Kohli, Garcia, and Pourahmadi}{Kohli
et~al.}{2016}]{kohligarciapourahmadi2016}
Kohli, P., T.~P. Garcia, and M.~Pourahmadi (2016).
\newblock Modeling the cholesky factors of covariance matrices of multivariate
longitudinal data.
\newblock {\em Journal of Multivariate Analysis\/}~{\em 145}, 87--100.
\bibitem[\protect\citeauthoryear{Lewis and Reinsel}{Lewis and
Reinsel}{1985}]{lewisreinsel1985}
Lewis, R. and G.~C. Reinsel (1985).
\newblock Prediction of multivariate time series by autoregressive model
fitting.
\newblock {\em Journal of Multivariate Analysis\/}~{\em 16}, 393--411.
\bibitem[\protect\citeauthoryear{Li, Phillips, and Gao}{Li
et~al.}{2020}]{lihillipsgao2017}
Li, D., P.~C.~B. Phillips, and J.~Gao (2020).
\newblock Kernel-based inference in time-varying coefficient cointegrating
regression.
\newblock {\em Journal of Econometrics\/}~{\em 215}, 607--632.
\bibitem[\protect\citeauthoryear{Mark, Ogaki, and Sul}{Mark
et~al.}{2005}]{markogakisul2005}
Mark, N.~C., M.~Ogaki, and D.~Sul (2005).
\newblock Dynamic seemingly unrelated cointegrating regressions.
\newblock {\em The Review of Economic Studies\/}~{\em 72}, 797--820.
\bibitem[\protect\citeauthoryear{McMurry and Politis}{McMurry and
Politis}{2010}]{mcmurypolitis2010}
McMurry, T.~L. and D.~N. Politis (2010).
\newblock Banded and tapered estimates for autocovariance matrices and the
linear process bootstrap.
\newblock {\em Journal of Time Series Analysis\/}~{\em 31}, 471--482.
\bibitem[\protect\citeauthoryear{Moon and Perron}{Moon and
Perron}{2005}]{moonperron2005}
Moon, H.~R. and B.~Perron (2005).
\newblock Efficient estimation of the seemingly unrelated regression
cointegration model and testing for purchasing power parity.
\newblock {\em Econometric Reviews\/}~{\em 23}, 293--323.
\bibitem[\protect\citeauthoryear{Newey and West}{Newey and
West}{1994}]{neweywest1994}
Newey, W.~K. and K.~D. West (1994).
\newblock Automatic lag selection in covariance matrix estimation.
\newblock {\em The Review of Economic Studies\/}~{\em 61}, 631--653.
\bibitem[\protect\citeauthoryear{Nordhaus}{Nordhaus}{2013}]{nordhaus2013}
Nordhaus, W.~D. (2013).
\newblock {\em The Climate Casino: Risk, Uncertainty, and Economics for a
Warming World}.
\newblock Yale University Press.
\bibitem[\protect\citeauthoryear{Nyblom and Harvey}{Nyblom and
Harvey}{2000}]{nyblomharvey2000}
Nyblom, J. and A.~Harvey (2000).
\newblock Tests of common stochastic trends.
\newblock {\em Econometric Theory\/}~{\em 16}, 176--199.
\bibitem[\protect\citeauthoryear{Park and Phillips}{Park and
Phillips}{1999}]{parkphillips1999}
Park, J.~Y. and P.~C.~B. Phillips (1999).
\newblock Asymptotics for nonlinear transformations of integrated time series.
\newblock {\em Econometric Theory\/}~{\em 15\/}(3), 269--298.
\bibitem[\protect\citeauthoryear{Park and Phillips}{Park and
Phillips}{2001}]{parkphillips2001}
Park, J.~Y. and P.~C.~B. Phillips (2001).
\newblock Nonlinear regressions with integrated time series.
\newblock {\em Econometrica\/}~{\em 69}, 117--161.
\bibitem[\protect\citeauthoryear{Phillips}{Phillips}{1991}]{phillips1991}
Phillips, P. C.~B. (1991).
\newblock Optimal inference in cointegrated systems.
\newblock {\em Econometrica\/}~{\em 59}, 283--306.
\bibitem[\protect\citeauthoryear{Phillips}{Phillips}{1995}]{phillips1995}
Phillips, P. C.~B. (1995).
\newblock Fully modified least squares and vector autoregression.
\newblock {\em Econometrica\/}~{\em 63}, 1023--1078.
\bibitem[\protect\citeauthoryear{Phillips and Hansen}{Phillips and
Hansen}{1990}]{phillipshansen1990}
Phillips, P. C.~B. and B.~E. Hansen (1990).
\newblock Statistical inference in instrumental variables regression with
{I}(1) processes.
\newblock {\em The Review of Economic Studies\/}~{\em 57}, 99--125.
\bibitem[\protect\citeauthoryear{Phillips and Park}{Phillips and
Park}{1988}]{phillipspark1988}
Phillips, P. C.~B. and J.~Y. Park (1988).
\newblock Asymptotic equivalence of ordinary least squares and generalized
least squares in regressions with integrated regressors.
\newblock {\em Journal of the American Statistical Association\/}~{\em 83},
111--115.
\bibitem[\protect\citeauthoryear{Phillips and Solo}{Phillips and
Solo}{1992}]{phillipssolo1992}
Phillips, P. C.~B. and V.~Solo (1992).
\newblock Asymptotics for linear processes.
\newblock {\em The Annals of Statistics\/}~{\em 20}, 971--1001.
\bibitem[\protect\citeauthoryear{Pourahmadi}{Pourahmadi}{1999}]{pourahmadi1999}
Pourahmadi, M. (1999).
\newblock Joint mean-covariance models with applications to longitudinal data:
Unconstrained parameterisation.
\newblock {\em Biometrika\/}~{\em 86}, 677--690.
\bibitem[\protect\citeauthoryear{Prais and Winsten}{Prais and
Winsten}{1954}]{praiswinston1954}
Prais, S.~J. and C.~B. Winsten (1954).
\newblock Trend estimators and serial correlation.
\newblock Cowles Foundation, Discussion Paper 383.
\bibitem[\protect\citeauthoryear{Romano and Wolf}{Romano and
Wolf}{2001}]{romanowolf2001}
Romano, J.~P. and M.~Wolf (2001).
\newblock Subsampling intervals in autoregressive models with linear time
trend.
\newblock {\em Econometrica\/}~{\em 69}, 1283--1314.
\bibitem[\protect\citeauthoryear{Saikkonen}{Saikkonen}{1992}]{saikkonen1992}
Saikkonen, P. (1992).
\newblock Estimation and testing of cointegrated systems by an autoregressive
approximation.
\newblock {\em Econometric Theory\/}~{\em 8}, 1--27.
\bibitem[\protect\citeauthoryear{Shin}{Shin}{1994}]{shin1994}
Shin, Y. (1994).
\newblock A residual-based test of the null of cointegration against the
alternative of no cointegration.
\newblock {\em Econometric Theory\/}~{\em 10}, 91--115.
\bibitem[\protect\citeauthoryear{Sims, Stock, and Watson}{Sims
et~al.}{1990}]{simsstockwatson1990}
Sims, C.~A., J.~H. Stock, and M.~W. Watson (1990).
\newblock Inference in linear time series mmodels with some unit roots.
\newblock {\em Econometrica\/}~{\em 58}, 113--144.
\bibitem[\protect\citeauthoryear{Stern}{Stern}{2004}]{stern2004}
Stern, D.~I. (2004).
\newblock The rise and fall of the environmental {K}uznets curve.
\newblock {\em World Development\/}~{\em 32}, 1419--1439.
\bibitem[\protect\citeauthoryear{Stern}{Stern}{2017}]{stern2017}
Stern, D.~I. (2017).
\newblock The environmental {K}uznets curve after 25 years.
\newblock {\em Journal of Bioeconomics\/}~{\em 19}, 7--28.
\bibitem[\protect\citeauthoryear{Tao}{Tao}{2012}]{tao2012}
Tao, T. (2012).
\newblock {\em Topics in Random Matrix Theory}.
\newblock American Mathematical Society.
\bibitem[\protect\citeauthoryear{Tj\o{}stheim}{Tj\o{}stheim}{2020}]{tjostheim2020}
Tj\o{}stheim, D. (2020).
\newblock Some notes on nonlinear cointegration: A partial review with some
novel perspectives.
\newblock {\em Econometric Reviews\/}~{\em 39}, 655--673.
\bibitem[\protect\citeauthoryear{Vogelsang and Wagner}{Vogelsang and
Wagner}{2014}]{vogelsangwagner2014}
Vogelsang, T.~J. and M.~Wagner (2014).
\newblock Integrated modified {OLS} estimation and fixed-b inference for
cointegrating regressions.
\newblock {\em Journal of Econometrics\/}~{\em 178}, 741--760.
\bibitem[\protect\citeauthoryear{Wagner, Grabarczyk, and Hong}{Wagner
et~al.}{2020}]{wagnergrabarczykhong2019}
Wagner, M., P.~Grabarczyk, and S.~H. Hong (2020).
\newblock Fully modified {OLS} estimation and inference for seemingly unrelated
cointegrating polynomial regressions and the environmental {K}uznets curve
for carbon dioxide emissions.
\newblock {\em Journal of Econometrics\/}~{\em 214}, 216--255.
\bibitem[\protect\citeauthoryear{Wagner and Hong}{Wagner and
Hong}{2016}]{wagnerhong2016}
Wagner, M. and S.~H. Hong (2016).
\newblock Cointegrating polynomial regressions: Fully modified {OLS} estimation
and inference.
\newblock {\em Econometric Theory\/}~{\em 32}, 1289--1315.
\bibitem[\protect\citeauthoryear{Wang}{Wang}{2015}]{wang2015}
Wang, Q. (2015).
\newblock {\em Limit Theorems for Nonlinear Cointegrating Regression}.
\newblock World Scientific.
\bibitem[\protect\citeauthoryear{Wang and Phillips}{Wang and
Phillips}{2009}]{wangphillips2009}
Wang, Q. and P.~C.~B. Phillips (2009).
\newblock Structural nonparametric cointegrating regression.
\newblock {\em Econometrica\/}~{\em 77}, 1901--1948.
\bibitem[\protect\citeauthoryear{Wei}{Wei}{1987}]{wei1987}
Wei, C.-Z. (1987).
\newblock Adaptive prediction by least squares predictors in stochastic
regression models with applications to time series.
\newblock {\em The Annals of Statistics\/}~{\em 15}, 1667--1682.
\bibitem[\protect\citeauthoryear{Wu and Pourahmadi}{Wu and
Pourahmadi}{2009}]{wupourahmadi2009}
Wu, W.~B. and M.~Pourahmadi (2009).
\newblock Banding sample autocovariance matrices of stationary processes.
\newblock {\em Statistica Sinica\/}~{\em 19}, 1755--1768.
\bibitem[\protect\citeauthoryear{Zellner}{Zellner}{1962}]{zellner1962}
Zellner, A. (1962).
\newblock An efficient method of estimating seemingly unrelated regressions and
tests for aggregation bias.
\newblock {\em Journal of the American Statistical Association\/}~{\em 57},
348--368.
\end{thebibliography}
\clearpage