Efficient Estimation by Fully Modified GLS with an Application to the Environmental Kuznets Curve

Efficient Estimation by Fully Modified GLS with an Application to the Environmental Kuznets Curve

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 wang2015 , the recent review article by 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 breitung2001, or nonlinear purchasing power parity due to transaction/transportation costs and trade barriers (e.g. 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 grossmankrueger1995.\footnote{There is no direct reference to Kuznets in the original paper by 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 parkphillips1999 and 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 wangphillips2009 or 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. wagnerhong2016 named this a cointegrating polynomial regression (CPR). The multivariate seemingly unrelated regressions extension is available in 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 markogakisul2005 and 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 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 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) 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) whereas a finite sample simulation study is presented in Section (ref). The empirical application can be found in Section (ref) where we look at the environmental Kuznets curve. Section (ref) 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)$.

The Model

As in wagnergrabarczykhong2019, we study a system of seemingly unrelated cointegrating polynomial regressions (SUCPR), that is

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 (ref) 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

(see Assumption (ref) for further details). Efficient estimation of the parameter vector $\vbeta$ now requires the use of generalized least squares (GLS). Our 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 pourahmadi1999. This extension was named the Modified Cholesky Block Decomposition (MCBD) by kimzimmerman2012 and kohligarciapourahmadi2016. The latter papers used the MCBD to parametrize the covariance matrix of multivariate longitudinal data. As in 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

and $\mS(0)=\E(\vu_t^{}\vu_t')$. The inverse of the covariance matrix $\mSigma_{\vu}$ is then given by

where $ \bm{\mathcal{S}}_{\vu} = \diag \Big(\mS(0),\mS(1),\ldots,\mS(T-1) \Big)$,

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 Banded Inverse Autocovariance Matrix (BIAM) as

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

The model of (ref) 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, we will assume $\mSigma_{\vu}^{-1}(q)$ to be known and focus on the following estimator:

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 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.

Asymptotic Theory

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 wagnergrabarczykhong2019. The following assumption will facilitate the development of the asymptotic theory.

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) ensure that the lag polynomials $\bm{\mathcal{A}}(L)$ and $\bm{\mathcal{D}}(L)$ are invertible (see for example Theorem 7.4.2 of hannandeistler2012), thereby showing that our Assumption (ref) is similar to the linear processes assumptions that are regularly adopted in the literature on nonlinear cointegration, cf. choisaikkonen2010, wagnerhong2016, and wagnergrabarczykhong2019. The assumption $\det\big(\bm{\mathcal{D}}(1)\big)\neq 0$ rules out cointegration among the components of $\{\vx_t\}$.

Under Assumption (ref)(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[

\right]$ where $\bm{B}_{\vzeta}$ denotes an $2n$-dimensional Brownian motion with covariance matrix $\left[

\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

where the Brownian motion $\bm{B}_{\vxi}(r)$ of dimension $2n$ has covariance matrix

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[

\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}^$.

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 (ref) can now be used to derive the limiting distribution of this infeasible GLS estimator. A insightful exposition of our results requires further notation.

Finally, we use $\bm{B}_{v_j}$ as shorthand notation for the $j^{th}$ component of $\bm{B}_v$.

The limiting result in (ref) 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 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. 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 (ref) reflect the presence of second order bias terms caused by serial correlation and endogeneity. In Section (ref), 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.

Consistent Estimation of $\mSigma_{\vu}^{-1}(q)$ and Feasible GLS

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. wupourahmadi2009 and 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 (ref). 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 chengingyu2015 and ICG2016 for univariate time series. For a multivariate time series, we define

$1\leq \ell\leq q$, and $\widehat{\mS}(0)=\frac{1}{T}\sum_{t=1}^{T}\widehat{\vu}_t^{}\widehat{\vu}_t'$. Similarly to (ref)-(ref), 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

Assumption (ref) 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) places constraints on the banding parameter $q_T$. First, Assumption (ref) 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 (ref), shows that we are fitting a vector autoregression (VAR) of increasing lag order to the residuals. Identical rate requirements are reported by 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.

Fully Modified Inference

The asymptotic results of Theorem (ref) 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 phillips1991, saikkonen1992's\ (saikkonen1992) dynamic least squares, and the integrated modified OLS and fixed-b approaches by vogelsangwagner2014. We rely on the fully modified (FM) approach advocated by phillipshansen1990 and 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, 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)). As such, we introduce the following assumption.

andrews1991 and 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) is easily satisfied by imposing suitable conditions on the kernel function and bandwidth parameter. We refer to phillips1995 and 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).\footnote{An overview of the procedure is given here. Section (ref) in the Supplement provides further details.} This approach resembles 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 (ref)). 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

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.

Definitions and limiting results for FM estimators are presented in Theorem (ref). 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.

The FM-GLS estimator is new to the seemingly unrelated CPR literature, whereas the FM-SOLS and FM-SUR estimators have recently appeared in wagnergrabarczykhong2019. Theorem (ref) 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) 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 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)). 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.

Testing the Null of Cointegration

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 nyblomharvey2000 with the sub-sampling approach found in choisaikkonen2010 and wagnerhong2016. We consider three test statistics. To treat all test statistics in a unified framework, we define

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. 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 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:

and

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 (ref) fits naturally into the FM-GLS estimation framework.

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 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). 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$.

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 wagnergrabarczykhong2019. All long-run covariance matrices are computed using a Bartlett kernel and the automatic bandwidth selection approach due to andrews1991. For FM-GLS, the banding parameter $q_T$ is selected using the subsampling and risk-minimization approach explained in section 5 from 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 (ref), whereas $K^{BIAM}$ employs the pre-filtered residuals from (ref). 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.

Monte Carlo Designs

We generate data according to a quadratic seemingly unrelated CPR. That is, we adopt the DGP in (ref) 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.

\noindentSetting A (Errors as in wagnergrabarczykhong2019): As a benchmark, we revisit the simulation setting in wagnergrabarczykhong2019 and generate innovations according to

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

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$.

\noindentSetting B (VARMA Errors): To further investigate the importance of serial correlation, we consider a second specification of the innovation process:

where $\veta_t$ and $\vepsi_{t}$ are generated as $\left[

\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 chang2004. That is, we take the following three steps:

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.

\noindentSetting C (Cointegration Tests): We continue to construct innovations according to Setting B. Moreover, we fix $\left[

\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(\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 quadratic seemingly unrelated CPR.

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.

Discussion of the Simulation Results

Tables (ref) and (ref) 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:

The subsequent set of simulations evolves around hypothesis testing, see Table (ref) and Figures (ref)-(ref). 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 wagnergrabarczykhong2019, and the Wald-FGLS test from Theorem (ref). We consider: (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 (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.

The simulation results for the KPSS-type of cointegration tests can be found in Table (ref). The general conclusions are as follows.

Empirical Application

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 stern2004 and 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 (nordhaus2013).

We builds upon and compare with 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, madison2018, contains the data on population size and real GDP. The data on $\text{CO}_2$ originates from 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 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:

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) to check this model specification (Table (ref)). All three tests reject the null of cointegration at a 5% level signalling inappropriateness of the quadratic formulation. Figure (ref) 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 wagnergrabarczykhong2019, we will continue the analysis using model (ref) 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 $$

. $$ 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 (ref) are reported in Table (ref). 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 wagnergrabarczykhong2019 have selected the current set of countries because they display the EKC behaviour. Also, our estimation results are slightly different from those in 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) 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.

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.

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.